diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/README b/examples/ELASTIC_T/BORN_MATRIX/Silicon/README new file mode 100644 index 0000000000..1d01f60477 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/README @@ -0,0 +1,37 @@ +This directory shows how to use the `fix born` command +to calculate the full matrix of elastic constants +for cubic diamond at finite temperature +running the Stillinger-Weber potential. + +The input script `in.elastic` can be run +directly from LAMMPS, or via a Python wrapper +script. + +to run directly from LAMMPS, use: + +mpirun -np 4 lmp.exe -in in.elastic + +This simulates an orthorhombic box with the cubic crystal axes +aligned with x, y, and z. +The default settings replicate the 1477~K benchmark of +Kluge, Ray, and Rahman (1986) that is Ref.[15] in: +Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265 + +The script contains many adjustable parameters that can be used +to generate different crystal structures, supercell sizes, +and sampling rates. + +to run via the Python wrapper, use: + +mpirun -np 4 python elastic.py + +This will first run the orthorhombic supercell as before, +follows by an equivalent simulation using a triclinic structure. +The script shows how the standard triclinic primitive cell for cubic diamond +can be rotated in to the LAMMPS upper triangular frame. The resultant +elastic constant matrix does not exhibit the standard symmetries of cubic crystals. +However, the matrix is then rotated back to the standard orientation +to recover the cubic symmetry form of the elastic matrix, +resulting in elastic constants that are the same for both +simulations, modulo statistical uncertainty. + diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic.py b/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic.py new file mode 100755 index 0000000000..1786525317 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python -i +# preceding line should have path for Python on your machine + +# elastic.py +# Purpose: demonstrate elastic constant calculation for +# two different crystal supercells, one with non-standard orientation +# +# Syntax: elastic.py +# uses in.elastic as LAMMPS input script + +from __future__ import print_function +import sys +from elastic_utils import * + +np.set_printoptions(precision = 3, suppress=True) + +# setup MPI, if wanted +me = 0 +# uncomment this if running in parallel via mpi4py +from mpi4py import MPI +me = MPI.COMM_WORLD.Get_rank() +nprocs = MPI.COMM_WORLD.Get_size() + +# cubic diamond lattice constants + +alat = 5.457 + +# define the cubic diamond orthorhombic supercell +# with 8 atoms + +basisstring = "" +origin = np.zeros(3) +bond = np.ones(3)*0.25 +b = origin +basisstring += "basis %g %g %g " % (b[0],b[1],b[2]) +b = bond +basisstring += "basis %g %g %g " % (b[0],b[1],b[2]) + +for i in range(3): + b = 2*bond + b[i] = 0 + basisstring += "basis %g %g %g " % (b[0],b[1],b[2]) + b += bond + basisstring += "basis %g %g %g " % (b[0],b[1],b[2]) + +hmat = np.eye(3) + +varlist = { + "a":alat, + "a1x":hmat[0,0], + "a2x":hmat[0,1], + "a2y":hmat[1,1], + "a3x":hmat[0,2], + "a3y":hmat[1,2], + "a3z":hmat[2,2], + "l":alat, + "basis":basisstring, + "nlat":3, +} + +cmdargs = gen_varargs(varlist) +cij_ortho = calculate_cij(cmdargs) + +# define the cubic diamond triclinic primitive cell +# with 2 atoms + +basisstring = "" +origin = np.zeros(3) +bond = np.ones(3)*0.25 +b = origin +basisstring += "basis %g %g %g " % (b[0],b[1],b[2]) +b = bond +basisstring += "basis %g %g %g " % (b[0],b[1],b[2]) + +hmat1 = np.array([[1, 1, 0], [0, 1, 1], [1, 0, 1]]).T/np.sqrt(2) + +# rotate primitive cell to LAMMPS orientation (uper triangular) + +qmat, rmat = np.linalg.qr(hmat1) +ss = np.diagflat(np.sign(np.diag(rmat))) +rot = ss @ qmat.T +hmat2 = rot @ hmat1 + +varlist = { + "a":alat, + "a1x":hmat2[0,0], + "a2x":hmat2[0,1], + "a2y":hmat2[1,1], + "a3x":hmat2[0,2], + "a3y":hmat2[1,2], + "a3z":hmat2[2,2], + "l":alat/2**0.5, + "basis":basisstring, + "nlat":5, +} + +cmdargs = gen_varargs(varlist) +cij_tri = calculate_cij(cmdargs) + +if me == 0: + print("\nPython output:") + print("C_ortho = \n",cij_ortho) + print() + print("C_tri = \n",cij_tri) + print() + + cij_tri_rot = rotate_cij(cij_tri, rot.T) + + print("C_tri(rotated back) = \n",cij_tri_rot) + print() + + print("C_ortho-C_tri = \n", cij_ortho-cij_tri_rot) + print() diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py b/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py new file mode 100644 index 0000000000..0daf62b1b0 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py @@ -0,0 +1,110 @@ +import numpy as np +from lammps import lammps, LAMMPS_INT, LMP_STYLE_GLOBAL, LMP_VAR_EQUAL, LMP_VAR_ATOM + +# method for rotating elastic constants + +def rotate_cij(cij, r): + + # K_1 + # R_11^2 R_12^2 R_13^2 + # R_21^2 R_22^2 R_23^2 + # R_31^2 R_32^2 R_33^2 + + k1 = r*r + + # K_2 + # R_12.R_13 R_13.R_11 R_11.R_12 + # R_22.R_23 R_23.R_21 R_21.R_22 + # R_32.R_33 R_33.R_31 R_31.R_32 + + k2 = np.array([ + [r[0][1]*r[0][2], r[0][2]*r[0][0], r[0][0]*r[0][1]], + [r[1][1]*r[1][2], r[1][2]*r[1][0], r[1][0]*r[1][1]], + [r[2][1]*r[2][2], r[2][2]*r[2][0], r[2][0]*r[2][1]], + ]) + + # K_3 + # R_21.R_31 R_22.R_32 R_23.R_33 + # R_31.R_11 R_32.R_12 R_33.R_13 + # R_11.R_21 R_12.R_22 R_13.R_23 + + k3 = np.array([ + [r[1][0]*r[2][0], r[1][1]*r[2][1], r[1][2]*r[2][2]], + [r[2][0]*r[0][0], r[2][1]*r[0][1], r[2][2]*r[0][2]], + [r[0][0]*r[1][0], r[0][1]*r[1][1], r[0][2]*r[1][2]], + ]) + + # K_4a + # R_22.R_33 R_23.R_31 R_21.R_32 + # R_32.R_13 R_33.R_11 R_31.R_12 + # R_12.R_23 R_13.R_21 R_11.R_22 + + k4a = np.array([ + [r[1][1]*r[2][2], r[1][2]*r[2][0], r[1][0]*r[2][1]], + [r[2][1]*r[0][2], r[2][2]*r[0][0], r[2][0]*r[0][1]], + [r[0][1]*r[1][2], r[0][2]*r[1][0], r[0][0]*r[1][1]], + ]) + + # K_4b + # R_23.R_32 R_21.R_33 R_22.R_31 + # R_33.R_12 R_31.R_13 R_32.R_11 + # R_13.R_22 R_11.R_23 R_12.R_21 + + k4b = np.array([ + [r[1][2]*r[2][1], r[1][0]*r[2][2], r[1][1]*r[2][0]], + [r[2][2]*r[0][1], r[2][0]*r[0][2], r[2][1]*r[0][0]], + [r[0][2]*r[1][1], r[0][0]*r[1][2], r[0][1]*r[1][0]], + ]) + + k = np.block([[k1, 2*k2],[k3, k4a+k4b]]) + cijrot = k.dot(cij.dot(k.T)) + return cijrot + +def calculate_cij(cmdargs): + lmp = lammps(cmdargs=cmdargs) + lmp.file("in.elastic") + + C11 = lmp.extract_variable("C11",None, LMP_VAR_EQUAL) + C22 = lmp.extract_variable("C22",None, LMP_VAR_EQUAL) + C33 = lmp.extract_variable("C33",None, LMP_VAR_EQUAL) + C44 = lmp.extract_variable("C44",None, LMP_VAR_EQUAL) + C55 = lmp.extract_variable("C55",None, LMP_VAR_EQUAL) + C66 = lmp.extract_variable("C66",None, LMP_VAR_EQUAL) + + C12 = lmp.extract_variable("C12",None, LMP_VAR_EQUAL) + C13 = lmp.extract_variable("C13",None, LMP_VAR_EQUAL) + C14 = lmp.extract_variable("C14",None, LMP_VAR_EQUAL) + C15 = lmp.extract_variable("C15",None, LMP_VAR_EQUAL) + C16 = lmp.extract_variable("C16",None, LMP_VAR_EQUAL) + + C23 = lmp.extract_variable("C23",None, LMP_VAR_EQUAL) + C24 = lmp.extract_variable("C24",None, LMP_VAR_EQUAL) + C25 = lmp.extract_variable("C25",None, LMP_VAR_EQUAL) + C26 = lmp.extract_variable("C26",None, LMP_VAR_EQUAL) + + C34 = lmp.extract_variable("C34",None, LMP_VAR_EQUAL) + C35 = lmp.extract_variable("C35",None, LMP_VAR_EQUAL) + C36 = lmp.extract_variable("C36",None, LMP_VAR_EQUAL) + + C45 = lmp.extract_variable("C45",None, LMP_VAR_EQUAL) + C46 = lmp.extract_variable("C46",None, LMP_VAR_EQUAL) + + C56 = lmp.extract_variable("C56",None, LMP_VAR_EQUAL) + + cij = np.array([ + [C11, C12, C13, C14, C15, C16], + [ 0, C22, C23, C24, C25, C26], + [ 0, 0, C33, C34, C35, C36], + [ 0, 0, 0, C44, C45, C46], + [ 0, 0, 0, 0, C55, C56], + [ 0, 0, 0, 0, 0, C66], + ]) + cij = np.triu(cij) + np.tril(cij.T, -1) + + return cij + +def gen_varargs(varlist): + cmdargs = [] + for key in varlist: + cmdargs += ["-var",key,str(varlist[key])] + return cmdargs diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in index 09aa85a489..292fb0728a 100644 --- a/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in @@ -1,7 +1,6 @@ # NOTE: This script can be modified for different atomic structures, # units, etc. See in.elastic for more info. # - # Define MD parameters # These can be modified by the user # These settings replicate the 1477~K benchmark of @@ -37,19 +36,34 @@ variable nrepeatborn equal floor(${nfreq}/${neveryborn}) # number of samples variable nequil equal 10*${nthermo} # length of equilibration run variable nrun equal 100*${nthermo} # length of equilibrated run -# generate the box and atom positions using a diamond lattice +# this generates a general triclinic cell +# conforming to LAMMPS cell (upper triangular) units metal +box tilt large -boundary p p p +# unit lattice vectors are +# a1 = (a1x 0 0) +# a2 = (a2x a2y 0) +# a3 = (a3x a3y a3z) -# this generates a standard 8-atom cubic cell +variable a1x index 1 +variable a2x index 0 +variable a2y index 1 +variable a3x index 0 +variable a3y index 0 +variable a3z index 1 +variable atmp equal $a +variable l index $a +variable basis index "basis 0 0 0 basis 0.25 0.25 0.25 basis 0 0.5 0.5 basis 0.25 0.75 0.75 basis 0.5 0 0.5 basis 0.75 0.25 0.75 basis 0.5 0.5 0 basis 0.75 0.75 0.25" +lattice custom ${l} & + a1 ${a1x} 0 0 & + a2 ${a2x} ${a2y} 0 & + a3 ${a3x} ${a3y} ${a3z} & + ${basis} & + spacing 1 1 1 -lattice diamond $a -region box prism 0 1 0 1 0 1 0 0 0 - -# this generates a 2-atom triclinic cell -#include tri.in +region box prism 0 ${a1x} 0 ${a2y} 0 ${a3z} ${a2x} ${a3x} ${a3y} create_box 1 box create_atoms 1 box @@ -57,3 +71,4 @@ mass 1 ${mass1} replicate ${nlat} ${nlat} ${nlat} velocity all create ${temp} 87287 + diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in deleted file mode 100644 index 20bfca16ec..0000000000 --- a/examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in +++ /dev/null @@ -1,26 +0,0 @@ -# this generates a 2-atom triclinic cell -# due to rotation on to x-axis, -# elastic constant analysis is not working yet - -# unit lattice vectors are -# a1 = (1 0 0) -# a2 = (1/2 sqrt3/2 0) -# a3 = (1/2 1/(2sqrt3) sqrt2/sqrt3) - -variable a1x equal 1 -variable a2x equal 1/2 -variable a2y equal sqrt(3)/2 -variable a3x equal 1/2 -variable a3y equal 1/(2*sqrt(3)) -variable a3z equal sqrt(2/3) -variable l equal $a/sqrt(2) - -lattice custom ${l} & - a1 ${a1x} 0 0 & - a2 ${a2x} ${a2y} 0.0 & - a3 ${a3x} ${a3y} ${a3z} & - basis 0 0 0 & - basis 0.25 0.25 0.25 & - spacing 1 1 1 - -region box prism 0 ${a1x} 0 ${a2y} 0 ${a3z} ${a2x} ${a3x} ${a3y}