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.

Additional dependencies

Disclaimers

  • NIST 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.

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 shutil
import datetime
from copy import deepcopy

# http://www.numpy.org/
import numpy as np

# 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 2020-09-22 using iprPy version 0.10.2

1.2. Default calculation setup

[2]:
# Specify calculation style
calc_style = 'relax_static'

# 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)

# Initialize connection to library
library = iprPy.Library(load=['lammps_potentials'])

2. Assign values for the calculation’s run parameters

2.1. Specify 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.

[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 is an atomman.lammps.Potential object (required).

[4]:
potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'

# Retrieve potential and parameter file(s)
potential = library.get_lammps_potential(id=potential_name, getfiles=True)

2.3. Load 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.

[5]:
# 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

2.4. 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).

[6]:
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

2.5. Specify calculation-specific run 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.

[7]:
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

3. Define calculation function(s) and generate template LAMMPS script(s)

3.1. minbox.template

[8]:
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.2. relax_static()

[9]:
def relax_static(lammps_command, system, potential, mpi_command=None,
                 p_xx=0.0, p_yy=0.0, p_zz=0.0, p_xy=0.0, p_xz=0.0, p_yz=0.0,
                 dispmult=0.0, etol=0.0, ftol=0.0,  maxiter=10000,
                 maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom'),
                 maxcycles=100, ctol=1e-10):
    """
    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:

        - **'relaxed_system'** (*float*) - The relaxed system.
        - **'E_coh'** (*float*) - The cohesive energy of the relaxed system.
        - **'measured_pxx'** (*float*) - The measured x tensile pressure of the
          relaxed system.
        - **'measured_pyy'** (*float*) - The measured y tensile pressure of the
          relaxed system.
        - **'measured_pzz'** (*float*) - The measured z tensile pressure of the
          relaxed system.
        - **'measured_pxy'** (*float*) - The measured xy shear pressure of the
          relaxed system.
        - **'measured_pxz'** (*float*) - The measured xz shear pressure of the
          relaxed system.
        - **'measured_pyz'** (*float*) - The measured yz shear pressure of the
          relaxed system.
    """
    # Build filedict if function was called from iprPy
    try:
        assert __name__ == pkg_name
        calc = iprPy.load_calculation(calculation_style)
        filedict = calc.filedict
    except:
        filedict = {}

    # 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,
                                  return_pair_info=True)
        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
        template_file = 'minbox.template'
        lammps_script = 'minbox.in'
        template = iprPy.tools.read_calc_file(template_file, filedict)
        with open(lammps_script, 'w') as f:
            f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))

        # Run LAMMPS and extract thermo data
        logfile = 'log-' + str(cycle) + '.lammps'
        output = lmp.run(lammps_command, lammps_script, 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_coh'] = 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

4. Run calculation function(s)

[10]:
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)
[11]:
results_dict.keys()
[11]:
dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'E_coh', 'lx', 'ly', 'lz', 'xy', 'xz', 'yz', 'measured_pxx', 'measured_pyy', 'measured_pzz', 'measured_pxy', 'measured_pxz', 'measured_pyz'])

5. Report results

5.1. Define units for outputting values

  • length_unit is the unit of length to display values in.

  • energy_unit is the unit of energy to display values in.

  • pressure_unit is the unit of pressure to display values in.

[12]:
length_unit = 'angstrom'
energy_unit = 'eV'
pressure_unit = 'GPa'

5.2. Print Ecoh and lattice constants of relaxed ucell

[13]:
print('Ecoh =', uc.get_in_units(results_dict['E_coh'], energy_unit), energy_unit)

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'])

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)
Ecoh = -4.449999998325555 eV
a = 3.51999946361013 angstrom
b = 3.5199994969059634 angstrom
c = 3.519999505001706 angstrom
alpha = 90.0
beta =  90.0
gamma = 90.0

5.3. Check final system pressures

[14]:
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 = -7.1585342727316e-06 GPa
Pyy = -8.1046080506064e-06 GPa
Pzz = -8.3345122489377e-06 GPa
Pyz = 3.6844223558997003e-10 GPa
Pxz = 1.6824819411886998e-10 GPa
Pxy = 3.9876281320881e-11 GPa

5.4. Show relaxed atomic configuration

[15]:
finalsystem = am.load('atom_dump', results_dict['dumpfile_final'],
                      symbols=results_dict['symbols_final'])
print(finalsystem)
avect =  [10.560,  0.000,  0.000]
bvect =  [ 0.000, 10.560,  0.000]
cvect =  [ 0.000,  0.000, 10.560]
origin = [-0.030, -0.030, -0.030]
natoms = 108
natypes = 1
symbols = ('Ni',)
pbc = [ True  True  True]
per-atom properties = ['atype', 'pos', 'atom_id', 'c_peatom']
     id   atype  pos[0]  pos[1]  pos[2]
      0       1  10.530  -0.030  10.530
      1       1  10.530   1.730   1.730
      2       1   1.730  -0.030   1.730
      3       1   1.730   1.730  10.530
      4       1   3.490  10.530  10.530
      5       1   3.490   1.730   1.730
      6       1   5.250  -0.030   1.730
      7       1   5.250   1.730  -0.030
      8       1   7.010  -0.030  -0.030
      9       1   7.010   1.730   1.730
     10       1   8.770  10.530   1.730
     11       1   8.770   1.730  -0.030
     12       1  -0.030   3.490  -0.030
     13       1  10.530   5.250   1.730
     14       1   1.730   3.490   1.730
     15       1   1.730   5.250  10.530
     16       1   3.490   3.490  -0.030
     17       1   3.490   5.250   1.730
     18       1   5.250   3.490   1.730
     19       1   5.250   5.250  10.530
     20       1   7.010   3.490  -0.030
     21       1   7.010   5.250   1.730
     22       1   8.770   3.490   1.730
     23       1   8.770   5.250  -0.030
     24       1  -0.030   7.010  10.530
     25       1  10.530   8.770   1.730
     26       1   1.730   7.010   1.730
     27       1   1.730   8.770  -0.030
     28       1   3.490   7.010  -0.030
     29       1   3.490   8.770   1.730
     30       1   5.250   7.010   1.730
     31       1   5.250   8.770  -0.030
     32       1   7.010   7.010  10.530
     33       1   7.010   8.770   1.730
     34       1   8.770   7.010   1.730
     35       1   8.770   8.770  -0.030
     36       1  10.530  10.530   3.490
     37       1  -0.030   1.730   5.250
     38       1   1.730  -0.030   5.250
     39       1   1.730   1.730   3.490
     40       1   3.490  -0.030   3.490
     41       1   3.490   1.730   5.250
     42       1   5.250  -0.030   5.250
     43       1   5.250   1.730   3.490
     44       1   7.010  -0.030   3.490
     45       1   7.010   1.730   5.250
     46       1   8.770  -0.030   5.250
     47       1   8.770   1.730   3.490
     48       1  10.530   3.490   3.490
     49       1  -0.030   5.250   5.250
     50       1   1.730   3.490   5.250
     51       1   1.730   5.250   3.490
     52       1   3.490   3.490   3.490
     53       1   3.490   5.250   5.250
     54       1   5.250   3.490   5.250
     55       1   5.250   5.250   3.490
     56       1   7.010   3.490   3.490
     57       1   7.010   5.250   5.250
     58       1   8.770   3.490   5.250
     59       1   8.770   5.250   3.490
     60       1  -0.030   7.010   3.490
     61       1  -0.030   8.770   5.250
     62       1   1.730   7.010   5.250
     63       1   1.730   8.770   3.490
     64       1   3.490   7.010   3.490
     65       1   3.490   8.770   5.250
     66       1   5.250   7.010   5.250
     67       1   5.250   8.770   3.490
     68       1   7.010   7.010   3.490
     69       1   7.010   8.770   5.250
     70       1   8.770   7.010   5.250
     71       1   8.770   8.770   3.490
     72       1  -0.030  10.530   7.010
     73       1  10.530   1.730   8.770
     74       1   1.730  -0.030   8.770
     75       1   1.730   1.730   7.010
     76       1   3.490  10.530   7.010
     77       1   3.490   1.730   8.770
     78       1   5.250  -0.030   8.770
     79       1   5.250   1.730   7.010
     80       1   7.010  -0.030   7.010
     81       1   7.010   1.730   8.770
     82       1   8.770  -0.030   8.770
     83       1   8.770   1.730   7.010
     84       1  -0.030   3.490   7.010
     85       1  -0.030   5.250   8.770
     86       1   1.730   3.490   8.770
     87       1   1.730   5.250   7.010
     88       1   3.490   3.490   7.010
     89       1   3.490   5.250   8.770
     90       1   5.250   3.490   8.770
     91       1   5.250   5.250   7.010
     92       1   7.010   3.490   7.010
     93       1   7.010   5.250   8.770
     94       1   8.770   3.490   8.770
     95       1   8.770   5.250   7.010
     96       1  -0.030   7.010   7.010
     97       1  -0.030   8.770   8.770
     98       1   1.730   7.010   8.770
     99       1   1.730   8.770   7.010
    100       1   3.490   7.010   7.010
    101       1   3.490   8.770   8.770
    102       1   5.250   7.010   8.770
    103       1   5.250   8.770   7.010
    104       1   7.010   7.010   7.010
    105       1   7.010   8.770   8.770
    106       1   8.770   7.010   8.770
    107       1   8.770   8.770   7.010