elastic_constants_static calculation style

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

Introduction

The elastic_constants_static calculation style computes the elastic constants, \(C_{ij}\), for a system by applying small strains and performing static energy minimizations of the initial and strained configurations. Three estimates of the elastic constants are returned: one for applying positive strains, one for applying negative strains, and a normalized estimate that averages the ± strains and the symmetric components of the \(C_{ij}\) tensor.

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

  • Unlike the previous LAMMPS_ELASTIC calculation, this calculation does not perform a box relaxation on the system prior to evaluating the elastic constants. This allows for the static elastic constants to be evaluated for systems that are relaxed to different pressures.

  • The elastic constants are estimated using small strains. Depending on the potential, the values for the elastic constants may vary with the size of the strain. This can come about either if the strain exceeds the linear elastic regime.

  • Some classical interatomic potentials have discontinuities in the fourth derivative of the energy function with respect to position. If the strained states straddle one of these discontinuities the resulting static elastic constants values will be nonsense.

Method and Theory

The calculation method used here for computing elastic constants is based on the method used in the ELASTIC demonstration script created by Steve Plimpton and distributed with LAMMPS.

The math in this section uses Voigt notation, where indicies i,j correspond to 1=xx, 2=yy, 3=zz, 4=yz, 5=xz, and 6=xy, and x, y and z are orthogonal box vectors.

A LAMMPS simulation performs thirteen energy/force minimizations

  • One for relaxing the initial system.

  • Twelve for relaxing systems in which a small strain of magnitude \(\Delta \epsilon\) is applied to the system in both the positive and negative directions of the six Voigt strain components, \(\pm \Delta \epsilon_{i}\).

The system virial pressures, \(P_{i}\), are recorded for each of the thirteen relaxed configurations. Two estimates for the \(C_{ij}\) matrix for the system are obtained as

\[C_{ij}^+ = - \frac{P_i(\Delta \epsilon_j) - P_i(0)}{\Delta \epsilon},\]
\[C_{ij}^- = - \frac{P_i(0) - P_i(-\Delta \epsilon_j)}{\Delta \epsilon}.\]

The negative out front comes from the fact that the system-wide stress state is \(\sigma_i = -P_i\). A normalized, average estimate is also obtained by averaging the positive and negative strain estimates, as well as the symmetric components of the tensor

\[C_{ij} = \frac{C_{ij}^+ + C_{ij}^- + C_{ji}^+ + C_{ji}^-}{4}.\]

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
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 = 'elastic_constants_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
with open('C:/Users/lmh1/Documents/potentials_nist_gov/password.txt') as f:
    user, pswd = f.read().strip().split()
library = iprPy.Library(load=['lammps_potentials'], username=user, password=pswd)

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 loaded from the database for the prototype.

[5]:
# Create ucell by loading prototype record
ucell = am.load('crystal', potential=potential, family='A1--Cu--fcc', database=library)

print(ucell)
avect =  [ 3.520,  0.000,  0.000]
bvect =  [ 0.000,  3.520,  0.000]
cvect =  [ 0.000,  0.000,  3.520]
origin = [ 0.000,  0.000,  0.000]
natoms = 4
natypes = 1
symbols = ('Ni',)
pbc = [ True  True  True]
per-atom properties = ['atype', 'pos']
     id   atype  pos[0]  pos[1]  pos[2]
      0       1   0.000   0.000   0.000
      1       1   0.000   1.760   1.760
      2       1   1.760   0.000   1.760
      3       1   1.760   1.760   0.000

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

  • strainrange specifies the \(\Delta \epsilon\) strain range to use in estimating \(C_{ij}\).

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

[7]:
strainrange = 1e-7
energytolerance = 1e-8
forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')
maxiterations = 10000
maxevaluations = 100000
maxatommotion = uc.set_in_units(0.01, 'angstrom')

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

3.1. cij.template

[8]:
with open('cij.template', 'w') as f:
    f.write("""# Performs simulations to statically evaluate elastic constants using small strains
# Based on the LAMMPS_ELASTIC script by Aidan Thompson (Sandia, athomps@sandia.gov)

box tilt large

<atomman_system_info>

change_box all triclinic

# Specify strain
variable strain equal <strainrange>

# Define minimization parameters
variable etol equal <etol>
variable ftol equal <ftol>
variable maxiter equal <maxiter>
variable maxeval equal <maxeval>
variable dmax equal <dmax>

# Specify variables of the initial configuration's dimensions
variable lx0 equal $(lx)
variable ly0 equal $(ly)
variable lz0 equal $(lz)

# Specify the thermo properties to calculate
variable peatom equal pe/atoms

# Read in potential and thermo information
include potential.in

# Relax initial configuration and save as restart
minimize ${etol} ${ftol} ${maxiter} ${maxeval}
write_restart initial.restart

# Apply -xx strain
clear
read_restart initial.restart
include potential.in

variable delta equal -${strain}*${lx0}
change_box all x delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +xx strain
clear
read_restart initial.restart
include potential.in

variable delta equal ${strain}*${lx0}
change_box all x delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -yy strain
clear
read_restart initial.restart
include potential.in

variable delta equal -${strain}*${ly0}
change_box all y delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +yy strain
clear
read_restart initial.restart
include potential.in

variable delta equal ${strain}*${ly0}
change_box all y delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -zz strain
clear
read_restart initial.restart
include potential.in

variable delta equal -${strain}*${lz0}
change_box all z delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +zz strain
clear
read_restart initial.restart
include potential.in

variable delta equal ${strain}*${lz0}
change_box all z delta 0 ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -yz strain
clear
read_restart initial.restart
include potential.in

variable delta equal -${strain}*${lz0}
change_box all yz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +yz strain
clear
read_restart initial.restart
include potential.in

variable delta equal ${strain}*${lz0}
change_box all yz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -xz strain
clear
read_restart initial.restart
include potential.in

variable delta equal -${strain}*${lz0}
change_box all xz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +xz strain
clear
read_restart initial.restart
include potential.in

variable delta equal ${strain}*${lz0}
change_box all xz delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply -xy strain
clear
read_restart initial.restart
include potential.in

variable delta equal -${strain}*${ly0}
change_box all xy delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}

# Apply +xy strain
clear
read_restart initial.restart
include potential.in

variable delta equal ${strain}*${ly0}
change_box all xy delta ${delta} remap units box
minimize ${etol} ${ftol} ${maxiter} ${maxeval}
""")

3.2. potential.template

[9]:
with open('potential.template', 'w') as f:
    f.write("""# NOTE: This script can be modified for different pair styles
# See in.elastic for more info.

# Choose potential
<atomman_pair_info>

# Setup minimization style
min_modify dmax ${dmax}

# Setup output
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e
""")

3.3. elastic_constants_static()

[10]:
def elastic_constants_static(lammps_command, system, potential, mpi_command=None,
                             strainrange=1e-6, etol=0.0, ftol=0.0, maxiter=10000,
                             maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom')):
    """
    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.
    strainrange : float, optional
        The small strain value to apply when calculating the elastic
        constants (default is 1e-6).
    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).

    Returns
    -------
    dict
        Dictionary of results consisting of keys:

        - **'a_lat'** (*float*) - The relaxed a lattice constant.
        - **'b_lat'** (*float*) - The relaxed b lattice constant.
        - **'c_lat'** (*float*) - The relaxed c lattice constant.
        - **'alpha_lat'** (*float*) - The alpha lattice angle.
        - **'beta_lat'** (*float*) - The beta lattice angle.
        - **'gamma_lat'** (*float*) - The gamma lattice angle.
        - **'E_coh'** (*float*) - The cohesive energy of the relaxed system.
        - **'stress'** (*numpy.array*) - The measured stress state of the
          relaxed system.
        - **'C_elastic'** (*atomman.ElasticConstants*) - The relaxed system's
          elastic constants.
        - **'system_relaxed'** (*atomman.System*) - 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 = {}

    # Convert hexagonal cells to orthorhombic to avoid LAMMPS tilt issues
    if am.tools.ishexagonal(system.box):
        system = system.rotate([[2,-1,-1,0], [0, 1, -1, 0], [0,0,0,1]])

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

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

    # Define lammps variables
    lammps_variables = {}
    system_info = system.dump('atom_data', f='init.dat',
                              units=potential.units,
                              atom_style=potential.atom_style)
    lammps_variables['atomman_system_info'] = system_info
    lammps_variables['atomman_pair_info'] = potential.pair_info(system.symbols)
    lammps_variables['strainrange'] = strainrange
    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'])

    # Fill in template files
    template_file = 'cij.template'
    lammps_script = 'cij.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, '<', '>'))

    template_file2 = 'potential.template'
    lammps_script2 = 'potential.in'
    template = iprPy.tools.read_calc_file(template_file2, filedict)
    with open(lammps_script2, 'w') as f:
        f.write(iprPy.tools.filltemplate(template, lammps_variables, '<', '>'))

    # Run LAMMPS
    output = lmp.run(lammps_command, lammps_script, mpi_command)

    # Pull out initial state
    thermo = output.simulations[0]['thermo']
    pxx0 = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
    pyy0 = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
    pzz0 = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
    pyz0 = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
    pxz0 = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
    pxy0 = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])

    # Negative strains
    cij_n = np.empty((6,6))
    for i in range(6):
        j = 1 + i * 2
        # Pull out strained state
        thermo = output.simulations[j]['thermo']
        pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
        pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
        pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
        pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
        pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
        pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])

        # Calculate cij_n using stress changes
        cij_n[i] = np.array([pxx - pxx0, pyy - pyy0, pzz - pzz0,
                             pyz - pyz0, pxz - pxz0, pxy - pxy0]) / strainrange

    # Positive strains
    cij_p = np.empty((6,6))
    for i in range(6):
        j = 2 + i * 2
        # Pull out strained state
        thermo = output.simulations[j]['thermo']
        pxx = uc.set_in_units(thermo.Pxx.values[-1], lammps_units['pressure'])
        pyy = uc.set_in_units(thermo.Pyy.values[-1], lammps_units['pressure'])
        pzz = uc.set_in_units(thermo.Pzz.values[-1], lammps_units['pressure'])
        pyz = uc.set_in_units(thermo.Pyz.values[-1], lammps_units['pressure'])
        pxz = uc.set_in_units(thermo.Pxz.values[-1], lammps_units['pressure'])
        pxy = uc.set_in_units(thermo.Pxy.values[-1], lammps_units['pressure'])

        # Calculate cij_p using stress changes
        cij_p[i] = np.array([pxx - pxx0, pyy - pyy0, pzz - pzz0,
                              pyz - pyz0, pxz - pxz0, pxy - pxy0]) / -strainrange

    # Average symmetric values
    cij = (cij_n + cij_p) / 2
    for i in range(6):
        for j in range(i):
            cij[i,j] = cij[j,i] = (cij[i,j] + cij[j,i]) / 2

    # Define results_dict
    results_dict = {}
    results_dict['raw_cij_negative'] = cij_n
    results_dict['raw_cij_positive'] = cij_p
    results_dict['C'] = am.ElasticConstants(Cij=cij)

    return results_dict

4. Run calculation function(s)

[11]:
results_dict = elastic_constants_static(lammps_command, system, potential,
                                        mpi_command = mpi_command,
                                        strainrange = strainrange,
                                        etol = energytolerance,
                                        ftol = forcetolerance,
                                        maxiter = maxiterations,
                                        maxeval = maxevaluations,
                                        dmax = maxatommotion)
[12]:
results_dict.keys()
[12]:
dict_keys(['raw_cij_negative', 'raw_cij_positive', 'C'])

5. Report results

5.1. Define units for outputting values

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

[13]:
pressure_unit = 'GPa'

5.2. Print raw Cij values for \(\pm \Delta \epsilon\) states

[14]:
print('Raw Cij for negative strains ('+pressure_unit+') =')
for Ci in uc.get_in_units(results_dict['raw_cij_negative'], pressure_unit):
    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))
print()

print('Raw Cij for positive strains ('+pressure_unit+') =')
for Ci in uc.get_in_units(results_dict['raw_cij_positive'], pressure_unit):
    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))
Raw Cij for negative strains (GPa) =
[ 247.8622  147.8285  147.8285    0.0000    0.0000   -0.0000]
[ 147.8285  247.8622  147.8285    0.0000    0.0000   -0.0000]
[ 147.8285  147.8285  247.8622   -0.0000    0.0000   -0.0000]
[  -0.0000    0.0001    0.0001  124.8381   -0.0000   -0.0000]
[   0.0001   -0.0000    0.0001   -0.0000  124.8381   -0.0000]
[   0.0001    0.0001   -0.0000   -0.0000   -0.0000  124.8381]

Raw Cij for positive strains (GPa) =
[ 247.8625  147.8283  147.8283   -0.0000    0.0000    0.0000]
[ 147.8283  247.8625  147.8283    0.0000   -0.0000    0.0000]
[ 147.8283  147.8283  247.8625    0.0000   -0.0000    0.0000]
[   0.0000   -0.0001   -0.0001  124.8381    0.0000    0.0000]
[  -0.0001    0.0000   -0.0001    0.0000  124.8381    0.0000]
[  -0.0001   -0.0001    0.0000    0.0000    0.0000  124.8381]

5.3. Print averaged and refined Cij values

[15]:
print('Cij ('+pressure_unit+') =')
for Ci in uc.get_in_units(results_dict['C'].Cij, pressure_unit):
    print('[%9.4f %9.4f %9.4f %9.4f %9.4f %9.4f]' % tuple(Ci))
Cij (GPa) =
[ 247.8624  147.8284  147.8284    0.0000    0.0000    0.0000]
[ 147.8284  247.8624  147.8284    0.0000    0.0000    0.0000]
[ 147.8284  147.8284  247.8624    0.0000    0.0000    0.0000]
[   0.0000    0.0000    0.0000  124.8381    0.0000    0.0000]
[   0.0000    0.0000    0.0000    0.0000  124.8381    0.0000]
[   0.0000    0.0000    0.0000    0.0000    0.0000  124.8381]