relax_box - Methodology and code

Python imports

[1]:
# Standard library imports
from pathlib import Path
import datetime
from copy import deepcopy
from math import floor
from typing import Optional

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

1.2. Display calculation description and theory

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

relax_box calculation style

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

Introduction

The relax_box calculation style refines the lattice parameters of an orthogonal system (crystal structure) by relaxing the box dimensions towards a given pressure. In refining the lattice parameter values, the box dimensions are allowed to relax, but the relative positions of the atoms within the box are held fixed.

This calculations provides a quick tool for obtaining lattice parameters for ideal crystal structures.

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. Method reworked to better treat triclinic systems.

Additional dependencies
Disclaimers
  • NIST disclaimers

  • With this method there is no guarantee that the resulting parameters are for a stable structure. Allowing internal relaxations may result in different values for some structures. Additionally, some transformation paths may be restricted from occurring due to symmetry, i.e. initially cubic structures may remain cubic instead of relaxing to a non-cubic structure.

Method and Theory

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.

An initial system (and corresponding unit cell system) is supplied with box dimensions, \(a_i^0\), close to the equilibrium values. A LAMMPS simulation is performed that evaluates the system’s pressures, \(P_{i}\), for the initial system as given, and subjected to twelve different strain states corresponding to one of \(\epsilon_{i}\) being given a value of \(\frac{\Delta \epsilon}{2}\), where \(\Delta \epsilon\) is the strain range parameter. Using the \(P_{i}\) values obtained from the strained states, the \(C_{ij}\) matrix for the system is estimated as

\[C_{ij} \approx - \frac{P_i(\epsilon_j=\frac{\Delta \epsilon}{2}) - P_i(\epsilon_j=-\frac{\Delta \epsilon}{2})}{\Delta \epsilon}.\]

The negative out front comes from the fact that the system-wide stress state is \(\sigma_i = -P_i\). Using \(C_{ij}\), an attempt is made to compute the elastic compliance matrix as \(S_{ij} = C_{ij}^{-1}\). If successful, new box dimensions are estimated using \(S_{ij}\), \(a_i^0\), and \(P_i\) for the unstrained system

\[a_i = \frac{a_i^0}{1 - (\sum_{j=1}^3{S_{ij} P_j})}.\]

The system is updated using the new box dimensions. The process is repeated until either \(a_i\) converge less than a specified tolerance, \(a_i\) diverge from \(a_i^0\) greater than some limit, or convergence is not reached after 100 iterations. If the calculation is successful, the final \(a_i\) dimensions are reported.

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

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_box(lammps_command: str,
              system: am.System,
              potential: lmp.Potential,
              mpi_command: Optional[str] = None,
              strainrange: float = 1e-6,
              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,
              tol: float = 1e-10,
              diverge_scale: float = 3.0)  -> dict:
    """
    Quickly refines static orthorhombic system by evaluating the elastic
    constants and the virial pressure.

    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).
    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).
    tol : float, optional
        The relative tolerance used to determine if the lattice constants have
        converged (default is 1e-10).
    diverge_scale : float, optional
        Factor to identify if the system's dimensions have diverged.  Divergence
        is identified if either any current box dimension is greater than the
        original dimension multiplied by diverge_scale, or if any current box
        dimension is less than the original dimension divided by diverge_scale.
        (Default is 3.0).

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

        - **'dumpfile_initial'** (*str*) - The name of the initial dump file
          created.
        - **'symbols_initial'** (*list*) - The symbols associated with the
          initial dump file.
        - **'dumpfile_final'** (*str*) - The name of the final dump file
          created.
        - **'symbols_final'** (*list*) - The symbols associated with the final
          dump file.
        - **'lx'** (*float*) - The relaxed lx box length.
        - **'ly'** (*float*) - The relaxed ly box length.
        - **'lz'** (*float*) - The relaxed lz box length.
        - **'xy'** (*float*) - The relaxed xy box tilt.
        - **'xz'** (*float*) - The relaxed xz box tilt.
        - **'yz'** (*float*) - The relaxed yz box tilt.
        - **'E_pot'** (*float*) - The potential energy per atom for the final
          configuration.
        - **'measured_pxx'** (*float*) - The measured x tensile pressure 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.

    Raises
    ------
    RuntimeError
        If system diverges or no convergence reached after 100 cycles.
    """

    # Flag for if values have converged
    converged = False

    # Define current and old systems
    system_current = deepcopy(system)
    system_old = None

    system.dump('atom_dump', f='initial.dump')

    for cycle in range(100):

        # Run LAMMPS and evaluate results based on system_old
        results = cij_run0(lammps_command, system_current, potential,
                           mpi_command=mpi_command, strainrange=strainrange,
                           cycle=cycle)
        pij = results['pij']
        Cij = results['C'].Cij
        system_new = update_box(system_current, results['C'], results['pij'],
                                p_xx, p_yy, p_zz, p_xy, p_xz, p_yz, tol)

        # Compare new and current to test for convergence
        if np.allclose(system_new.box.vects,
                       system_current.box.vects,
                       rtol=tol, atol=0):
            converged = True
            break

        # Compare old and new to test for double-value convergence
        elif system_old is not None and np.allclose(system_new.box.vects,
                                                    system_old.box.vects,
                                                    rtol=tol, atol=0):

            # Update current to average of old and new
            system_current.box_set(a = (system_new.box.a+system_old.box.a) / 2.,
                                   b = (system_new.box.b+system_old.box.b) / 2.,
                                   c = (system_new.box.c+system_old.box.c) / 2.,
                                   scale=True)

            # Calculate Cij for the averaged system
            results = cij_run0(lammps_command, system_current, potential,
                               mpi_command=mpi_command, strainrange=strainrange,
                               cycle=cycle)
            system_new = update_box(system_current, results['C'], results['pij'],
                                    p_xx, p_yy, p_zz, p_xy, p_xz, p_yz, tol)
            converged = True
            break

        # Test for divergence
        elif system_new.box.a < system.box.a / diverge_scale:
            raise RuntimeError('Divergence of box dimensions')
        elif system_new.box.a > system.box.a * diverge_scale:
            raise RuntimeError('Divergence of box dimensions')
        elif system_new.box.b < system.box.b / diverge_scale:
            raise RuntimeError('Divergence of box dimensions')
        elif system_new.box.b > system.box.b * diverge_scale:
            raise RuntimeError('Divergence of box dimensions')
        elif system_new.box.c < system.box.c / diverge_scale:
            raise RuntimeError('Divergence of box dimensions')
        elif system_new.box.c > system.box.c * diverge_scale:
            raise RuntimeError('Divergence of box dimensions')
        elif results['E_pot'] == 0.0:
            raise RuntimeError('Divergence: potential energy is 0')

        # If not converged or diverged, current -> old and new -> current
        else:
            system_old, system_current = system_current, system_new

    # Return values when converged
    if converged:
        system_new.dump('atom_dump', f='final.dump')

        # Build results_dict
        results_dict = {}
        results_dict['dumpfile_initial'] = 'initial.dump'
        results_dict['symbols_initial'] = system.symbols
        results_dict['dumpfile_final'] = 'final.dump'
        results_dict['symbols_final'] = system.symbols

        results_dict['lx'] = system_new.box.lx
        results_dict['ly'] = system_new.box.ly
        results_dict['lz'] = system_new.box.lz
        results_dict['xy'] = system_new.box.xy
        results_dict['xz'] = system_new.box.xz
        results_dict['yz'] = system_new.box.yz

        results_dict['E_pot'] = results['E_pot']
        results_dict['measured_pxx'] = results['pij'][0,0]
        results_dict['measured_pyy'] = results['pij'][1,1]
        results_dict['measured_pzz'] = results['pij'][2,2]
        results_dict['measured_pxy'] = results['pij'][0,1]
        results_dict['measured_pxz'] = results['pij'][0,2]
        results_dict['measured_pyz'] = results['pij'][1,2]

        return results_dict
    else:
        raise RuntimeError('Failed to converge after 100 cycles')

2.2. cij_run0()

[5]:
def cij_run0(lammps_command: str,
             system: am.System,
             potential: lmp.Potential,
             mpi_command: Optional[str] = None,
             strainrange: float = 1e-6,
             cycle: int = 0) -> dict:
    """
    Runs cij_run0.in LAMMPS script to evaluate the elastic constants,
    pressure and potential energy of the current system.

    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).
    cycle : int, optional
        Indicates the iteration cycle of quick_a_Cij().  This is used to
        uniquely save the LAMMPS input and output files.

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

        - **'E_pot'** (*float*) - The potential energy per atom for the
          supplied system.
        - **'pressure'** (*numpy.array*) - The measured pressure state of the
          supplied system.
        - **'C_elastic'** (*atomman.ElasticConstants*) - The supplied system's
          elastic constants.
    """

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

    # Define lammps variables
    lammps_variables = {}
    system_info = system.dump('atom_data', f='init.dat',
                              potential=potential)
    restart_info = potential.pair_restart_info('initial.restart', system.symbols)
    lammps_variables['pair_data_info'] = system_info
    lammps_variables['pair_restart_info'] = restart_info
    lammps_variables['strainrange'] = strainrange

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

    # Run lammps
    output = lmp.run(lammps_command, script_name=lammps_script,
                     mpi_command=mpi_command,
                     logfile=f'cij-{cycle}-log.lammps')
    thermo = output.flatten('all').thermo

    # Extract LAMMPS thermo data. Each term ranges i=0-12 where i=0 is undeformed
    # The remaining values are for -/+ strain pairs in the six unique directions
    lx = uc.set_in_units(thermo.Lx, lammps_units['length'])
    ly = uc.set_in_units(thermo.Ly, lammps_units['length'])
    lz = uc.set_in_units(thermo.Lz, lammps_units['length'])
    xy = uc.set_in_units(thermo.Xy, lammps_units['length'])
    xz = uc.set_in_units(thermo.Xz, lammps_units['length'])
    yz = uc.set_in_units(thermo.Yz, lammps_units['length'])

    pxx = uc.set_in_units(thermo.Pxx, lammps_units['pressure'])
    pyy = uc.set_in_units(thermo.Pyy, lammps_units['pressure'])
    pzz = uc.set_in_units(thermo.Pzz, lammps_units['pressure'])
    pxy = uc.set_in_units(thermo.Pxy, lammps_units['pressure'])
    pxz = uc.set_in_units(thermo.Pxz, lammps_units['pressure'])
    pyz = uc.set_in_units(thermo.Pyz, lammps_units['pressure'])

    pe = uc.set_in_units(thermo.PotEng / system.natoms, lammps_units['energy'])

    # Extract the pressure tensor
    pij = np.array([[pxx[0], pxy[0], pxz[0]],
                    [pxy[0], pyy[0], pyz[0]],
                    [pxz[0], pyz[0], pzz[0]]])

    # Set the six non-zero strain values
    strains = np.array([ (lx[2] -  lx[1])  / lx[0],
                         (ly[4] -  ly[3])  / ly[0],
                         (lz[6] -  lz[5])  / lz[0],
                         (yz[8] -  yz[7])  / lz[0],
                         (xz[10] - xz[9])  / lz[0],
                         (xy[12] - xy[11]) / ly[0] ])

    # Calculate cij using stress changes associated with each non-zero strain
    cij = np.empty((6,6))
    for i in range(6):
        delta_stress = np.array([ pxx[2*i+1]-pxx[2*i+2],
                                  pyy[2*i+1]-pyy[2*i+2],
                                  pzz[2*i+1]-pzz[2*i+2],
                                  pyz[2*i+1]-pyz[2*i+2],
                                  pxz[2*i+1]-pxz[2*i+2],
                                  pxy[2*i+1]-pxy[2*i+2] ])

        cij[i] = delta_stress / strains[i]

    for i in range(6):
        for j in range(i):
            cij[i,j] = cij[j,i] = (cij[i,j] + cij[j,i]) / 2

    C = am.ElasticConstants(Cij=cij)

    results = {}
    results['E_pot'] = pe[0]
    results['pij'] = pij
    results['C'] = C

    return results

2.3. update_box()

[6]:
def update_box(system: am.System,
               C: am.ElasticConstants,
               pij: np.ndarray,
               target_pxx: float = 0.0,
               target_pyy: float = 0.0,
               target_pzz: float = 0.0,
               target_pxy: float = 0.0,
               target_pxz: float = 0.0,
               target_pyz: float = 0.0,
               tol: float = 1e-10) -> am.System:
    """
    Generates a new system with an updated box that attempts to reach the target
    pressure. The new box dimensions are estimated by assuming linear elasticity
    and using the pressure and elastic constants of the current system.

    Parameters
    ----------
    system : atomman.System
        The system to update
    C : atomman.ElasticConstants
        The computed elastic constants for the system.
    pij : numpy.NDArray
        The 3x3 array of pressure tensors computed for the system.
    target_pxx : float, optional
        The value to relax the x tensile pressure component to. Default is
        0.0.
    target_pyy : float, optional
        The value to relax the y tensile pressure component to. Default is
        0.0.
    target_pzz : float, optional
        The value to relax the z tensile pressure component to. Default is
        0.0).
    target_pyz : float, optional
        The value to relax the yz shear pressure component to. Default is
        0.0).
    target_pxz : float, optional
        The value to relax the xz shear pressure component to. Default is
        0.0).
    target_pyz : float, optional
        The value to relax the xy shear pressure component to. Default is
        0.0).
    tol : float, optional
        The target tolerance.  Any strains less than this will be ignored.
        Default value is 1e-10.

    Returns
    -------
    atomman.System
        The System with updated box dimensions.
    """

    # Build the target pij array
    target_pij = np.array([[target_pxx, target_pxy, target_pxz],
                           [target_pxy, target_pyy, target_pyz],
                           [target_pxz, target_pyz, target_pzz]])

    # Adjust pij by the target
    pij = pij - target_pij

    # The stress state is the negative of the pressure state
    σij = -1 * pij

    # Compute the strain associated with the system relative to the target
    ϵij = np.einsum('ijkl,kl->ij', C.Sijkl, σij)
    ϵij[np.abs(ϵij) <= tol] = 0.0

    # Compute new box dimensions by removing the strain
    lx = system.box.lx - ϵij[0,0] * system.box.lx
    ly = system.box.ly - ϵij[1,1] * system.box.ly
    lz = system.box.lz - ϵij[2,2] * system.box.lz
    yz = system.box.yz - 2*ϵij[1,2] * system.box.lz
    xz = system.box.xz - 2*ϵij[0,2] * system.box.lz
    xy = system.box.xy - 2*ϵij[0,1] * system.box.ly

    if lx <= 0.0 or ly <= 0.0 or lz <= 0.0:
        raise RuntimeError('Divergence of box dimensions to <= 0')

    # Duplicate system and change dimensions
    system_new = deepcopy(system)
    system_new.box_set(lx=lx, ly=ly, lz=lz, yz=yz, xz=xz, xy=xy, scale=True)
    system_new.wrap()
    return system_new

2.4. cij_run0.template file

[7]:
with open('cij_run0.template', 'w') as f:
    f.write("""# Evaluates P and Cij using small strains to guess how to relax the box.
# NOTE: Cij predictions may be poor for some crystal structures: use
# elastic_constants_static calculation instead.

box tilt large

<pair_data_info>

change_box all triclinic

# Specify strain
variable strain equal <strainrange>

# 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
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Evaluate and save the initial configuration
run 0
write_restart initial.restart

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply -xx strain
variable delta equal -${strain}*${lx0}
change_box all x delta 0 ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply +xx strain
variable delta equal ${strain}*${lx0}
change_box all x delta 0 ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply -yy strain
variable delta equal -${strain}*${ly0}
change_box all y delta 0 ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply +yy strain
variable delta equal ${strain}*${ly0}
change_box all y delta 0 ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply -zz strain
variable delta equal -${strain}*${lz0}
change_box all z delta 0 ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply +zz strain
variable delta equal ${strain}*${lz0}
change_box all z delta 0 ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply -yz strain
variable delta equal -${strain}*${lz0}
change_box all yz delta ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply +yz strain
variable delta equal ${strain}*${lz0}
change_box all yz delta ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply -xz strain
variable delta equal -${strain}*${lz0}
change_box all xz delta ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply +xz strain
variable delta equal ${strain}*${lz0}
change_box all xz delta ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply -xy strain
variable delta equal -${strain}*${ly0}
change_box all xy delta ${delta} remap units box
run 0

# Reset simulation
clear
box tilt large
<pair_restart_info>
thermo_style custom step lx ly lz yz xz xy pxx pyy pzz pyz pxz pxy v_peatom pe
thermo_modify format float %.13e


# Apply +xy strain
variable delta equal ${strain}*${ly0}
change_box all xy delta ${delta} remap units box
run 0""")

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 = None

# 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 = [3, 3, 3]

# Generate system by supersizing ucell
system = ucell.supersize(*sizemults)
print('# of atoms in system =', system.natoms)
# of atoms in system = 108

3.5. Calculation-specific parameters

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

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

  • convergence_tol is the relative tolerance to use in identifying if the lattice constants have converged.

  • divergence_scale is a factor for identifying if the lattice constants have diverged from the original guess. Divergence is identified if \(a > a^0 d\) or \(a < a^0 / d\), where d is divergence_scale.

[12]:
strainrange = 1e-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')
convergence_tol = 1e-11
divergence_scale = 3.

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_box(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,
                         strainrange = strainrange,
                         tol = convergence_tol,
                         diverge_scale = divergence_scale)
print(results_dict.keys())
dict_keys(['dumpfile_initial', 'symbols_initial', 'dumpfile_final', 'symbols_final', 'lx', 'ly', 'lz', 'xy', 'xz', 'yz', 'E_pot', 'measured_pxx', 'measured_pyy', 'measured_pzz', 'measured_pxy', 'measured_pxz', 'measured_pyz'])

4.2. Report results

Values returned in the results_dict:

  • ‘dumpfile_initial’ (str) - The name of the initial dump file created.

  • ‘symbols_initial’ (list) - The symbols associated with the initial dump file.

  • ‘dumpfile_final’ (str) - The name of the final dump file created.

  • ‘symbols_final’ (list) - The symbols associated with the final dump file.

  • ‘lx’ (float) - The relaxed lx box length.

  • ‘ly’ (float) - The relaxed ly box length.

  • ‘lz’ (float) - The relaxed lz box length.

  • ‘xy’ (float) - The relaxed xy box tilt.

  • ‘xz’ (float) - The relaxed xz box tilt.

  • ‘yz’ (float) - The relaxed yz box tilt.

  • ‘E_pot’ (float) - The potential energy per atom for the final configuration.

  • ‘measured_pxx’ (float) - The measured x tensile pressure 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.

[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'])
initial.dump
('Ni',)
final.dump
('Ni',)
[15]:
length_unit = 'angstrom'
energy_unit = 'eV'

# Show the per atom potential energy
print('E_pot =', uc.get_in_units(results_dict['E_pot'], energy_unit), energy_unit)

# Construct a Box from the returned system dimensions
box = am.Box(lx=results_dict['lx'], ly=results_dict['ly'], lz=results_dict['lz'],
             xy=results_dict['xy'], xz=results_dict['xz'], yz=results_dict['yz'])

# Retrieve lattice constants by dividing by sizemults
print('a =', uc.get_in_units(box.a / sizemults[0], length_unit), length_unit)
print('b =', uc.get_in_units(box.b / sizemults[1], length_unit), length_unit)
print('c =', uc.get_in_units(box.c / sizemults[2], length_unit), length_unit)
print('alpha =', box.alpha)
print('beta = ', box.beta)
print('gamma =', box.gamma)
E_pot = -4.449999998348981 eV
a = 3.519999437539791 angstrom
b = 3.519999437539791 angstrom
c = 3.5199994375397874 angstrom
alpha = 90.0
beta =  90.0
gamma = 90.0
[16]:
pressure_unit = 'GPa'

# Show the computed pressure tensor
print('Pxx =', uc.get_in_units(results_dict['measured_pxx'], pressure_unit), pressure_unit)
print('Pyy =', uc.get_in_units(results_dict['measured_pyy'], pressure_unit), pressure_unit)
print('Pzz =', uc.get_in_units(results_dict['measured_pzz'], pressure_unit), pressure_unit)
print('Pyz =', uc.get_in_units(results_dict['measured_pyz'], pressure_unit), pressure_unit)
print('Pxz =', uc.get_in_units(results_dict['measured_pxz'], pressure_unit), pressure_unit)
print('Pxy =', uc.get_in_units(results_dict['measured_pxy'], pressure_unit), pressure_unit)
Pxx = 9.7049382646425e-11 GPa
Pyy = 9.704980181819598e-11 GPa
Pzz = 9.704606703547402e-11 GPa
Pyz = -4.3238889965929e-16 GPa
Pxz = -7.1938939200957e-16 GPa
Pxy = -7.0239594180462e-16 GPa