Rearranged the directory tree for ELASTIC_T and BORN_MATRIX
This commit is contained in:
Aidan Thompson
3
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/born.out
Normal file
3
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/born.out
Normal file
@ -0,0 +1,3 @@
|
||||
# Time-averaged data for fix CB
|
||||
# TimeStep c_born[1] c_born[2] c_born[3] c_born[4] c_born[5] c_born[6] c_born[7] c_born[8] c_born[9] c_born[10] c_born[11] c_born[12] c_born[13] c_born[14] c_born[15] c_born[16] c_born[17] c_born[18] c_born[19] c_born[20] c_born[21]
|
||||
100000 9548.45 9583.86 9603.94 5344.22 5340.32 5332.84 5332.84 5340.32 6.76151 5.615 -4.57733 5344.22 2.92646 2.31369 -1.40512 -0.756072 2.34424 -1.68733 -1.68733 2.31369 6.76151
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
3.36316 1.87373 1.87607 -0.00346 -0.00172 -0.00104
|
||||
1.87373 3.36170 1.87425 0.00443 0.00033 0.00014
|
||||
1.87607 1.87425 3.36573 0.00143 0.00155 0.00127
|
||||
-0.00346 0.00443 0.00143 1.87425 0.00127 0.00033
|
||||
-0.00172 0.00033 0.00155 0.00127 1.87607 -0.00346
|
||||
-0.00104 0.00014 0.00127 0.00033 -0.00346 1.87373
|
||||
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py
Normal file
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py
Normal file
@ -0,0 +1,118 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
|
||||
import sys
|
||||
import numpy as np
|
||||
|
||||
def reduce_Born(Cf):
|
||||
C = np.zeros((6,6), dtype=np.float64)
|
||||
C[0,0] = Cf[0]
|
||||
C[1,1] = Cf[1]
|
||||
C[2,2] = Cf[2]
|
||||
C[3,3] = Cf[3]
|
||||
C[4,4] = Cf[4]
|
||||
C[5,5] = Cf[5]
|
||||
C[0,1] = Cf[6]
|
||||
C[0,2] = Cf[7]
|
||||
C[0,3] = Cf[8]
|
||||
C[0,4] = Cf[9]
|
||||
C[0,5] = Cf[10]
|
||||
C[1,2] = Cf[11]
|
||||
C[1,3] = Cf[12]
|
||||
C[1,4] = Cf[13]
|
||||
C[1,5] = Cf[14]
|
||||
C[2,3] = Cf[15]
|
||||
C[2,4] = Cf[16]
|
||||
C[2,5] = Cf[17]
|
||||
C[3,4] = Cf[18]
|
||||
C[3,5] = Cf[19]
|
||||
C[4,5] = Cf[20]
|
||||
C = np.where(C,C,C.T)
|
||||
return C
|
||||
|
||||
def compute_delta():
|
||||
D = np.zeros((3,3,3,3))
|
||||
for a in range(3):
|
||||
for b in range(3):
|
||||
for m in range(3):
|
||||
for n in range(3):
|
||||
D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
|
||||
d = np.zeros((6,6))
|
||||
d[0,0] = D[0,0,0,0]
|
||||
d[1,1] = D[1,1,1,1]
|
||||
d[2,2] = D[2,2,2,2]
|
||||
d[3,3] = D[1,2,1,2]
|
||||
d[4,4] = D[0,2,0,2]
|
||||
d[5,5] = D[0,1,0,1]
|
||||
d[0,1] = D[0,0,1,1]
|
||||
d[0,2] = D[0,0,2,2]
|
||||
d[0,3] = D[0,0,1,2]
|
||||
d[0,4] = D[0,0,0,2]
|
||||
d[0,5] = D[0,0,0,1]
|
||||
d[1,2] = D[1,1,2,2]
|
||||
d[1,3] = D[1,1,1,2]
|
||||
d[1,4] = D[1,1,0,2]
|
||||
d[1,5] = D[1,1,0,1]
|
||||
d[2,3] = D[2,2,1,2]
|
||||
d[2,4] = D[2,2,0,2]
|
||||
d[2,5] = D[2,2,0,1]
|
||||
d[3,4] = D[1,2,0,2]
|
||||
d[3,5] = D[1,2,0,1]
|
||||
d[4,5] = D[0,2,0,1]
|
||||
d = np.where(d,d,d.T)
|
||||
return d
|
||||
|
||||
|
||||
def write_matrix(C, filename):
|
||||
with open(filename, 'w') as f:
|
||||
f.write("# Cij Matrix from post process computation\n")
|
||||
for i in C:
|
||||
f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
|
||||
i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
|
||||
)
|
||||
)
|
||||
return
|
||||
|
||||
def main():
|
||||
|
||||
N = 500
|
||||
vol = 27.047271**3 * 10**-30 # m^3
|
||||
T = 60 # K
|
||||
kb = 1.380649 * 10**-23 # J/K
|
||||
kbT = T*kb # J
|
||||
kcalmol2J = 4183.9954/(6.022*10**23)
|
||||
|
||||
born = np.loadtxt('born.out')
|
||||
stre = np.loadtxt('vir.out')
|
||||
stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa
|
||||
try:
|
||||
mean_born = np.mean(born[:, 1:], axis=0)
|
||||
except IndexError:
|
||||
mean_born = born[1:]
|
||||
|
||||
CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa
|
||||
Cs = vol/kbT*np.cov(stre[:,1:].T)
|
||||
Ct = N*kbT/vol * compute_delta()
|
||||
C = CB - Cs + Ct
|
||||
write_matrix(CB, 'born_matrix.out')
|
||||
write_matrix(Cs, 'stre_matrix.out')
|
||||
write_matrix(Ct, 'temp_matrix.out')
|
||||
write_matrix(C, 'full_matrix.out')
|
||||
C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
print(C*10**-9)
|
||||
print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
|
||||
|
||||
return
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
try:
|
||||
main()
|
||||
except KeyboardInterrupt:
|
||||
raise SystemExit("User interruption.")
|
||||
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
2.18161 1.13726 1.16596 -0.01607 -0.02637 0.00291
|
||||
1.13726 2.20242 1.16714 0.00386 -0.05820 0.02644
|
||||
1.16596 1.16714 2.24704 -0.00354 -0.00368 0.02714
|
||||
-0.01607 0.00386 -0.00354 1.43706 0.00210 0.01003
|
||||
-0.02637 -0.05820 -0.00368 0.00210 1.37530 0.01401
|
||||
0.00291 0.02644 0.02714 0.01003 0.01401 1.42403
|
||||
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov
Normal file
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov
Normal file
@ -0,0 +1,154 @@
|
||||
# Numerical difference calculation
|
||||
# of Born matrix
|
||||
|
||||
# Note that because of cubic symmetry and central forces, we have:
|
||||
# C11, pure axial == positive mean value: 1,2,3
|
||||
# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
|
||||
# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
|
||||
# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
units real
|
||||
variable nsteps index 100000 # length of run
|
||||
variable nthermo index 1000 # thermo output interval
|
||||
variable nlat equal 5 # size of box
|
||||
variable T equal 60. # Temperature in K
|
||||
variable rho equal 5.405 # Lattice spacing in A
|
||||
|
||||
atom_style atomic
|
||||
|
||||
lattice fcc ${rho}
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
|
||||
mass * 39.948
|
||||
|
||||
velocity all create ${T} 87287 loop geom
|
||||
velocity all zero linear
|
||||
|
||||
pair_style lj/cut 12.
|
||||
pair_coeff 1 1 0.238067 3.405
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
variable vol equal vol
|
||||
thermo 100
|
||||
fix aL all ave/time 1 1 1 v_vol ave running
|
||||
fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
|
||||
|
||||
run 20000
|
||||
|
||||
unfix NPT
|
||||
|
||||
variable newL equal "f_aL^(1./3.)"
|
||||
change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
|
||||
|
||||
unfix aL
|
||||
|
||||
reset_timestep 0
|
||||
|
||||
# Conversion variables
|
||||
variable kb equal 1.38065e-23 # J/K
|
||||
variable Myvol equal "vol*10^-30" # Volume in m^3
|
||||
variable kbt equal "v_kb*v_T"
|
||||
variable Nat equal atoms
|
||||
variable Rhokbt equal "v_kbt*v_Nat/v_Myvol"
|
||||
variable at2Pa equal 101325
|
||||
variable kcalmol2J equal "4183.9954/(6.022e23)"
|
||||
variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
|
||||
variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms
|
||||
variable Pa2GPa equal 1e-9
|
||||
|
||||
# Born compute giving <C^b> terms
|
||||
compute born all born/matrix
|
||||
# The six virial stress component to compute <C^fl>
|
||||
compute VIR all pressure NULL virial
|
||||
variable s1 equal "-c_VIR[1]*v_at2Pa"
|
||||
variable s2 equal "-c_VIR[2]*v_at2Pa"
|
||||
variable s3 equal "-c_VIR[3]*v_at2Pa"
|
||||
variable s6 equal "-c_VIR[4]*v_at2Pa"
|
||||
variable s5 equal "-c_VIR[5]*v_at2Pa"
|
||||
variable s4 equal "-c_VIR[6]*v_at2Pa"
|
||||
variable press equal press
|
||||
|
||||
|
||||
# Average of Born term and vector to store stress
|
||||
# for post processing
|
||||
fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
|
||||
fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
|
||||
fix APR all ave/time 1 1 1 v_press ave running
|
||||
fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
|
||||
|
||||
thermo ${nthermo}
|
||||
thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
|
||||
thermo_modify line multi
|
||||
|
||||
fix 1 all nvt temp $T $T 100
|
||||
|
||||
run ${nsteps}
|
||||
|
||||
# Compute vector averages
|
||||
# Note the indice switch.
|
||||
# LAMMPS convention is NOT the Voigt notation.
|
||||
variable aves1 equal "ave(f_VEC[1])"
|
||||
variable aves2 equal "ave(f_VEC[2])"
|
||||
variable aves3 equal "ave(f_VEC[3])"
|
||||
variable aves4 equal "ave(f_VEC[6])"
|
||||
variable aves5 equal "ave(f_VEC[5])"
|
||||
variable aves6 equal "ave(f_VEC[4])"
|
||||
|
||||
# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
|
||||
# is numerically instable. Here we go through the <(s-<s>)^2>
|
||||
# definition.
|
||||
|
||||
# Computing difference relative to average values
|
||||
variable ds1 vector "f_VEC[1]-v_aves1"
|
||||
variable ds2 vector "f_VEC[2]-v_aves2"
|
||||
variable ds3 vector "f_VEC[3]-v_aves3"
|
||||
variable ds4 vector "f_VEC[4]-v_aves4"
|
||||
variable ds5 vector "f_VEC[5]-v_aves5"
|
||||
variable ds6 vector "f_VEC[6]-v_aves6"
|
||||
|
||||
# Squaring and averaging
|
||||
variable dds1 vector "v_ds1*v_ds1"
|
||||
variable dds2 vector "v_ds2*v_ds2"
|
||||
variable dds3 vector "v_ds3*v_ds3"
|
||||
variable vars1 equal "ave(v_dds1)"
|
||||
variable vars2 equal "ave(v_dds2)"
|
||||
variable vars3 equal "ave(v_dds3)"
|
||||
variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
|
||||
variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
|
||||
variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
|
||||
|
||||
variable dds12 vector "v_ds1*v_ds2"
|
||||
variable dds13 vector "v_ds1*v_ds3"
|
||||
variable dds23 vector "v_ds2*v_ds3"
|
||||
variable vars12 equal "ave(v_dds12)"
|
||||
variable vars13 equal "ave(v_dds13)"
|
||||
variable vars23 equal "ave(v_dds23)"
|
||||
variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
|
||||
variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
|
||||
variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
|
||||
|
||||
variable dds4 vector "v_ds4*v_ds4"
|
||||
variable dds5 vector "v_ds5*v_ds5"
|
||||
variable dds6 vector "v_ds6*v_ds6"
|
||||
variable vars4 equal "ave(v_dds4)"
|
||||
variable vars5 equal "ave(v_dds5)"
|
||||
variable vars6 equal "ave(v_dds6)"
|
||||
variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
|
||||
variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
|
||||
variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
|
||||
|
||||
variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
|
||||
variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
|
||||
variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
|
||||
|
||||
print """
|
||||
C11 = ${aC11}
|
||||
C12 = ${aC12}
|
||||
C44 = ${aC44}
|
||||
"""
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
1.22342 0.73647 0.71011 0.01261 0.02465 -0.00395
|
||||
0.73647 1.20115 0.70711 0.00057 0.05854 -0.02630
|
||||
0.71011 0.70711 1.16055 0.00497 0.00524 -0.02587
|
||||
0.01261 0.00057 0.00497 0.45813 -0.00083 -0.00969
|
||||
0.02465 0.05854 0.00524 -0.00083 0.52170 -0.01747
|
||||
-0.00395 -0.02630 -0.02587 -0.00969 -0.01747 0.47064
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
0.04187 0.00000 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.04187 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.04187 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.02093 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.02093 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.00000 0.02093
|
||||
100003
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/vir.out
Normal file
100003
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/vir.out
Normal file
File diff suppressed because it is too large
Load Diff
3
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/born.out
Normal file
3
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/born.out
Normal file
@ -0,0 +1,3 @@
|
||||
# Time-averaged data for fix CB
|
||||
# TimeStep c_born[1] c_born[2] c_born[3] c_born[4] c_born[5] c_born[6] c_born[7] c_born[8] c_born[9] c_born[10] c_born[11] c_born[12] c_born[13] c_born[14] c_born[15] c_born[16] c_born[17] c_born[18] c_born[19] c_born[20] c_born[21]
|
||||
100000 9504.48 9542.2 9561.01 5333.83 5329.43 5322.42 5322.45 5329.46 6.62317 5.91783 -4.83884 5333.63 2.95757 2.42725 -0.949434 -0.946756 1.95764 -1.6134 -1.48155 2.30577 6.76235
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
3.35855 1.86892 1.87139 0.00233 0.00218 -0.00179
|
||||
1.86892 3.37104 1.87285 0.00112 0.00085 -0.00007
|
||||
1.87139 1.87285 3.37707 -0.00058 0.00038 -0.00057
|
||||
0.00233 0.00112 -0.00058 1.88326 -0.00039 0.00065
|
||||
0.00218 0.00085 0.00038 -0.00039 1.88229 0.00242
|
||||
-0.00179 -0.00007 -0.00057 0.00065 0.00242 1.87968
|
||||
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py
Normal file
118
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py
Normal file
@ -0,0 +1,118 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
|
||||
import sys
|
||||
import numpy as np
|
||||
|
||||
def reduce_Born(Cf):
|
||||
C = np.zeros((6,6), dtype=np.float64)
|
||||
C[0,0] = Cf[0]
|
||||
C[1,1] = Cf[1]
|
||||
C[2,2] = Cf[2]
|
||||
C[3,3] = Cf[3]
|
||||
C[4,4] = Cf[4]
|
||||
C[5,5] = Cf[5]
|
||||
C[0,1] = Cf[6]
|
||||
C[0,2] = Cf[7]
|
||||
C[0,3] = Cf[8]
|
||||
C[0,4] = Cf[9]
|
||||
C[0,5] = Cf[10]
|
||||
C[1,2] = Cf[11]
|
||||
C[1,3] = Cf[12]
|
||||
C[1,4] = Cf[13]
|
||||
C[1,5] = Cf[14]
|
||||
C[2,3] = Cf[15]
|
||||
C[2,4] = Cf[16]
|
||||
C[2,5] = Cf[17]
|
||||
C[3,4] = Cf[18]
|
||||
C[3,5] = Cf[19]
|
||||
C[4,5] = Cf[20]
|
||||
C = np.where(C,C,C.T)
|
||||
return C
|
||||
|
||||
def compute_delta():
|
||||
D = np.zeros((3,3,3,3))
|
||||
for a in range(3):
|
||||
for b in range(3):
|
||||
for m in range(3):
|
||||
for n in range(3):
|
||||
D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
|
||||
d = np.zeros((6,6))
|
||||
d[0,0] = D[0,0,0,0]
|
||||
d[1,1] = D[1,1,1,1]
|
||||
d[2,2] = D[2,2,2,2]
|
||||
d[3,3] = D[1,2,1,2]
|
||||
d[4,4] = D[0,2,0,2]
|
||||
d[5,5] = D[0,1,0,1]
|
||||
d[0,1] = D[0,0,1,1]
|
||||
d[0,2] = D[0,0,2,2]
|
||||
d[0,3] = D[0,0,1,2]
|
||||
d[0,4] = D[0,0,0,2]
|
||||
d[0,5] = D[0,0,0,1]
|
||||
d[1,2] = D[1,1,2,2]
|
||||
d[1,3] = D[1,1,1,2]
|
||||
d[1,4] = D[1,1,0,2]
|
||||
d[1,5] = D[1,1,0,1]
|
||||
d[2,3] = D[2,2,1,2]
|
||||
d[2,4] = D[2,2,0,2]
|
||||
d[2,5] = D[2,2,0,1]
|
||||
d[3,4] = D[1,2,0,2]
|
||||
d[3,5] = D[1,2,0,1]
|
||||
d[4,5] = D[0,2,0,1]
|
||||
d = np.where(d,d,d.T)
|
||||
return d
|
||||
|
||||
|
||||
def write_matrix(C, filename):
|
||||
with open(filename, 'w') as f:
|
||||
f.write("# Cij Matrix from post process computation\n")
|
||||
for i in C:
|
||||
f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
|
||||
i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
|
||||
)
|
||||
)
|
||||
return
|
||||
|
||||
def main():
|
||||
|
||||
N = 500
|
||||
vol = 27.047271**3 * 10**-30 # m^3
|
||||
T = 60 # K
|
||||
kb = 1.380649 * 10**-23 # J/K
|
||||
kbT = T*kb # J
|
||||
kcalmol2J = 4183.9954/(6.022*10**23)
|
||||
|
||||
born = np.loadtxt('born.out')
|
||||
stre = np.loadtxt('vir.out')
|
||||
stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa
|
||||
try:
|
||||
mean_born = np.mean(born[:, 1:], axis=0)
|
||||
except IndexError:
|
||||
mean_born = born[1:]
|
||||
|
||||
CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa
|
||||
Cs = vol/kbT*np.cov(stre[:,1:].T)
|
||||
Ct = N*kbT/vol * compute_delta()
|
||||
C = CB - Cs + Ct
|
||||
write_matrix(CB, 'born_matrix.out')
|
||||
write_matrix(Cs, 'stre_matrix.out')
|
||||
write_matrix(Ct, 'temp_matrix.out')
|
||||
write_matrix(C, 'full_matrix.out')
|
||||
C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
|
||||
eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
|
||||
eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
|
||||
print(C*10**-9)
|
||||
print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
|
||||
|
||||
return
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
try:
|
||||
main()
|
||||
except KeyboardInterrupt:
|
||||
raise SystemExit("User interruption.")
|
||||
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
2.29922 1.17439 1.20854 -0.01653 -0.01684 0.01188
|
||||
1.17439 2.20673 1.21718 -0.00781 -0.00753 0.00867
|
||||
1.20854 1.21718 2.30804 0.01535 -0.01596 0.00426
|
||||
-0.01653 -0.00781 0.01535 1.47647 -0.01355 -0.01601
|
||||
-0.01684 -0.00753 -0.01596 -0.01355 1.37905 0.01975
|
||||
0.01188 0.00867 0.00426 -0.01601 0.01975 1.40170
|
||||
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov
Normal file
154
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov
Normal file
@ -0,0 +1,154 @@
|
||||
# Numerical difference calculation
|
||||
# of Born matrix
|
||||
|
||||
# Note that because of cubic symmetry and central forces, we have:
|
||||
# C11, pure axial == positive mean value: 1,2,3
|
||||
# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
|
||||
# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
|
||||
# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
|
||||
|
||||
# adjustable parameters
|
||||
|
||||
units real
|
||||
variable nsteps index 100000 # length of run
|
||||
variable nthermo index 1000 # thermo output interval
|
||||
variable nlat equal 5 # size of box
|
||||
variable T equal 60. # Temperature in K
|
||||
variable rho equal 5.405 # Lattice spacing in A
|
||||
|
||||
atom_style atomic
|
||||
|
||||
lattice fcc ${rho}
|
||||
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
|
||||
mass * 39.948
|
||||
|
||||
velocity all create ${T} 87287 loop geom
|
||||
velocity all zero linear
|
||||
|
||||
pair_style lj/cut 12.
|
||||
pair_coeff 1 1 0.238067 3.405
|
||||
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 delay 0 check no
|
||||
|
||||
variable vol equal vol
|
||||
thermo 100
|
||||
fix aL all ave/time 1 1 1 v_vol ave running
|
||||
fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
|
||||
|
||||
run 20000
|
||||
|
||||
unfix NPT
|
||||
|
||||
variable newL equal "f_aL^(1./3.)"
|
||||
change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
|
||||
|
||||
unfix aL
|
||||
|
||||
reset_timestep 0
|
||||
|
||||
# Conversion variables
|
||||
variable kb equal 1.38065e-23 # J/K
|
||||
variable Myvol equal "vol*10^-30" # Volume in m^3
|
||||
variable kbt equal "v_kb*v_T"
|
||||
variable Nat equal atoms
|
||||
variable Rhokbt equal "v_kbt*v_Nat/v_Myvol"
|
||||
variable at2Pa equal 101325
|
||||
variable kcalmol2J equal "4183.9954/(6.022e23)"
|
||||
variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
|
||||
variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms
|
||||
variable Pa2GPa equal 1e-9
|
||||
|
||||
# Born compute giving <C^b> terms
|
||||
# The six virial stress component to compute <C^fl>
|
||||
compute VIR all pressure NULL virial
|
||||
compute born all born/matrix numdiff 1e-6 VIR
|
||||
variable s1 equal "-c_VIR[1]*v_at2Pa"
|
||||
variable s2 equal "-c_VIR[2]*v_at2Pa"
|
||||
variable s3 equal "-c_VIR[3]*v_at2Pa"
|
||||
variable s6 equal "-c_VIR[4]*v_at2Pa"
|
||||
variable s5 equal "-c_VIR[5]*v_at2Pa"
|
||||
variable s4 equal "-c_VIR[6]*v_at2Pa"
|
||||
variable press equal press
|
||||
|
||||
|
||||
# Average of Born term and vector to store stress
|
||||
# for post processing
|
||||
fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
|
||||
fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
|
||||
fix APR all ave/time 1 1 1 v_press ave running
|
||||
fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
|
||||
|
||||
thermo ${nthermo}
|
||||
thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
|
||||
thermo_modify line multi
|
||||
|
||||
fix 1 all nvt temp $T $T 100
|
||||
|
||||
run ${nsteps}
|
||||
|
||||
# Compute vector averages
|
||||
# Note the indice switch.
|
||||
# LAMMPS convention is NOT the Voigt notation.
|
||||
variable aves1 equal "ave(f_VEC[1])"
|
||||
variable aves2 equal "ave(f_VEC[2])"
|
||||
variable aves3 equal "ave(f_VEC[3])"
|
||||
variable aves4 equal "ave(f_VEC[6])"
|
||||
variable aves5 equal "ave(f_VEC[5])"
|
||||
variable aves6 equal "ave(f_VEC[4])"
|
||||
|
||||
# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
|
||||
# is numerically instable. Here we go through the <(s-<s>)^2>
|
||||
# definition.
|
||||
|
||||
# Computing difference relative to average values
|
||||
variable ds1 vector "f_VEC[1]-v_aves1"
|
||||
variable ds2 vector "f_VEC[2]-v_aves2"
|
||||
variable ds3 vector "f_VEC[3]-v_aves3"
|
||||
variable ds4 vector "f_VEC[4]-v_aves4"
|
||||
variable ds5 vector "f_VEC[5]-v_aves5"
|
||||
variable ds6 vector "f_VEC[6]-v_aves6"
|
||||
|
||||
# Squaring and averaging
|
||||
variable dds1 vector "v_ds1*v_ds1"
|
||||
variable dds2 vector "v_ds2*v_ds2"
|
||||
variable dds3 vector "v_ds3*v_ds3"
|
||||
variable vars1 equal "ave(v_dds1)"
|
||||
variable vars2 equal "ave(v_dds2)"
|
||||
variable vars3 equal "ave(v_dds3)"
|
||||
variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
|
||||
variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
|
||||
variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
|
||||
|
||||
variable dds12 vector "v_ds1*v_ds2"
|
||||
variable dds13 vector "v_ds1*v_ds3"
|
||||
variable dds23 vector "v_ds2*v_ds3"
|
||||
variable vars12 equal "ave(v_dds12)"
|
||||
variable vars13 equal "ave(v_dds13)"
|
||||
variable vars23 equal "ave(v_dds23)"
|
||||
variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
|
||||
variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
|
||||
variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
|
||||
|
||||
variable dds4 vector "v_ds4*v_ds4"
|
||||
variable dds5 vector "v_ds5*v_ds5"
|
||||
variable dds6 vector "v_ds6*v_ds6"
|
||||
variable vars4 equal "ave(v_dds4)"
|
||||
variable vars5 equal "ave(v_dds5)"
|
||||
variable vars6 equal "ave(v_dds6)"
|
||||
variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
|
||||
variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
|
||||
variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
|
||||
|
||||
variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
|
||||
variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
|
||||
variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
|
||||
|
||||
print """
|
||||
C11 = ${aC11}
|
||||
C12 = ${aC12}
|
||||
C44 = ${aC44}
|
||||
"""
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
1.10120 0.69454 0.66285 0.01885 0.01902 -0.01367
|
||||
0.69454 1.20617 0.65567 0.00893 0.00839 -0.00873
|
||||
0.66285 0.65567 1.11090 -0.01593 0.01634 -0.00483
|
||||
0.01885 0.00893 -0.01593 0.42772 0.01316 0.01666
|
||||
0.01902 0.00839 0.01634 0.01316 0.52416 -0.01733
|
||||
-0.01367 -0.00873 -0.00483 0.01666 -0.01733 0.49891
|
||||
@ -0,0 +1,7 @@
|
||||
# Cij Matrix from post process computation
|
||||
0.04187 0.00000 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.04187 0.00000 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.04187 0.00000 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.02093 0.00000 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.02093 0.00000
|
||||
0.00000 0.00000 0.00000 0.00000 0.00000 0.02093
|
||||
100003
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/vir.out
Normal file
100003
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/vir.out
Normal file
File diff suppressed because it is too large
Load Diff
13
examples/ELASTIC_T/BORN_MATRIX/Argon/README.md
Normal file
13
examples/ELASTIC_T/BORN_MATRIX/Argon/README.md
Normal file
@ -0,0 +1,13 @@
|
||||
This repository is a test case for the compute born/matrix. It provides short
|
||||
scripts creating argon fcc crystal and computing the Born term using the two
|
||||
methods described in the documentation.
|
||||
|
||||
In the __Analytical__ directory the terms are computed using the analytical
|
||||
derivation of the Born term for the lj/cut pair style.
|
||||
|
||||
In the __Numdiff__ directory, the Born term is evaluated through small
|
||||
numerical differences of the stress tensor. This method can be used with any
|
||||
interaction potential.
|
||||
|
||||
Both script show examples on how to compute the full Cij elastic stiffness
|
||||
tensor in LAMMPS.
|
||||
@ -2,7 +2,7 @@
|
||||
# See in.elastic for more info.
|
||||
|
||||
# we must undefine any fix ave/* fix before using reset_timestep
|
||||
if "$(is_defined(fix,avp)" then "unfix avp"
|
||||
if "$(is_defined(fix,avp))" then "unfix avp"
|
||||
reset_timestep 0
|
||||
|
||||
# Choose potential
|
||||
Reference in New Issue
Block a user