relax_dynamic - 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
import random

# 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 2023-07-31 using iprPy version 0.11.6

1. Load calculation and view description

1.1. Load the calculation

[2]:
# Load the calculation being demoed
calculation = iprPy.load_calculation('relax_dynamic')

1.2. Display calculation description and theory

[3]:
# Display main docs and theory
display(Markdown(calculation.maindoc))
display(Markdown(calculation.theorydoc))

relax_dynamic calculation style

Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.

Introduction

The relax_dynamic calculation style dynamically relaxes an atomic configuration for a specified number of timesteps. Upon completion, the mean, \(\langle X \rangle\), and standard deviation, \(\sigma_X\), of all thermo properties, \(X\), are computed for a specified range of times. This method is meant to measure equilibrium properties of bulk materials, both at zero K and at various temperatures.

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-11: Notebook updated to reflect version 0.11. Restart capability added in.

Additional dependencies
Disclaimers
  • NIST disclaimers

  • The calculation reports the standard deviation, \(\sigma_X\) of the measured properties not the standard error of the mean, \(\sigma_{\langle X \rangle}\). The two are related to each other according to \(\sigma_{\langle X \rangle} = \sigma_X \sqrt{\frac{C}{N}}\), where \(N\) is the number of samples taken of \(X\), and \(C\) is a statistical inefficiency due to the autocorrelation of the measurements with time. Obtaining a proper estimate of \(\sigma_{\langle X \rangle}\) requires either estimating \(C\) from the raw thermo data (not done here), or only taking measurements sporadically to ensure the samples are independent.

  • Good (low error) results requires running large simulations for a long time. The reasons for this are:

    • Systems have to be large enough to avoid issues with fluctuations across the periodic boundaries.

    • Runs must first let the systems equilibrate before meaningful measurements can be taken.

    • The standard deviation, \(\sigma\), of thermo properties is proportional to the number of atoms, \(N_a\) as \(\sigma \propto \frac{1}{\sqrt{N_a}}\).

    • The standard error, \(\sigma_x\) of thermo properties is proportional to the number of samples taken, \(N\) as \(\sigma_x \propto \frac{1}{\sqrt{N}}\).

Method and Theory

An initial system (and corresponding unit cell system) is supplied with box dimensions, \(a_i^0\), close to the equilibrium values. A LAMMPS simulation then integrates the atomic positions and velocities for a specified number of timesteps.

The calculation script allows for the use of different integration methods:

  • nve integrates atomic positions without changing box dimensions or the system’s total energy.

  • npt integrates atomic positions and applies Nose-Hoover style thermostat and barostat (equilibriate to specified T and P).

  • nvt integrates atomic positions and applies Nose-Hoover style thermostat (equilibriate to specified T).

  • nph integrates atomic positions and applies Nose-Hoover style barostat (equilibriate to specified P).

  • nve+l integrates atomic positions and applies Langevin style thermostat (equilibriate to specified T).

  • nph+l integrates atomic positions and applies Nose-Hoover style barostat and Langevin style thermostat (equilibriate to specified T and P).

Notes on the different control schemes:

  • The Nose-Hoover barostat works by rescaling the box dimensions according to the measured system pressures.

  • The Nose-Hoover thermostat works by rescaling the atomic velocities according to the measured system temperature (kinetic energy). Cannot be used with a temperature of 0 K.

  • The Langevin thermostat works by modifying the forces on all atoms with both a dampener and a random temperature dependent fluctuation. Used at 0 K, only the force dampener is applied.

Notes on run parameter values. The proper time to reach equilibrium (equilsteps), and sample frequency to ensure uncorrelated measurements (thermosteps) is simulation dependent. They can be influenced by the potential, timestep size, crystal structure, integration method, presence of defects, etc. The default values of equilsteps = 20,000 and thermosteps = 100 are based on general rule-of-thumb estimates for bulk crystals and EAM potentials, and may or may not be adequate.

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_dynamic()

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_dynamic(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,
                  temperature: float = 0.0,
                  integrator: Optional[str] = None,
                  runsteps: int = 220000,
                  thermosteps: int = 100,
                  dumpsteps: Optional[int] = None,
                  restartsteps: Optional[int] = None,
                  equilsteps: int = 20000,
                  randomseed: Optional[int] = None) -> dict:
    """
    Performs a full dynamic relax on a given system at the given temperature
    to the specified pressure state.

    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.
    symbols : list of str
        The list of element-model symbols for the Potential that correspond to
        system's atypes.
    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).
    temperature : float, optional
        The temperature to relax at (default is 0.0).
    runsteps : int, optional
        The number of integration steps to perform (default is 220000).
    integrator : str or None, optional
        The integration method to use. Options are 'npt', 'nvt', 'nph',
        'nve', 'nve+l', 'nph+l'. The +l options use Langevin thermostat.
        (Default is None, which will use 'nph+l' for temperature == 0, and
        'npt' otherwise.)
    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 None,
        which sets dumpsteps equal to runsteps).
    restartsteps : int or None, optional
        Restart files will be saved every this many steps (default is None,
        which sets restartsteps equal to runsteps).
    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:

        - **'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.
        - **'nsamples'** (*int*) - The number of thermodynamic samples included
          in the mean and standard deviation estimates.  Can also be used to
          estimate standard error values assuming that the thermo step size is
          large enough (typically >= 100) to assume the samples to be
          independent.
        - **'E_pot'** (*float*) - The mean measured potential energy.
        - **'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.
        - **'temp'** (*float*) - The mean measured temperature.
        - **'E_pot_std'** (*float*) - The standard deviation in the measured
          potential energy values.
        - **'measured_pxx_std'** (*float*) - The standard deviation in the
          measured x tensile pressure of the relaxed system.
        - **'measured_pyy_std'** (*float*) - The standard deviation in the
          measured y tensile pressure of the relaxed system.
        - **'measured_pzz_std'** (*float*) - The standard deviation in the
          measured z tensile pressure of the relaxed system.
        - **'measured_pxy_std'** (*float*) - The standard deviation in the
          measured xy shear pressure of the relaxed system.
        - **'measured_pxz_std'** (*float*) - The standard deviation in the
          measured xz shear pressure of the relaxed system.
        - **'measured_pyz_std'** (*float*) - The standard deviation in the
          measured yz shear pressure of the relaxed system.
        - **'temp_std'** (*float*) - The standard deviation in the measured
          temperature values.
    """

    # Get lammps units
    lammps_units = lmp.style.unit(potential.units)

    #Get lammps version date
    lammps_date = lmp.checkversion(lammps_command)['date']

    # Handle default values
    if dumpsteps is None:
        dumpsteps = runsteps
    if restartsteps is None:
        restartsteps = runsteps

    # Define lammps variables
    lammps_variables = {}

    # Dump initial system as data and build LAMMPS inputs
    system_info = system.dump('atom_data', f='init.dat',
                              potential=potential)
    lammps_variables['atomman_system_pair_info'] = system_info

    # Generate LAMMPS inputs for restarting
    system_info2 = potential.pair_restart_info('*.restart', system.symbols)
    lammps_variables['atomman_pair_restart_info'] = system_info2

    # Integrator lines for main run
    integ_info = integrator_info(integrator=integrator,
                                 p_xx=p_xx, p_yy=p_yy, p_zz=p_zz,
                                 p_xy=p_xy, p_xz=p_xz, p_yz=p_yz,
                                 temperature=temperature,
                                 randomseed=randomseed,
                                 units=potential.units)
    lammps_variables['integrator_info'] = integ_info

    # Integrator lines for restarts
    integ_info2 = integrator_info(integrator=integrator,
                                 p_xx=p_xx, p_yy=p_yy, p_zz=p_zz,
                                 p_xy=p_xy, p_xz=p_xz, p_yz=p_yz,
                                 temperature=temperature,
                                 velocity_temperature=0.0,
                                 randomseed=randomseed,
                                 units=potential.units)
    lammps_variables['integrator_restart_info'] = integ_info2

    # Other run settings
    lammps_variables['thermosteps'] = thermosteps
    lammps_variables['runsteps'] = runsteps
    lammps_variables['dumpsteps'] = dumpsteps
    lammps_variables['restartsteps'] = restartsteps

    # Set compute stress/atom based on LAMMPS version
    if lammps_date < datetime.date(2014, 2, 12):
        lammps_variables['stressterm'] = ''
    else:
        lammps_variables['stressterm'] = 'NULL'

    # Set dump_keys based on atom_style
    if potential.atom_style in ['charge']:
        lammps_variables['dump_keys'] = 'id type q xu yu zu c_pe c_ke &\n'
        lammps_variables['dump_keys'] += 'c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]'
    else:
        lammps_variables['dump_keys'] = 'id type xu yu zu c_pe c_ke &\n'
        lammps_variables['dump_keys'] += 'c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]'

    # 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'] = f'"%d %d{12 * " %.13e"}"'
        else:
            lammps_variables['dump_modify_format'] = f'"%d %d{11 * " %.13e"}"'
    else:
        lammps_variables['dump_modify_format'] = 'float %.13e'

    # Write lammps input script
    lammps_script = 'full_relax.in'
    template = read_calc_file('iprPy.calculation.relax_dynamic',
                              'full_relax.template')
    with open(lammps_script, 'w') as f:
        f.write(filltemplate(template, lammps_variables, '<', '>'))

    # Write lammps restart input script
    restart_script = 'full_relax_restart.in'
    template = read_calc_file('iprPy.calculation.relax_dynamic',
                              'full_relax_restart.template')
    with open(restart_script, 'w') as f:
        f.write(filltemplate(template, lammps_variables, '<', '>'))

    # Run lammps
    output = lmp.run(lammps_command, script_name=lammps_script,
                     restart_script_name=restart_script,
                     mpi_command=mpi_command, screen=False)

    # Extract LAMMPS thermo data.
    results = {}
    thermo = output.flatten()['thermo']

    results['dumpfile_initial'] = '0.dump'
    results['symbols_initial'] = system.symbols

    # Load relaxed system from dump file
    last_dump_file = f'{thermo.Step.values[-1]}.dump'
    results['dumpfile_final'] = last_dump_file
    system = am.load('atom_dump', last_dump_file, symbols=system.symbols)
    results['symbols_final'] = system.symbols

    # Only consider values where Step >= equilsteps
    thermo = thermo[thermo.Step >= equilsteps]
    results['nsamples'] = len(thermo)

    # Get cohesive energy estimates
    natoms = system.natoms
    results['E_pot'] = uc.set_in_units(thermo.PotEng.mean() / natoms, lammps_units['energy'])
    results['E_pot_std'] = uc.set_in_units(thermo.PotEng.std() / natoms, lammps_units['energy'])

    results['E_total'] = uc.set_in_units(thermo.TotEng.mean() / natoms, lammps_units['energy'])
    results['E_total_std'] = uc.set_in_units(thermo.TotEng.std() / natoms, lammps_units['energy'])

    results['lx'] = uc.set_in_units(thermo.Lx.mean(), lammps_units['length'])
    results['lx_std'] = uc.set_in_units(thermo.Lx.std(), lammps_units['length'])
    results['ly'] = uc.set_in_units(thermo.Ly.mean(), lammps_units['length'])
    results['ly_std'] = uc.set_in_units(thermo.Ly.std(), lammps_units['length'])
    results['lz'] = uc.set_in_units(thermo.Lz.mean(), lammps_units['length'])
    results['lz_std'] = uc.set_in_units(thermo.Lz.std(), lammps_units['length'])
    results['xy'] = uc.set_in_units(thermo.Xy.mean(), lammps_units['length'])
    results['xy_std'] = uc.set_in_units(thermo.Xy.std(), lammps_units['length'])
    results['xz'] = uc.set_in_units(thermo.Xz.mean(), lammps_units['length'])
    results['xz_std'] = uc.set_in_units(thermo.Xz.std(), lammps_units['length'])
    results['yz'] = uc.set_in_units(thermo.Yz.mean(), lammps_units['length'])
    results['yz_std'] = uc.set_in_units(thermo.Yz.std(), lammps_units['length'])

    results['measured_pxx'] = uc.set_in_units(thermo.Pxx.mean(), lammps_units['pressure'])
    results['measured_pxx_std'] = uc.set_in_units(thermo.Pxx.std(), lammps_units['pressure'])
    results['measured_pyy'] = uc.set_in_units(thermo.Pyy.mean(), lammps_units['pressure'])
    results['measured_pyy_std'] = uc.set_in_units(thermo.Pyy.std(), lammps_units['pressure'])
    results['measured_pzz'] = uc.set_in_units(thermo.Pzz.mean(), lammps_units['pressure'])
    results['measured_pzz_std'] = uc.set_in_units(thermo.Pzz.std(), lammps_units['pressure'])
    results['measured_pxy'] = uc.set_in_units(thermo.Pxy.mean(), lammps_units['pressure'])
    results['measured_pxy_std'] = uc.set_in_units(thermo.Pxy.std(), lammps_units['pressure'])
    results['measured_pxz'] = uc.set_in_units(thermo.Pxz.mean(), lammps_units['pressure'])
    results['measured_pxz_std'] = uc.set_in_units(thermo.Pxz.std(), lammps_units['pressure'])
    results['measured_pyz'] = uc.set_in_units(thermo.Pyz.mean(), lammps_units['pressure'])
    results['measured_pyz_std'] = uc.set_in_units(thermo.Pyz.std(), lammps_units['pressure'])
    results['temp'] = thermo.Temp.mean()
    results['temp_std'] = thermo.Temp.std()

    return results

2.2. integrator_info()

[5]:
def integrator_info(integrator: 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,
                    temperature: float = 0.0,
                    velocity_temperature: Optional[float] = None,
                    randomseed: Optional[int] = None,
                    units: str = 'metal') -> str:
    """
    Generates LAMMPS commands for velocity creation and fix integrators.

    Parameters
    ----------
    integrator : str or None, optional
        The integration method to use. Options are 'npt', 'nvt', 'nph',
        'nve', 'nve+l', 'nph+l'. The +l options use Langevin thermostat.
        (Default is None, which will use 'nph+l' for temperature == 0, and
        'npt' otherwise.)
    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).
    temperature : float, optional
        The temperature to relax at (default is 0.0).
    velocity_temperature : float or None, optional
        The temperature to use for generating initial atomic velocities with
        integrators containing thermostats.  If None (default) then it will
        be set to 2 * temperature + 1.  If 0.0 then the velocities will not be
        (re)set.
    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.)
    units : str, optional
        The LAMMPS units style to use (default is 'metal').

    Returns
    -------
    str
        The generated LAMMPS input lines for velocity create and fix
        integration commands.
    """

    # Get lammps units
    lammps_units = lmp.style.unit(units)

    # Convert pressures to lammps units
    p_xx = uc.get_in_units(p_xx, lammps_units['pressure'])
    p_yy = uc.get_in_units(p_yy, lammps_units['pressure'])
    p_zz = uc.get_in_units(p_zz, lammps_units['pressure'])
    p_xy = uc.get_in_units(p_xy, lammps_units['pressure'])
    p_xz = uc.get_in_units(p_xz, lammps_units['pressure'])
    p_yz = uc.get_in_units(p_yz, lammps_units['pressure'])

    # Check temperature and set default integrator
    if temperature == 0.0:
        if integrator is None: integrator = 'nph+l'
        assert integrator not in ['npt', 'nvt'], 'npt and nvt cannot run at 0 K'
    elif temperature > 0:
        if integrator is None: integrator = 'npt'
    else:
        raise ValueError('Temperature must be positive')

    # Set dampening parameters based on timestep values
    temperature_damp = 100 * lmp.style.timestep(units)
    pressure_damp = 1000 * lmp.style.timestep(units)

    # Set default randomseed
    if randomseed is None:
        randomseed = random.randint(1, 900000000)

    # Set default velocity_temperature
    if velocity_temperature is None:
        velocity_temperature = 2.0 * temperature + 1

    # Build info_lines
    info_lines = []

    if integrator == 'npt':

        if velocity_temperature > 0.0:
            info_lines.append(f'velocity all create {velocity_temperature} {randomseed}')

        info_lines.extend([
            f'fix npt all npt temp {temperature} {temperature} {temperature_damp} &',
            f'                x {p_xx} {p_xx} {pressure_damp} &',
            f'                y {p_yy} {p_yy} {pressure_damp} &',
            f'                z {p_zz} {p_zz} {pressure_damp} &',
            f'                xy {p_xy} {p_xy} {pressure_damp} &',
            f'                xz {p_xz} {p_xz} {pressure_damp} &',
            f'                yz {p_yz} {p_yz} {pressure_damp}'
        ])


    elif integrator == 'nvt':

        if velocity_temperature > 0.0:
            info_lines.append(f'velocity all create {velocity_temperature} {randomseed}')

        info_lines.extend([
            f'fix nvt all nvt temp {temperature} {temperature} {temperature_damp}'
        ])

    elif integrator == 'nph':

        info_lines.extend([
            f'fix nph all nph x {p_xx} {p_xx} {pressure_damp} &',
            f'                y {p_yy} {p_yy} {pressure_damp} &',
            f'                z {p_zz} {p_zz} {pressure_damp} &',
            f'                xy {p_xy} {p_xy} {pressure_damp} &',
            f'                xz {p_xz} {p_xz} {pressure_damp} &',
            f'                yz {p_yz} {p_yz} {pressure_damp}'
        ])

    elif integrator == 'nve':

        info_lines.extend([
            'fix nve all nve'
        ])

    elif integrator == 'nve+l':

        if velocity_temperature > 0.0:
            info_lines.append(f'velocity all create {velocity_temperature} {randomseed}')

        info_lines.extend([
            'fix nve all nve',
            f'fix langevin all langevin {temperature} {temperature} {temperature_damp} {randomseed}'
        ])

    elif integrator == 'nph+l':

        info_lines.extend([
            f'fix nph all nph x {p_xx} {p_xx} {pressure_damp} &',
            f'                y {p_yy} {p_yy} {pressure_damp} &',
            f'                z {p_zz} {p_zz} {pressure_damp} &',
            f'                xy {p_xy} {p_xy} {pressure_damp} &',
            f'                xz {p_xz} {p_xz} {pressure_damp} &',
            f'                yz {p_yz} {p_yz} {pressure_damp}',
            f'fix langevin all langevin {temperature} {temperature} {temperature_damp} {randomseed}'
        ])

    else:
        raise ValueError('Invalid integrator style')

    return '\n'.join(info_lines)

2.3. full_relax.template file

[6]:
with open('full_relax.template', 'w') as f:
    f.write("""#LAMMPS input script that performs a simple dynamic integration

box tilt large

<atomman_system_pair_info>

change_box all triclinic

compute pe all pe/atom
compute ke all ke/atom
compute stress all stress/atom <stressterm>

thermo <thermosteps>
thermo_style custom step temp pe ke etotal lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy
thermo_modify format float %.13e
timestep 0.001

<integrator_info>

dump dumpit all custom <dumpsteps> *.dump <dump_keys>
dump_modify dumpit format <dump_modify_format>
restart <restartsteps> *.restart

run <runsteps> upto""")

2.4. full_relax_restart.template file

[7]:
with open('full_relax_restart.template', 'w') as f:
    f.write("""#LAMMPS input script that performs a simple dynamic integration

box tilt large

<atomman_pair_restart_info>

change_box all triclinic

compute pe all pe/atom
compute ke all ke/atom
compute stress all stress/atom <stressterm>

thermo <thermosteps>
thermo_style custom step temp pe ke etotal lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy
thermo_modify format float %.13e

<integrator_restart_info>

dump dumpit all custom <dumpsteps> *.dump <dump_keys>
dump_modify dumpit format <dump_modify_format>
restart <restartsteps> *.restart

run <runsteps> upto""")

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.

[8]:
lammps_command = 'lmp'
mpi_command = 'C:/Program Files/MPICH2/bin/mpiexec.exe -localonly 6'

# Optional: check that LAMMPS works and show its version
print(f'LAMMPS version = {am.lammps.checkversion(lammps_command)["version"]}')
LAMMPS version = 15 Sep 2022

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

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

[10]:
# 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).

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

3.5. Calculation-specific parameters

  • pressure_xx gives the xx component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).

  • pressure_yy gives the yy component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).

  • pressure_zz gives the zz component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).

  • pressure_xy gives the xy component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).

  • pressure_xz gives the xz component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).

  • pressure_yz gives the yz component of the pressure to equilibriate the system to (npt, nph, and nph+l styles).

  • temperature gives the temperature to equilibriate the system to (nvt, npt, nve+l, and nph+l styles).

  • integrator specifies the integrator style to use. Default value is ‘nph+l’ for temperature = 0, and ‘npt’ otherwise.

  • runsteps is the total number of integration timesteps to perform. Default value is 220000.

  • thermosteps specifies to output thermo values every this many timesteps. Default value is 100.

  • dumpsteps specifies to output dump files every this many timesteps. Default value is runsteps (only first and last steps are outputted as dump files).

  • dumpsteps specifies to output restart files every this many timesteps. Default value is runsteps (only last step is outputted as a restart file).

  • equilsteps is the number of timesteps to equilibriate the system for. Only thermo values associated with timesteps greater than equilsteps will be included in the mean and standard deviation calculations. Default value is 20000.

  • randomseed specifies a random number seed used to generate the initial atomic velocities and the Langevin thermostat fluctuations. Default value generates a new random integer every time.

[12]:
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')
temperature = 300.0
integrator = 'npt'
runsteps = 220000
thermosteps = 100
dumpsteps = runsteps
equilsteps = 20000
randomseed = None

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.

[13]:
results_dict = relax_dynamic(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,
                             temperature = temperature,
                             runsteps = runsteps,
                             integrator = integrator,
                             thermosteps = thermosteps,
                             dumpsteps = dumpsteps,
                             equilsteps = equilsteps,
                             randomseed = randomseed)
print(results_dict.keys())
dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'nsamples', 'E_pot', 'E_pot_std', 'E_total', 'E_total_std', 'lx', 'lx_std', 'ly', 'ly_std', 'lz', 'lz_std', 'xy', 'xy_std', 'xz', 'xz_std', 'yz', 'yz_std', 'measured_pxx', 'measured_pxx_std', 'measured_pyy', 'measured_pyy_std', 'measured_pzz', 'measured_pzz_std', 'measured_pxy', 'measured_pxy_std', 'measured_pxz', 'measured_pxz_std', 'measured_pyz', 'measured_pyz_std', 'temp', 'temp_std'])

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.

  • ‘nsamples’ (int) - The number of thermodynamic samples included in the mean and standard deviation estimates. Can also be used to estimate standard error values assuming that the thermo step size is large enough (typically >= 100) to assume the samples to be independent.

  • ‘E_pot’ (float) - The mean measured potential energy.

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

  • ‘temp’ (float) - The mean measured temperature.

  • ‘E_pot_std’ (float) - The standard deviation in the measured potential energy values.

  • ‘measured_pxx_std’ (float) - The standard deviation in the measured x tensile pressure of the relaxed system.

  • ‘measured_pyy_std’ (float) - The standard deviation in the measured y tensile pressure of the relaxed system.

  • ‘measured_pzz_std’ (float) - The standard deviation in the measured z tensile pressure of the relaxed system.

  • ‘measured_pxy_std’ (float) - The standard deviation in the measured xy shear pressure of the relaxed system.

  • ‘measured_pxz_std’ (float) - The standard deviation in the measured xz shear pressure of the relaxed system.

  • ‘measured_pyz_std’ (float) - The standard deviation in the measured yz shear pressure of the relaxed system.

  • ‘temp_std’ (float) - The standard deviation in the measured temperature values.

[14]:
# 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'])
0.dump
('Ni',)
220000.dump
('Ni',)
[15]:
print('number of samples =', results_dict['nsamples'])

# Calculate sqrt(N) for transforming standard deviations to standard errors
sqrt_n = results_dict['nsamples']**0.5
number of samples = 2001
[16]:
energy_unit = 'eV'

print('Per-atom potential energy and standard error:')
print('E_pot =', uc.get_in_units(results_dict['E_pot'], energy_unit),
      '+-', uc.get_in_units(results_dict['E_pot_std'], energy_unit) / sqrt_n, energy_unit)
Per-atom potential energy and standard error:
E_pot = -4.413036859648026 +- 1.1020441612326114e-05 eV
[17]:
length_unit = 'angstrom'

print('Box lengths, tilts and standard errors:')
print('lx =', uc.get_in_units(results_dict['lx'] / sizemults[0], length_unit),
      '+-', uc.get_in_units(results_dict['lx_std'] / sizemults[0] / sqrt_n, length_unit), length_unit)
print('ly =', uc.get_in_units(results_dict['ly'] / sizemults[1], length_unit),
      '+-', uc.get_in_units(results_dict['ly_std'] / sizemults[1] / sqrt_n, length_unit), length_unit)
print('lz =', uc.get_in_units(results_dict['lz'] / sizemults[2], length_unit),
      '+-', uc.get_in_units(results_dict['lz_std'] / sizemults[2] / sqrt_n, length_unit), length_unit)
print('xy =', uc.get_in_units(results_dict['xy'] / sizemults[1], length_unit),
      '+-', uc.get_in_units(results_dict['xy_std'] / sizemults[1] / sqrt_n, length_unit), length_unit)
print('xz =', uc.get_in_units(results_dict['xz'] / sizemults[2], length_unit),
      '+-', uc.get_in_units(results_dict['xz_std'] / sizemults[2] / sqrt_n, length_unit), length_unit)
print('yz =', uc.get_in_units(results_dict['yz'] / sizemults[2], length_unit),
      '+-', uc.get_in_units(results_dict['yz_std'] / sizemults[2] / sqrt_n, length_unit), length_unit)
Box lengths, tilts and standard errors:
lx = 3.5331129564611565 +- 5.773669391414223e-05 angstrom
ly = 3.53296915195529 +- 5.979159319290978e-05 angstrom
lz = 3.5329970355875737 +- 6.0218416870231855e-05 angstrom
xy = -2.312045931302908e-05 +- 6.654549358172896e-05 angstrom
xz = 2.8157721680541136e-06 +- 6.374138106096007e-05 angstrom
yz = 6.608261957953486e-05 +- 5.655036762675057e-05 angstrom
[18]:
pressure_unit = 'GPa'

# Show the computed pressure tensor
print('Pxx =', uc.get_in_units(results_dict['measured_pxx'], pressure_unit),
      '+-', uc.get_in_units(results_dict['measured_pxx_std'], pressure_unit) / sqrt_n, pressure_unit)
print('Pyy =', uc.get_in_units(results_dict['measured_pyy'], pressure_unit),
      '+-', uc.get_in_units(results_dict['measured_pyy_std'], pressure_unit) / sqrt_n, pressure_unit)
print('Pzz =', uc.get_in_units(results_dict['measured_pzz'], pressure_unit),
      '+-', uc.get_in_units(results_dict['measured_pzz_std'], pressure_unit) / sqrt_n, pressure_unit)
print('Pyz =', uc.get_in_units(results_dict['measured_pyz'], pressure_unit),
      '+-', uc.get_in_units(results_dict['measured_pyz_std'], pressure_unit) / sqrt_n, pressure_unit)
print('Pxz =', uc.get_in_units(results_dict['measured_pxz'], pressure_unit),
      '+-', uc.get_in_units(results_dict['measured_pxz_std'], pressure_unit) / sqrt_n, pressure_unit)
print('Pxy =', uc.get_in_units(results_dict['measured_pxy'], pressure_unit),
      '+-', uc.get_in_units(results_dict['measured_pxy_std'], pressure_unit) / sqrt_n, pressure_unit)
Pxx = -0.0006119571960764553 +- 0.0038901984574203647 GPa
Pyy = -1.7729015128904964e-05 +- 0.003945718913997667 GPa
Pzz = -0.0010041150769341074 +- 0.004009725790232544 GPa
Pyz = -0.0005549426289240447 +- 0.0024020374408581903 GPa
Pxz = -0.00071677380151718 +- 0.002710929768537196 GPa
Pxy = -0.00011167112626227435 +- 0.002844613710616187 GPa