Added Python wrapper that handles rotations to and from LAMMPS frame
This commit is contained in:
Aidan Thompson
37
examples/ELASTIC_T/BORN_MATRIX/Silicon/README
Normal file
37
examples/ELASTIC_T/BORN_MATRIX/Silicon/README
Normal file
@ -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.
|
||||
|
||||
113
examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic.py
Executable file
113
examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic.py
Executable file
@ -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()
|
||||
110
examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py
Normal file
110
examples/ELASTIC_T/BORN_MATRIX/Silicon/elastic_utils.py
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
|
||||
|
||||
@ -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}
|
||||
Reference in New Issue
Block a user