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.

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.

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 datetime
import random
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_dynamic'

# 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_mpi'
#mpi_command = None
mpi_command = 'C:/Program Files/MPICH2/bin/mpiexec.exe -localonly 6'

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 = [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.5. Specify calculation-specific run 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).

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

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

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

3.1. full_relax.template

[8]:
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 <runsteps> *.restart

run <runsteps> upto""")

3.2. integrator_info()

[9]:
def integrator_info(integrator=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, temperature=0.0, randomseed=None,
                    units='metal'):
    """
    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).
    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)
    Px = uc.get_in_units(p_xx, lammps_units['pressure'])
    Py = uc.get_in_units(p_yy, lammps_units['pressure'])
    Pz = uc.get_in_units(p_zz, lammps_units['pressure'])
    Pxy = uc.get_in_units(p_xy, lammps_units['pressure'])
    Pxz = uc.get_in_units(p_xz, lammps_units['pressure'])
    Pyz = uc.get_in_units(p_yz, lammps_units['pressure'])
    T = temperature

    # 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 default randomseed
    if randomseed is None: randomseed = random.randint(1, 900000000)

    if integrator == 'npt':
        start_temp = T*2.+1
        Tdamp = 100 * lmp.style.timestep(units)
        Pdamp = 1000 * lmp.style.timestep(units)
        int_info = '\n'.join([
                'velocity all create %f %i' % (start_temp, randomseed),
                'fix npt all npt temp %f %f %f &' % (T, T, Tdamp),
                '                x %f %f %f &' % (Px, Px, Pdamp),
                '                y %f %f %f &' % (Py, Py, Pdamp),
                '                z %f %f %f &' % (Pz, Pz, Pdamp),
                '                xy %f %f %f &' % (Pxy, Pxy, Pdamp),
                '                xz %f %f %f &' % (Pxz, Pxz, Pdamp),
                '                yz %f %f %f' % (Pyz, Pyz, Pdamp),
                ])

    elif integrator == 'nvt':
        start_temp = T*2.+1
        Tdamp = 100 * lmp.style.timestep(units)
        int_info = '\n'.join([
                'velocity all create %f %i' % (start_temp, randomseed),
                'fix nvt all nvt temp %f %f %f' % (T, T, Tdamp),
                ])

    elif integrator == 'nph':
        Pdamp = 1000 * lmp.style.timestep(units)
        int_info = '\n'.join([
                'fix nph all nph x %f %f %f &' % (Px, Px, Pdamp),
                '                y %f %f %f &' % (Py, Py, Pdamp),
                '                z %f %f %f &' % (Pz, Pz, Pdamp),
                '                xy %f %f %f &' % (Pxy, Pxy, Pdamp),
                '                xz %f %f %f &' % (Pxz, Pxz, Pdamp),
                '                yz %f %f %f' % (Pyz, Pyz, Pdamp),
                ])

    elif integrator == 'nve':
        int_info = 'fix nve all nve'

    elif integrator == 'nve+l':
        start_temp = T*2.+1
        Tdamp = 100 * lmp.style.timestep(units)
        int_info = '\n'.join([
                'velocity all create %f %i' % (start_temp, randomseed),
                'fix nve all nve',
                'fix langevin all langevin %f %f %f %i' % (T, T, Tdamp,
                                                           randomseed),
                ])

    elif integrator == 'nph+l':
        start_temp = T*2.+1
        Tdamp = 100 * lmp.style.timestep(units)
        Pdamp = 1000 * lmp.style.timestep(units)
        int_info = '\n'.join([
                'fix nph all nph x %f %f %f &' % (Px, Px, Pdamp),
                '                y %f %f %f &' % (Py, Py, Pdamp),
                '                z %f %f %f &' % (Pz, Pz, Pdamp),
                '                xy %f %f %f &' % (Pxy, Pxy, Pdamp),
                '                xz %f %f %f &' % (Pxz, Pxz, Pdamp),
                '                yz %f %f %f' % (Pyz, Pyz, Pdamp),
                'fix langevin all langevin %f %f %f %i' % (T, T, Tdamp,
                                                           randomseed),
                ])
    else:
        raise ValueError('Invalid integrator style')

    return int_info

3.3. full_relax()

[10]:
def relax_dynamic(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,
                  temperature=0.0, integrator=None, runsteps=220000,
                  thermosteps=100, dumpsteps=None, equilsteps=20000,
                  randomseed=None):
    """
    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).
    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:

        - **'relaxed_system'** (*float*) - The relaxed system.
        - **'E_coh'** (*float*) - The mean measured cohesive 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_coh_std'** (*float*) - The standard deviation in the measured
          cohesive 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.
    """
    # 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']

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

    # 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

    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
    lammps_variables['thermosteps'] = thermosteps
    lammps_variables['runsteps'] = runsteps
    lammps_variables['dumpsteps'] = dumpsteps

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

    # Write lammps input script
    template_file = 'full_relax.template'
    lammps_script = 'full_relax.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
    output = lmp.run(lammps_command, lammps_script, mpi_command)

    # Extract LAMMPS thermo data.
    results = {}
    thermo = output.simulations[0]['thermo']

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

    # Load relaxed system from dump file
    last_dump_file = str(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_coh'] = uc.set_in_units(thermo.PotEng.mean() / natoms, lammps_units['energy'])
    results['E_coh_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

4. Run calculation function(s)

[11]:
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)
[12]:
results_dict.keys()
[12]:
dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'nsamples', 'E_coh', 'E_coh_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'])

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.

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

5.2. Display number of samples

[14]:
print('# of samples = ', results_dict['nsamples'])
# of samples =  2001

5.3. Print mean and standard deviations of measured values

[15]:
print('Measured thermo data with standard deviation:')
print('Ecoh = %9f +- %9f %s' % (uc.get_in_units(results_dict['E_coh'], energy_unit),
                                    uc.get_in_units(results_dict['E_coh_std'], energy_unit),
                                    energy_unit))
print('lx   = %9f +- %9f %s' % (uc.get_in_units(results_dict['lx'], length_unit),
                                    uc.get_in_units(results_dict['lx_std'], length_unit),
                                    length_unit))
print('ly   = %9f +- %9f %s' % (uc.get_in_units(results_dict['ly'], length_unit),
                                    uc.get_in_units(results_dict['ly_std'], length_unit),
                                    length_unit))
print('lz   = %9f +- %9f %s' % (uc.get_in_units(results_dict['lz'], length_unit),
                                    uc.get_in_units(results_dict['lz_std'], length_unit),
                                    length_unit))
print('xy   = %9f +- %9f %s' % (uc.get_in_units(results_dict['xy'], length_unit),
                                    uc.get_in_units(results_dict['xy_std'], length_unit),
                                    length_unit))
print('xz   = %9f +- %9f %s' % (uc.get_in_units(results_dict['xz'], length_unit),
                                    uc.get_in_units(results_dict['xz_std'], length_unit),
                                    length_unit))
print('yz   = %9f +- %9f %s' % (uc.get_in_units(results_dict['yz'], length_unit),
                                    uc.get_in_units(results_dict['yz_std'], length_unit),
                                    length_unit))
print('T    = %9f +- %9f %s' % (results_dict['temp'],results_dict['temp_std'], 'K'))
print('Pxx  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pxx'], pressure_unit),
                                    uc.get_in_units(results_dict['measured_pxx_std'], pressure_unit),
                                    pressure_unit))
print('Pyy  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pyy'], pressure_unit),
                                    uc.get_in_units(results_dict['measured_pyy_std'], pressure_unit),
                                    pressure_unit))
print('Pzz  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pzz'], pressure_unit),
                                    uc.get_in_units(results_dict['measured_pzz_std'], pressure_unit),
                                    pressure_unit))
print('Pxy  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pxy'], pressure_unit),
                                    uc.get_in_units(results_dict['measured_pxy_std'], pressure_unit),
                                    pressure_unit))
print('Pxz  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pxz'], pressure_unit),
                                    uc.get_in_units(results_dict['measured_pxz_std'], pressure_unit),
                                    pressure_unit))
print('Pyz  = %9f +- %9f %s' % (uc.get_in_units(results_dict['measured_pyz'], pressure_unit),
                                    uc.get_in_units(results_dict['measured_pyz_std'], pressure_unit),
                                    pressure_unit))
Measured thermo data with standard deviation:
Ecoh = -4.413022 +-  0.000507 eV
lx   = 35.330025 +-  0.025861 angstrom
ly   = 35.329559 +-  0.026324 angstrom
lz   = 35.330773 +-  0.026982 angstrom
xy   = -0.000836 +-  0.029263 angstrom
xz   = -0.000028 +-  0.027548 angstrom
yz   =  0.000586 +-  0.024387 angstrom
T    = 299.963948 +-  3.788839 K
Pxx  =  0.000186 +-  0.190085 GPa
Pyy  =  0.000212 +-  0.188994 GPa
Pzz  =  0.000974 +-  0.190964 GPa
Pxy  =  0.000299 +-  0.125513 GPa
Pxz  =  0.000353 +-  0.118209 GPa
Pyz  =  0.000150 +-  0.103573 GPa

5.4. Compute lattice constant with standard error

NOTE: This step makes two assumptions

  1. The crystal structure is cubic and remains cubic after relaxation. Check values above to verify this.

  2. The thermosteps parameter is large enough that the measurements are not correlated. If thermosteps ≥ 100 this is likely a sound assumption.

[16]:
a = results_dict['lx'] / sizemults[0]
b = results_dict['ly'] / sizemults[1]
c = results_dict['lz'] / sizemults[2]

a_std = results_dict['lx_std'] / sizemults[1]
b_std = results_dict['ly_std'] / sizemults[2]
c_std = results_dict['lz_std'] / sizemults[2]

a_mean = (a + b + c) / 3
a_combined_std = ((a_std**2 + b_std**2 + c_std**2
                   + (a - a_mean)**2 + (b - a_mean)**2 + (c - a_mean)**2) / 3)**0.5
a_standard_error = a_combined_std * (3 * results_dict['nsamples'])**-0.5

print('Cubic lattice constant with standard error:')
print('a = %9f +- %9f %s' % (uc.get_in_units(a_mean, length_unit),
                             uc.get_in_units(a_standard_error, length_unit),
                             length_unit))
Cubic lattice constant with standard error:
a =  3.533012 +-  0.000034 angstrom

5.5. Load final configuration and show box

[17]:
finalsystem = am.load('atom_dump', results_dict['dumpfile_final'],
                      symbols=results_dict['symbols_final'])
print(finalsystem.box)
avect =  [35.319,  0.000,  0.000]
bvect =  [-0.027, 35.360,  0.000]
cvect =  [ 0.026, -0.004, 35.285]
origin = [-0.159, -0.180, -0.143]