relax_static - Methodology and code
Python imports
[1]:
# Standard library imports
from pathlib import Path
import shutil
import datetime
from copy import deepcopy
from math import floor
from typing import Optional, Tuple
# http://www.numpy.org/
import numpy as np
# https://ipython.org/
from IPython.display import display, Markdown
# https://github.com/usnistgov/atomman
import atomman as am
import atomman.lammps as lmp
import atomman.unitconvert as uc
from atomman.tools import filltemplate
# https://github.com/usnistgov/iprPy
import iprPy
from iprPy.tools import read_calc_file
print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)
Notebook last executed on 2022-03-07 using iprPy version 0.11.2
1. Load calculation and view description
1.1. Load the calculation
[2]:
# Load the calculation being demoed
calculation = iprPy.load_calculation('relax_static')
1.2. Display calculation description and theory
[3]:
# Display main docs and theory
display(Markdown(calculation.maindoc))
display(Markdown(calculation.theorydoc))
relax_static calculation style
Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.
Introduction
The relax_static calculation style uses static energy/force minimizations to relax the atomic positions and box dimensions of a system to a specified pressure.
Version notes
2018-07-09: Notebook added.
2019-07-30: Description updated and small changes due to iprPy version.
2020-05-22: Version 0.10 update - potentials now loaded from database.
2020-09-22: Setup and parameter definition streamlined.
2022-03-04: Notebook updated to reflect version 0.11.
Additional dependencies
Disclaimers
The minimization algorithm will drive the system to a local minimum, which may not be the global minimum. There is no guarantee that the resulting structure is dynamically stable, and it is possible that the relaxation of certain dimensions may be constrained to move together during the minimization preventing a full relaxation.
Method and Theory
This method uses the LAMMPS minimization plus box_relax commands to simultaneously relax both the atomic positions and the system’s box dimensions towards a local minimum. The LAMMPS documentation of the box_relax command notes that the complete minimization algorithm is not well defined which may prevent a complete relaxation during a single run. To overcome this limitation, the calculation script continuously restarts the minimization until the box dimensions from one run to the next remain within a specified tolerance.
2. Define calculation functions and generate files
This section defines the calculation functions and associated resource files exactly as they exist inside the iprPy package. This allows for the code used to be directly visible and modifiable by anyone looking to see how it works.
2.1. relax_static()
This is the primary function for the calculation. The version of this function built in iprPy can be accessed by calling the calc() method of an object of the associated calculation class.
[4]:
def relax_static(lammps_command: str,
system: am.System,
potential: lmp.Potential,
mpi_command: Optional[str] = None,
p_xx: float = 0.0,
p_yy: float = 0.0,
p_zz: float = 0.0,
p_xy: float = 0.0,
p_xz: float = 0.0,
p_yz: float = 0.0,
dispmult: float = 0.0,
etol: float = 0.0,
ftol: float = 0.0,
maxiter: int = 100000,
maxeval: int = 1000000,
dmax: float = uc.set_in_units(0.01, 'angstrom'),
maxcycles: int = 100,
ctol: float = 1e-10) -> dict:
"""
Repeatedly runs the ELASTIC example distributed with LAMMPS until box
dimensions converge within a tolerance.
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.
mpi_command : str, optional
The MPI command for running LAMMPS in parallel. If not given, LAMMPS
will run serially.
p_xx : float, optional
The value to relax the x tensile pressure component to (default is
0.0).
p_yy : float, optional
The value to relax the y tensile pressure component to (default is
0.0).
p_zz : float, optional
The value to relax the z tensile pressure component to (default is
0.0).
p_xy : float, optional
The value to relax the xy shear pressure component to (default is
0.0).
p_xz : float, optional
The value to relax the xz shear pressure component to (default is
0.0).
p_yz : float, optional
The value to relax the yz shear pressure component to (default is
0.0).
dispmult : float, optional
Multiplier for applying a random displacement to all atomic positions
prior to relaxing. Default value is 0.0.
etol : float, optional
The energy tolerance for the structure minimization. This value is
unitless. (Default is 0.0).
ftol : float, optional
The force tolerance for the structure minimization. This value is in
units of force. (Default is 0.0).
maxiter : int, optional
The maximum number of minimization iterations to use (default is 10000).
maxeval : int, optional
The maximum number of minimization evaluations to use (default is
100000).
dmax : float, optional
The maximum distance in length units that any atom is allowed to relax
in any direction during a single minimization iteration (default is
0.01 Angstroms).
pressure_unit : str, optional
The unit of pressure to calculate the elastic constants in (default is
'GPa').
maxcycles : int, optional
The maximum number of times the minimization algorithm is called.
Default value is 100.
ctol : float, optional
The relative tolerance used to determine if the lattice constants have
converged (default is 1e-10).
Returns
-------
dict
Dictionary of results consisting of keys:
- **'dumpfile_initial'** (*str*) - The name of the initial dump file
created.
- **'symbols_initial'** (*list*) - The symbols associated with the
initial dump file.
- **'dumpfile_final'** (*str*) - The name of the final dump file
created.
- **'symbols_final'** (*list*) - The symbols associated with the final
dump file.
- **'lx'** (*float*) - The relaxed lx box length.
- **'ly'** (*float*) - The relaxed ly box length.
- **'lz'** (*float*) - The relaxed lz box length.
- **'xy'** (*float*) - The relaxed xy box tilt.
- **'xz'** (*float*) - The relaxed xz box tilt.
- **'yz'** (*float*) - The relaxed yz box tilt.
- **'E_pot'** (*float*) - The potential energy per atom for the final
configuration.
- **'measured_pxx'** (*float*) - The measured x tensile pressure
component for the final configuration.
- **'measured_pyy'** (*float*) - The measured y tensile pressure
component for the final configuration.
- **'measured_pzz'** (*float*) - The measured z tensile pressure
component for the final configuration.
- **'measured_pxy'** (*float*) - The measured xy shear pressure
component for the final configuration.
- **'measured_pxz'** (*float*) - The measured xz shear pressure
component for the final configuration.
- **'measured_pyz'** (*float*) - The measured yz shear pressure
component for the final configuration.
"""
# Get lammps units
lammps_units = lmp.style.unit(potential.units)
# Get lammps version date
lammps_date = lmp.checkversion(lammps_command)['date']
# Save initial configuration as a dump file
system.dump('atom_dump', f='initial.dump')
# Apply small random distortions to atoms
system.atoms.pos += dispmult * np.random.rand(*system.atoms.pos.shape) - dispmult / 2
# Initialize parameters
old_vects = system.box.vects
converged = False
# Run minimizations up to maxcycles times
for cycle in range(maxcycles):
# Define lammps variables
lammps_variables = {}
system_info = system.dump('atom_data', f='init.dat',
potential=potential)
lammps_variables['atomman_system_pair_info'] = system_info
lammps_variables['p_xx'] = uc.get_in_units(p_xx, lammps_units['pressure'])
lammps_variables['p_yy'] = uc.get_in_units(p_yy, lammps_units['pressure'])
lammps_variables['p_zz'] = uc.get_in_units(p_zz, lammps_units['pressure'])
lammps_variables['p_xy'] = uc.get_in_units(p_xy, lammps_units['pressure'])
lammps_variables['p_xz'] = uc.get_in_units(p_xz, lammps_units['pressure'])
lammps_variables['p_yz'] = uc.get_in_units(p_yz, lammps_units['pressure'])
lammps_variables['etol'] = etol
lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])
lammps_variables['maxiter'] = maxiter
lammps_variables['maxeval'] = maxeval
lammps_variables['dmax'] = uc.get_in_units(dmax, lammps_units['length'])
# Set dump_keys based on atom_style
if potential.atom_style in ['charge']:
lammps_variables['dump_keys'] = 'id type q x y z c_peatom'
else:
lammps_variables['dump_keys'] = 'id type x y z c_peatom'
# Set dump_modify_format based on lammps_date
if lammps_date < datetime.date(2016, 8, 3):
if potential.atom_style in ['charge']:
lammps_variables['dump_modify_format'] = '"%d %d %.13e %.13e %.13e %.13e %.13e"'
else:
lammps_variables['dump_modify_format'] = '"%d %d %.13e %.13e %.13e %.13e"'
else:
lammps_variables['dump_modify_format'] = 'float %.13e'
# Write lammps input script
lammps_script = 'minbox.in'
template = read_calc_file('iprPy.calculation.relax_static',
'minbox.template')
with open(lammps_script, 'w') as f:
f.write(filltemplate(template, lammps_variables, '<', '>'))
# Run LAMMPS and extract thermo data
logfile = 'log-' + str(cycle) + '.lammps'
output = lmp.run(lammps_command, script_name=lammps_script,
mpi_command=mpi_command, logfile=logfile)
thermo = output.simulations[0]['thermo']
# Clean up dump files
Path('0.dump').unlink()
last_dump_file = str(thermo.Step.values[-1]) + '.dump'
renamed_dump_file = 'relax_static-' + str(cycle) + '.dump'
shutil.move(last_dump_file, renamed_dump_file)
# Load relaxed system
system = am.load('atom_dump', renamed_dump_file, symbols=system.symbols)
# Test if box dimensions have converged
if np.allclose(old_vects, system.box.vects, rtol=ctol, atol=0):
converged = True
break
else:
old_vects = system.box.vects
# Check for convergence
if converged is False:
raise RuntimeError('Failed to converge after ' + str(maxcycles) + ' cycles')
# Zero out near-zero tilt factors
lx = system.box.lx
ly = system.box.ly
lz = system.box.lz
xy = system.box.xy
xz = system.box.xz
yz = system.box.yz
if np.isclose(xy/ly, 0.0, rtol=0.0, atol=1e-10):
xy = 0.0
if np.isclose(xz/lz, 0.0, rtol=0.0, atol=1e-10):
xz = 0.0
if np.isclose(yz/lz, 0.0, rtol=0.0, atol=1e-10):
yz = 0.0
system.box.set(lx=lx, ly=ly, lz=lz, xy=xy, xz=xz, yz=yz)
system.wrap()
# Build results_dict
results_dict = {}
results_dict['dumpfile_initial'] = 'initial.dump'
results_dict['symbols_initial'] = system.symbols
results_dict['dumpfile_final'] = renamed_dump_file
results_dict['symbols_final'] = system.symbols
results_dict['E_pot'] = uc.set_in_units(thermo.PotEng.values[-1] / system.natoms,
lammps_units['energy'])
results_dict['lx'] = uc.set_in_units(lx, lammps_units['length'])
results_dict['ly'] = uc.set_in_units(ly, lammps_units['length'])
results_dict['lz'] = uc.set_in_units(lz, lammps_units['length'])
results_dict['xy'] = uc.set_in_units(xy, lammps_units['length'])
results_dict['xz'] = uc.set_in_units(xz, lammps_units['length'])
results_dict['yz'] = uc.set_in_units(yz, lammps_units['length'])
results_dict['measured_pxx'] = uc.set_in_units(thermo.Pxx.values[-1],
lammps_units['pressure'])
results_dict['measured_pyy'] = uc.set_in_units(thermo.Pyy.values[-1],
lammps_units['pressure'])
results_dict['measured_pzz'] = uc.set_in_units(thermo.Pzz.values[-1],
lammps_units['pressure'])
results_dict['measured_pxy'] = uc.set_in_units(thermo.Pxy.values[-1],
lammps_units['pressure'])
results_dict['measured_pxz'] = uc.set_in_units(thermo.Pxz.values[-1],
lammps_units['pressure'])
results_dict['measured_pyz'] = uc.set_in_units(thermo.Pyz.values[-1],
lammps_units['pressure'])
return results_dict
2.2. minbox.template file
[5]:
with open('minbox.template', 'w') as f:
f.write("""# LAMMPS input script that performs an energy minimization and box relaxation
box tilt large
<atomman_system_pair_info>
change_box all triclinic
thermo_style custom step lx ly lz xy xz yz pxx pyy pzz pxy pxz pyz pe
thermo_modify format float %.13e
compute peatom all pe/atom
dump dumpit all custom <maxeval> *.dump <dump_keys>
dump_modify dumpit format <dump_modify_format>
fix boxrelax all box/relax x <p_xx> y <p_yy> z <p_zz> xy <p_xy> xz <p_xz> yz <p_yz>
min_modify dmax <dmax>
minimize <etol> <ftol> <maxiter> <maxeval>""")
3. Specify input parameters
3.1. System-specific paths
lammps_command is the LAMMPS command to use (required).
mpi_command MPI command for running LAMMPS in parallel. A value of None will run simulations serially.
[6]:
lammps_command = 'lmp_serial'
mpi_command = None
3.2. Interatomic potential
potential_name gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.
potential is an atomman.lammps.Potential object (required).
[7]:
potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'
# Retrieve potential and parameter file(s) using atomman
potential = am.load_lammps_potential(id=potential_name, getfiles=True)
3.3. Initial unit cell system
ucell is an atomman.System representing a fundamental unit cell of the system (required). Here, this is generated using the load parameters and symbols.
[8]:
# Create ucell by loading prototype record
ucell = am.load('prototype', 'A1--Cu--fcc', symbols='Ni', a=3.5)
print(ucell)
avect = [ 3.500, 0.000, 0.000]
bvect = [ 0.000, 3.500, 0.000]
cvect = [ 0.000, 0.000, 3.500]
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.750 1.750
2 1 1.750 0.000 1.750
3 1 1.750 1.750 0.000
3.4. System modifications
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).
[9]:
sizemults = [3, 3, 3]
# Generate system by supersizing ucell
system = ucell.supersize(*sizemults)
print('# of atoms in system =', system.natoms)
# of atoms in system = 108
3.5. Calculation-specific parameters
pressure_xx gives the xx component of the pressure to equilibriate the system to.
pressure_yy gives the yy component of the pressure to equilibriate the system to.
pressure_zz gives the zz component of the pressure to equilibriate the system to.
pressure_xy gives the xy component of the pressure to equilibriate the system to.
pressure_xz gives the xz component of the pressure to equilibriate the system to.
pressure_yz gives the yz component of the pressure to equilibriate the system to.
displacementkick specifies a multiplier for a random shift of atomic positions to apply prior to relaxation. This is in length units.
energytolerance is the energy tolerance to use during the minimizations. This is unitless.
forcetolerance is the force tolerance to use during the minimizations. This is in energy/length units.
maxiterations is the maximum number of minimization iterations to use.
maxevaluations is the maximum number of minimization evaluations to use.
maxatommotion is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units.
maxcycles is the maximum number of minimization runs (cycles) to perform.
cycletolerance is the relative tolerance to use in identifying if the lattice constants have converged from one cycle to the next.
[10]:
pressure_xx = uc.set_in_units(0.0, 'GPa')
pressure_yy = uc.set_in_units(0.0, 'GPa')
pressure_zz = uc.set_in_units(0.0, 'GPa')
pressure_xy = uc.set_in_units(0.0, 'GPa')
pressure_xz = uc.set_in_units(0.0, 'GPa')
pressure_yz = uc.set_in_units(0.0, 'GPa')
displacementkick = uc.set_in_units(0.00001, 'angstrom')
energytolerance = 1e-8
forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')
maxiterations = 10000
maxevaluations = 100000
maxatommotion = uc.set_in_units(0.01, 'angstrom')
maxcycles = 100
cycletolerance = 1e-7
4. Run calculation and view results
4.1. Run calculation
All primary calculation method functions take a series of inputs and return a dictionary of outputs.
[11]:
results_dict = relax_static(lammps_command, system, potential,
mpi_command = mpi_command,
p_xx = pressure_xx,
p_yy = pressure_yy,
p_zz = pressure_zz,
p_xy = pressure_xy,
p_xz = pressure_xz,
p_yz = pressure_yz,
dispmult = displacementkick,
etol = energytolerance,
ftol = forcetolerance,
maxiter = maxiterations,
maxeval = maxevaluations,
dmax = maxatommotion,
maxcycles = maxcycles,
ctol = cycletolerance)
print(results_dict.keys())
dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'E_pot', 'lx', 'ly', 'lz', 'xy', 'xz', 'yz', 'measured_pxx', 'measured_pyy', 'measured_pzz', 'measured_pxy', 'measured_pxz', 'measured_pyz'])
4.2. Report results
Values returned in the results_dict:
‘dumpfile_initial’ (str) - The name of the initial dump file created.
‘symbols_initial’ (list) - The symbols associated with the initial dump file.
‘dumpfile_final’ (str) - The name of the final dump file created.
‘symbols_final’ (list) - The symbols associated with the final dump file.
‘lx’ (float) - The relaxed lx box length.
‘ly’ (float) - The relaxed ly box length.
‘lz’ (float) - The relaxed lz box length.
‘xy’ (float) - The relaxed xy box tilt.
‘xz’ (float) - The relaxed xz box tilt.
‘yz’ (float) - The relaxed yz box tilt.
‘E_pot’ (float) - The potential energy per atom for the final configuration.
‘measured_pxx’ (float) - The measured x tensile pressure component for the final configuration.
‘measured_pyy’ (float) - The measured y tensile pressure component for the final configuration.
‘measured_pzz’ (float) - The measured z tensile pressure component for the final configuration.
‘measured_pxy’ (float) - The measured xy shear pressure component for the final configuration.
‘measured_pxz’ (float) - The measured xz shear pressure component for the final configuration.
‘measured_pyz’ (float) - The measured yz shear pressure component for the final configuration.
[12]:
# Show initial and final dump files
print(results_dict['dumpfile_initial'])
print(results_dict['symbols_initial'])
print(results_dict['dumpfile_final'])
print(results_dict['symbols_final'])
initial.dump
('Ni',)
relax_static-16.dump
('Ni',)
[13]:
length_unit = 'angstrom'
energy_unit = 'eV'
# Show the per atom potential energy
print('E_pot =', uc.get_in_units(results_dict['E_pot'], energy_unit), energy_unit)
# Construct a Box from the returned system dimensions
box = am.Box(lx=results_dict['lx'], ly=results_dict['ly'], lz=results_dict['lz'],
xy=results_dict['xy'], xz=results_dict['xz'], yz=results_dict['yz'])
# Retrieve lattice constants by dividing by sizemults
print('a =', uc.get_in_units(box.a / sizemults[0], length_unit), length_unit)
print('b =', uc.get_in_units(box.b / sizemults[1], length_unit), length_unit)
print('c =', uc.get_in_units(box.c / sizemults[2], length_unit), length_unit)
print('alpha =', box.alpha)
print('beta = ', box.beta)
print('gamma =', box.gamma)
E_pot = -4.449999998328333 eV
a = 3.5199995783500473 angstrom
b = 3.5199994658760616 angstrom
c = 3.5199994012564813 angstrom
alpha = 90.0
beta = 90.0
gamma = 90.0
[14]:
pressure_unit = 'GPa'
# Show the computed pressure tensor
print('Pxx =', uc.get_in_units(results_dict['measured_pxx'], pressure_unit), pressure_unit)
print('Pyy =', uc.get_in_units(results_dict['measured_pyy'], pressure_unit), pressure_unit)
print('Pzz =', uc.get_in_units(results_dict['measured_pzz'], pressure_unit), pressure_unit)
print('Pyz =', uc.get_in_units(results_dict['measured_pyz'], pressure_unit), pressure_unit)
print('Pxz =', uc.get_in_units(results_dict['measured_pxz'], pressure_unit), pressure_unit)
print('Pxy =', uc.get_in_units(results_dict['measured_pxy'], pressure_unit), pressure_unit)
Pxx = -9.578224572466802e-06 GPa
Pyy = -6.381935436363099e-06 GPa
Pzz = -4.5453178638912e-06 GPa
Pyz = -1.8889229415301003e-10 GPa
Pxz = -1.2387278167398997e-10 GPa
Pxy = 1.0204743784825e-10 GPa