point_defect_diffusion calculation style
Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.
Description updated: 2019-07-26
Introduction
The point_defect_diffusion calculation style estimates the diffusion rate of a point defect at a specified temperature. A system is created with a single point defect, and subjected to a long time molecular dynamics simulation. The mean square displacement for the defect is computed, and used to estimate the diffusion constant.
Disclaimers
The calculation estimates the defect’s diffusion by computing the mean square displacement of all atoms in the system. This is useful for estimating rates associated with vacancies and self-interstitials as the process is not confined to a single atom’s motion. However, this makes the calculation ill-suited to measuring diffusion of substitutional impurities as it does not individually track each atom’s position throughout.
Method and Theory
First, a defect system is constructed by adding a single point defect (or defect cluster) to an initially bulk system using the atomman.defect.point() function.
A LAMMPS simulation is then performed on the defect system. The simulation consists of two separate runs
NVT equilibrium run: The system is allowed to equilibrate at the given temperature using nvt integration.
NVE measurement run: The system is then evolved using nve integration, and the total mean square displacement of all atoms is measured as a function of time.
Between the two runs, the atomic velocities are scaled such that the average temperature of the nve run is closer to the target temperature.
The mean square displacement of the defect, \(\left< \Delta r_{ptd}^2 \right>\) is then estimated using the mean square displacement of the atoms \(\left< \Delta r_{i}^2 \right>\). Under the assumption that all diffusion is associated with the single point defect, the defect’s mean square displacement can be taken as the summed square displacement of the atoms
where \(N\) is the number of atoms in the system. The diffusion constant is then estimated by linearly fitting the change in mean square displacement with time
where d is the number of dimensions included.
Demonstration
1. Setup
1.1. Library imports
Import libraries needed by the calculation. The external libraries used are:
[1]:
# Standard library imports
from pathlib import Path
import os
import sys
import uuid
import shutil
import random
import datetime
from copy import deepcopy
# http://www.numpy.org/
import numpy as np
# https://github.com/usnistgov/DataModelDict
from DataModelDict import DataModelDict as DM
# https://github.com/usnistgov/atomman
import atomman as am
import atomman.lammps as lmp
import atomman.unitconvert as uc
# https://github.com/usnistgov/iprPy
import iprPy
print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)
Notebook last executed on 2019-07-29 using iprPy version 0.9.0
1.2. Default calculation setup
[2]:
# Specify calculation style
calc_style = 'point_defect_diffusion'
# If workingdir is already set, then do nothing (already in correct folder)
try:
workingdir = workingdir
# Change to workingdir if not already there
except:
workingdir = Path('calculationfiles', calc_style)
if not workingdir.is_dir():
workingdir.mkdir(parents=True)
os.chdir(workingdir)
2. Assign values for the calculation’s run parameters
2.1. Specify system-specific paths
lammps_command (required) is the LAMMPS command to use.
mpi_command MPI command for running LAMMPS in parallel. A value of None will run simulations serially.
[3]:
lammps_command = 'lmp_serial'
mpi_command = None
2.2. Load interatomic potential
potential_name gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.
potential_file gives the path to the potential_LAMMPS reference record to use. Here, this parameter is automatically generated using potential_name and librarydir.
potential_dir gives the path for the folder containing the artifacts associated with the potential (i.e. eam.alloy file). Here, this parameter is automatically generated using potential_name and librarydir.
potential is an atomman.lammps.Potential object (required). Here, this parameter is automatically generated from potential_file and potential_dir.
[4]:
potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'
# Define potential_file and potential_dir using librarydir and potential_name
potential_file = Path(iprPy.libdir, 'potential_LAMMPS', f'{potential_name}.json')
potential_dir = Path(iprPy.libdir, 'potential_LAMMPS', potential_name)
# Initialize Potential object using potential_file and potential_dir.
potential = lmp.Potential(potential_file, potential_dir)
print('Successfully loaded potential', potential)
Successfully loaded potential 1999--Mishin-Y--Ni--LAMMPS--ipr1
2.3. Load initial unit cell system
prototype_name gives the name of the crystal_prototype reference record in the iprPy library to load.
symbols is a list of the potential’s elemental model symbols to associate with the unique atom types of the loaded system.
box_parameters is a list of the a, b, c lattice constants to assign to the loaded file.
load_file gives the path to the atomic configuration file to load for the ucell system. Here, this is generated automatically using prototype_name and librarydir.
load_style specifies the format of load_file. Here, this is automatically set for crystal_prototype records.
load_options specifies any other keyword options for properly loading the load_file. Here, this is automatically set for crystal_prototype records.
ucell is an atomman.System representing a fundamental unit cell of the system (required). Here, this is generated using the load parameters and symbols.
[5]:
prototype_name = 'A1--Cu--fcc'
symbols = ['Ni']
box_parameters = uc.set_in_units([3.52, 3.52, 3.52], 'angstrom')
# Define load_file using librarydir and prototype_name
load_file = Path(iprPy.libdir, 'crystal_prototype', f'{prototype_name}.json')
# Define load_style and load_options for crystal_prototype records
load_style = 'system_model'
load_options = {}
# Create ucell by loading prototype record
ucell = am.load(load_style, load_file, symbols=symbols, **load_options)
# Rescale ucell using box_parameters
ucell.box_set(a=box_parameters[0], b=box_parameters[1], c=box_parameters[2], scale=True)
print(ucell)
avect = [ 3.520, 0.000, 0.000]
bvect = [ 0.000, 3.520, 0.000]
cvect = [ 0.000, 0.000, 3.520]
origin = [ 0.000, 0.000, 0.000]
natoms = 4
natypes = 1
symbols = ('Ni',)
pbc = [ True True True]
per-atom properties = ['atype', 'pos']
id atype pos[0] pos[1] pos[2]
0 1 0.000 0.000 0.000
1 1 0.000 1.760 1.760
2 1 1.760 0.000 1.760
3 1 1.760 1.760 0.000
2.4. Specify the defect parameters
pointdefect_name gives the name of a point-defect reference record in the iprPy library containing point defect input parameters.
pointdefect_file gives the path to a point_defect reference containing point defect input parameters. Here, this is built automatically using pointdefect_name and librarydir.
point_kwargs (required) is a dictionary or list of dictonaries containing parameters for generating the defect. Here, values are extracted from pointdefect_file. Allowed keywords are:
ptd_type indicates which defect type to generate: ‘v’ for vacancy, ‘i’ for interstitial, ‘s’ for substitutional, or ‘db’ for dumbbell.
atype is the atom type to assign to the defect atom (‘i’, ‘s’, ‘db’ ptd_types).
pos specifies the position for adding the defect atom (all ptd_types).
ptd_id specifies the id of an atom in the initial system where the defect is to be added. Alternative to using pos (‘v’, ‘s’, ‘db’ ptd_types).
db_vect gives the vector associated with the dumbbell interstitial to generate (‘db’ ptd_type).
scale indicates if pos and db_vect are in absolute (False) or box-relative (True) coordinates. Default is False.
atol is the absolute tolerance for position-based searching. Default is 1e-3 angstroms.
[6]:
pointdefect_name = 'A1--Cu--fcc--vacancy'
#pointdefect_name = 'A1--Cu--fcc--1nn-divacancy'
#pointdefect_name = 'A1--Cu--fcc--2nn-divacancy'
#pointdefect_name = 'A1--Cu--fcc--100-dumbbell'
#pointdefect_name = 'A1--Cu--fcc--110-dumbbell'
#pointdefect_name = 'A1--Cu--fcc--111-dumbbell'
#pointdefect_name = 'A1--Cu--fcc--octahedral-interstitial'
#pointdefect_name = 'A1--Cu--fcc--tetrahedral-interstitial'
#pointdefect_name = 'A1--Cu--fcc--crowdion-interstitial'
# Define pointdefect_file using librarydir and pointdefect_name
pointdefect_file = Path(iprPy.libdir, 'point_defect', f'{pointdefect_name}.json')
# Parse pointdefect_file using iprPy.input.interpret()
defectinputs = {'ucell':ucell, 'pointdefect_file':pointdefect_file}
iprPy.input.subset('pointdefect').interpret(defectinputs)
# Extract point_kwargs
point_kwargs = defectinputs['point_kwargs']
print('point_kwargs =')
point_kwargs
point_kwargs =
[6]:
[{'ptd_type': 'v', 'pos': array([0., 0., 0.]), 'scale': False}]
2.5. Modify system
sizemults list of three integers specifying how many times the ucell vectors of \(a\), \(b\) and \(c\) are replicated in creating system.
system is an atomman.System to perform the scan on (required).
[7]:
sizemults = [10, 10, 10]
# Generate system by supersizing ucell
system = ucell.supersize(*sizemults)
print('# of atoms in system =', system.natoms)
# of atoms in system = 4000
2.6. Specify calculation-specific run parameters
temperature temperature in Kelvin at which to run the MD integration scheme at. Default value is ‘0’.
runsteps specifies how many timesteps to integrate the system. Default value is 200000.
thermosteps specifies how often LAMMPS prints the system-wide thermo data. Default value is runsteps/1000, or 1 if runsteps is less than 1000.
dumpsteps specifies how often LAMMPS saves the atomic configuration to a LAMMPS dump file. Default value is runsteps, meaning only the first and last states are saved.
equilsteps specifies how many timesteps are ignored as equilibration time when computing the mean box parameters. Default value is 10000.
randomseed provides a random number seed to generating the initial atomic velocities. Default value gives a random number as the seed.
[8]:
temperature = 1900
runsteps = 200000
thermosteps = 100
dumpsteps = runsteps
equilsteps = 20000
randomseed = None
3. Define calculation function(s) and generate template LAMMPS script(s)
3.1. diffusion.template
[9]:
with open('diffusion.template', 'w') as f:
f.write("""# LAMMPS input script for dynamic msd computation
box tilt large
<atomman_system_info>
<atomman_pair_info>
# Assign simulation parameter values
variable temperature equal <temperature>
variable randomseed equal <randomseed>
variable thermosteps equal <thermosteps>
variable timestep equal <timestep>
variable equilsteps equal <equilsteps>
variable dumpsteps equal <dumpsteps>
variable runsteps equal <runsteps>
variable twotemp equal 2*${temperature}
variable damptemp equal 100*${timestep}
# Specify property computes
compute peatom all pe/atom
compute msd all msd com yes
# Define thermo data
thermo ${thermosteps}
thermo_style custom step temp pe pxx pyy pzz c_msd[1] c_msd[2] c_msd[3] c_msd[4]
thermo_modify format float %.13e
# Specify timestep
timestep ${timestep}
# Create velocities and equilibrate system using nvt
velocity all create ${twotemp} ${randomseed}
fix 1 all nvt temp ${temperature} ${temperature} ${damptemp}
run ${equilsteps}
unfix 1
<dump_info>
# Scale velocities to wanted temperature and run nve
velocity all scale ${temperature}
reset_timestep 0
fix 2 all nve
run ${runsteps}""")
3.2. pointdiffusion()
[10]:
def pointdiffusion(lammps_command, system, potential, point_kwargs,
mpi_command=None, temperature=300,
runsteps=200000, thermosteps=None, dumpsteps=0,
equilsteps=20000, randomseed=None):
"""
Evaluates the diffusion rate of a point defect at a given temperature. This
method will run two simulations: an NVT run at the specified temperature to
equilibrate the system, then an NVE run to measure the defect's diffusion
rate. The diffusion rate is evaluated using the mean squared displacement of
all atoms in the system, and using the assumption that diffusion is only due
to the added defect(s).
Parameters
----------
lammps_command :str
Command for running LAMMPS.
system : atomman.System
The system to perform the calculation on.
potential : atomman.lammps.Potential
The LAMMPS implemented potential to use.
point_kwargs : dict or list of dict
One or more dictionaries containing the keyword arguments for
the atomman.defect.point() function to generate specific point
defect configuration(s).
mpi_command : str, optional
The MPI command for running LAMMPS in parallel. If not given, LAMMPS
will run serially.
temperature : float, optional
The temperature to run at (default is 300.0).
runsteps : int, optional
The number of integration steps to perform (default is 200000).
thermosteps : int, optional
Thermo values will be reported every this many steps (default is
100).
dumpsteps : int or None, optional
Dump files will be saved every this many steps (default is 0,
which does not output dump files).
equilsteps : int, optional
The number of timesteps at the beginning of the simulation to
exclude when computing average values (default is 20000).
randomseed : int or None, optional
Random number seed used by LAMMPS in creating velocities and with
the Langevin thermostat. (Default is None which will select a
random int between 1 and 900000000.)
Returns
-------
dict
Dictionary of results consisting of keys:
- **'natoms'** (*int*) - The number of atoms in the system.
- **'temp'** (*float*) - The mean measured temperature.
- **'pxx'** (*float*) - The mean measured normal xx pressure.
- **'pyy'** (*float*) - The mean measured normal yy pressure.
- **'pzz'** (*float*) - The mean measured normal zz pressure.
- **'Epot'** (*numpy.array*) - The mean measured total potential
energy.
- **'temp_std'** (*float*) - The standard deviation in the measured
temperature values.
- **'pxx_std'** (*float*) - The standard deviation in the measured
normal xx pressure values.
- **'pyy_std'** (*float*) - The standard deviation in the measured
normal yy pressure values.
- **'pzz_std'** (*float*) - The standard deviation in the measured
normal zz pressure values.
- **'Epot_std'** (*float*) - The standard deviation in the measured
total potential energy values.
- **'dx'** (*float*) - The computed diffusion constant along the
x-direction.
- **'dy'** (*float*) - The computed diffusion constant along the
y-direction.
- **'dz'** (*float*) - The computed diffusion constant along the
y-direction.
- **'d'** (*float*) - The total computed diffusion constant.
"""
try:
# Get script's location if __file__ exists
script_dir = Path(__file__).parent
except:
# Use cwd otherwise
script_dir = Path.cwd()
# Add defect(s) to the initially perfect system
if not isinstance(point_kwargs, (list, tuple)):
point_kwargs = [point_kwargs]
for pkwargs in point_kwargs:
system = am.defect.point(system, **pkwargs)
# Get lammps units
lammps_units = lmp.style.unit(potential.units)
#Get lammps version date
lammps_date = lmp.checkversion(lammps_command)['date']
# Check that temperature is greater than zero
if temperature <= 0.0:
raise ValueError('Temperature must be greater than zero')
# Handle default values
if thermosteps is None:
thermosteps = runsteps // 1000
if thermosteps == 0:
thermosteps = 1
if dumpsteps is None:
dumpsteps = runsteps
if randomseed is None:
randomseed = random.randint(1, 900000000)
# Define lammps variables
lammps_variables = {}
system_info = system.dump('atom_data', f='initial.dat',
units=potential.units,
atom_style=potential.atom_style)
lammps_variables['atomman_system_info'] = system_info
lammps_variables['atomman_pair_info'] = potential.pair_info(system.symbols)
lammps_variables['temperature'] = temperature
lammps_variables['runsteps'] = runsteps
lammps_variables['equilsteps'] = equilsteps
lammps_variables['thermosteps'] = thermosteps
lammps_variables['dumpsteps'] = dumpsteps
lammps_variables['randomseed'] = randomseed
lammps_variables['timestep'] = lmp.style.timestep(potential.units)
# Set dump_info
if dumpsteps == 0:
lammps_variables['dump_info'] = ''
else:
lammps_variables['dump_info'] = '\n'.join([
'',
'# Define dump files',
'dump dumpit all custom ${dumpsteps} *.dump id type x y z c_peatom',
'dump_modify dumpit format <dump_modify_format>',
'',
])
# Set dump_modify_format based on lammps_date
if lammps_date < datetime.date(2016, 8, 3):
lammps_variables['dump_modify_format'] = '"%d %d %.13e %.13e %.13e %.13e"'
else:
lammps_variables['dump_modify_format'] = 'float %.13e'
# Write lammps input script
template_file = Path(script_dir, 'diffusion.template')
lammps_script = 'diffusion.in'
with open(template_file) as f:
template = f.read()
with open(lammps_script, 'w') as f:
f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))
# Run lammps
output = lmp.run(lammps_command, 'diffusion.in', mpi_command)
# Extract LAMMPS thermo data.
thermo = output.simulations[1]['thermo']
temps = thermo.Temp.values
pxxs = uc.set_in_units(thermo.Pxx.values, lammps_units['pressure'])
pyys = uc.set_in_units(thermo.Pyy.values, lammps_units['pressure'])
pzzs = uc.set_in_units(thermo.Pzz.values, lammps_units['pressure'])
potengs = uc.set_in_units(thermo.PotEng.values, lammps_units['energy'])
steps = thermo.Step.values
# Read user-defined thermo data
if output.lammps_date < datetime.date(2016, 8, 1):
msd_x = uc.set_in_units(thermo['msd[1]'].values,
lammps_units['length']+'^2')
msd_y = uc.set_in_units(thermo['msd[2]'].values,
lammps_units['length']+'^2')
msd_z = uc.set_in_units(thermo['msd[3]'].values,
lammps_units['length']+'^2')
msd = uc.set_in_units(thermo['msd[4]'].values,
lammps_units['length']+'^2')
else:
msd_x = uc.set_in_units(thermo['c_msd[1]'].values,
lammps_units['length']+'^2')
msd_y = uc.set_in_units(thermo['c_msd[2]'].values,
lammps_units['length']+'^2')
msd_z = uc.set_in_units(thermo['c_msd[3]'].values,
lammps_units['length']+'^2')
msd = uc.set_in_units(thermo['c_msd[4]'].values,
lammps_units['length']+'^2')
# Initialize results dict
results = {}
results['natoms'] = system.natoms
# Get mean and std for temperature, pressure, and potential energy
results['temp'] = np.mean(temps)
results['temp_std'] = np.std(temps)
results['pxx'] = np.mean(pxxs)
results['pxx_std'] = np.std(pxxs)
results['pyy'] = np.mean(pyys)
results['pyy_std'] = np.std(pyys)
results['pzz'] = np.mean(pzzs)
results['pzz_std'] = np.std(pzzs)
results['Epot'] = np.mean(potengs)
results['Epot_std'] = np.std(potengs)
# Convert steps to times
times = steps * uc.set_in_units(lammps_variables['timestep'], lammps_units['time'])
# Estimate diffusion rates
# MSD_ptd = natoms * MSD_atoms (if one defect in system)
# MSD = 2 * ndim * D * t --> D = MSD/t / (2 * ndim)
mx = np.polyfit(times, system.natoms * msd_x, 1)[0]
my = np.polyfit(times, system.natoms * msd_y, 1)[0]
mz = np.polyfit(times, system.natoms * msd_z, 1)[0]
m = np.polyfit(times, system.natoms * msd, 1)[0]
results['dx'] = mx / 2
results['dy'] = my / 2
results['dz'] = mz / 2
results['d'] = m / 6
return results
4. Run calculation function(s)
4.1. Generate point defect system and evaluate the energy
[11]:
results_dict = pointdiffusion(lammps_command, system, potential, point_kwargs,
mpi_command = mpi_command,
temperature = temperature,
runsteps = runsteps,
thermosteps = thermosteps,
dumpsteps = dumpsteps,
equilsteps = equilsteps,
randomseed = randomseed)
[12]:
results_dict.keys()
[12]:
dict_keys(['natoms', 'temp', 'temp_std', 'pxx', 'pxx_std', 'pyy', 'pyy_std', 'pzz', 'pzz_std', 'Epot', 'Epot_std', 'dx', 'dy', 'dz', 'd'])
5. Report results
5.1 Define units for outputting values
length2_pertime_unit is the units to display the computed diffusion constants in.
[13]:
length2_pertime_unit = 'm^2/s'
5.2 Print directional diffusion rates, \(D_x\), \(D_y\), and \(D_z\) and total diffusion rate, \(D\)
[14]:
print('Dx =', uc.get_in_units(results_dict['dx'], length2_pertime_unit), length2_pertime_unit)
print('Dy =', uc.get_in_units(results_dict['dy'], length2_pertime_unit), length2_pertime_unit)
print('Dz =', uc.get_in_units(results_dict['dz'], length2_pertime_unit), length2_pertime_unit)
print('D = ', uc.get_in_units(results_dict['d'], length2_pertime_unit), length2_pertime_unit)
Dx = 1.289386392173063e-10 m^2/s
Dy = 1.8253069489194243e-10 m^2/s
Dz = 1.2757396287269496e-10 m^2/s
D = 1.4634776566064903e-10 m^2/s