-
Notifications
You must be signed in to change notification settings - Fork 1.7k
/
in.elastic
210 lines (166 loc) · 6.53 KB
/
in.elastic
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
# Compute elastic constant tensor for a crystal at finite temperature
#
# WARNING: This script attempts to measure small differences in
# statistical averages generated by finite temperature MD simulations.
# This creates many more opportunities for errors to occur, compared
# to the zero temperature version provided in examples/ELASTIC.
# Hence, it is always a good idea to start with the zero temperature
# calculation, before attempting to use this script.
#
# Written by Aidan Thompson (Sandia, athomps@sandia.gov)
#
# This script uses the following three include files.
#
# init.mod (must be modified for different crystal structures)
# Define units, MD parameters, deformation parameters,
# and initial configuration of the atoms and simulation cell.
#
#
# potential.mod (must be modified for different pair styles)
# Define pair style and other attributes
# not stored in restart file
#
#
# displace.mod (displace.mod should not need to be modified)
# Perform positive and negative box displacements
# in direction ${dir} and size ${up}.
# It uses the resultant changes
# in stress to compute one
# row of the elastic stiffness tensor
#
# Inputs variables:
# dir = the Voigt deformation component
# (1,2,3,4,5,6)
# Global constants:
# up = the deformation magnitude (strain units)
# cfac = conversion from LAMMPS pressure units to
# output units for elastic constants
#
#
# To run this on a different system, it should only be necessary to
# modify the files init.mod and potential.mod. In order to calculate
# the elastic constants correctly, care must be taken to specify
# the correct units in init.mod (units, cfac and cunits). It is also
# important to verify that the MD sampling of stress components
# is generating accurate statistical averages.
# One indication of this is that the elastic constants are insensitive
# to the choice of the variable ${up} in init.mod. Another is to
# check for finite size effects.
#
include init.mod
# Compute initial state
variable thermostat equal 1
include potential.mod
run ${nequil}
if "${adiabatic} == 1" &
then "variable thermostat equal 0" &
else "variable thermostat equal 1"
print ${thermostat}
include potential.mod
run ${nrun}
variable pxx0 equal f_avp[1]
variable pyy0 equal f_avp[2]
variable pzz0 equal f_avp[3]
variable pxy0 equal f_avp[4]
variable pxz0 equal f_avp[5]
variable pyz0 equal f_avp[6]
variable tmp equal lx
variable lx0 equal ${tmp}
variable tmp equal ly
variable ly0 equal ${tmp}
variable tmp equal lz
variable lz0 equal ${tmp}
# These formulas define the derivatives w.r.t. strain components
# Constants uses $, variables use v_
variable d1 equal -(v_pxx1-${pxx0})/(v_delta/v_len0)*${cfac}
variable d2 equal -(v_pyy1-${pyy0})/(v_delta/v_len0)*${cfac}
variable d3 equal -(v_pzz1-${pzz0})/(v_delta/v_len0)*${cfac}
variable d4 equal -(v_pyz1-${pyz0})/(v_delta/v_len0)*${cfac}
variable d5 equal -(v_pxz1-${pxz0})/(v_delta/v_len0)*${cfac}
variable d6 equal -(v_pxy1-${pxy0})/(v_delta/v_len0)*${cfac}
# Write restart
write_restart restart.equil
# uxx Perturbation
variable dir equal 1
include displace.mod
# uyy Perturbation
variable dir equal 2
include displace.mod
# uzz Perturbation
variable dir equal 3
include displace.mod
# uyz Perturbation
variable dir equal 4
include displace.mod
# uxz Perturbation
variable dir equal 5
include displace.mod
# uxy Perturbation
variable dir equal 6
include displace.mod
# Output final values
variable C11all equal ${C11}
variable C22all equal ${C22}
variable C33all equal ${C33}
variable C12all equal 0.5*(${C12}+${C21})
variable C13all equal 0.5*(${C13}+${C31})
variable C23all equal 0.5*(${C23}+${C32})
variable C44all equal ${C44}
variable C55all equal ${C55}
variable C66all equal ${C66}
variable C14all equal 0.5*(${C14}+${C41})
variable C15all equal 0.5*(${C15}+${C51})
variable C16all equal 0.5*(${C16}+${C61})
variable C24all equal 0.5*(${C24}+${C42})
variable C25all equal 0.5*(${C25}+${C52})
variable C26all equal 0.5*(${C26}+${C62})
variable C34all equal 0.5*(${C34}+${C43})
variable C35all equal 0.5*(${C35}+${C53})
variable C36all equal 0.5*(${C36}+${C63})
variable C45all equal 0.5*(${C45}+${C54})
variable C46all equal 0.5*(${C46}+${C64})
variable C56all equal 0.5*(${C56}+${C65})
# Average moduli for cubic crystals
variable C11cubic equal (${C11all}+${C22all}+${C33all})/3.0
variable C12cubic equal (${C12all}+${C13all}+${C23all})/3.0
variable C44cubic equal (${C44all}+${C55all}+${C66all})/3.0
variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0
variable shearmodulus1 equal ${C44cubic}
variable shearmodulus2 equal (${C11cubic}-${C12cubic})/2.0
variable poissonratio equal 1.0/(1.0+${C11cubic}/${C12cubic})
# For Stillinger-Weber silicon, the analytical results
# are known to be (E. R. Cowley, 1988):
# C11 = 151.4 GPa
# C12 = 76.4 GPa
# C44 = 56.4 GPa
print "========================================="
print "Components of the Elastic Constant Tensor"
print "========================================="
print "Elastic Constant C11all = ${C11all} ${cunits}"
print "Elastic Constant C22all = ${C22all} ${cunits}"
print "Elastic Constant C33all = ${C33all} ${cunits}"
print "Elastic Constant C12all = ${C12all} ${cunits}"
print "Elastic Constant C13all = ${C13all} ${cunits}"
print "Elastic Constant C23all = ${C23all} ${cunits}"
print "Elastic Constant C44all = ${C44all} ${cunits}"
print "Elastic Constant C55all = ${C55all} ${cunits}"
print "Elastic Constant C66all = ${C66all} ${cunits}"
print "Elastic Constant C14all = ${C14all} ${cunits}"
print "Elastic Constant C15all = ${C15all} ${cunits}"
print "Elastic Constant C16all = ${C16all} ${cunits}"
print "Elastic Constant C24all = ${C24all} ${cunits}"
print "Elastic Constant C25all = ${C25all} ${cunits}"
print "Elastic Constant C26all = ${C26all} ${cunits}"
print "Elastic Constant C34all = ${C34all} ${cunits}"
print "Elastic Constant C35all = ${C35all} ${cunits}"
print "Elastic Constant C36all = ${C36all} ${cunits}"
print "Elastic Constant C45all = ${C45all} ${cunits}"
print "Elastic Constant C46all = ${C46all} ${cunits}"
print "Elastic Constant C56all = ${C56all} ${cunits}"
print "========================================="
print "Average properties for a cubic crystal"
print "========================================="
print "Bulk Modulus = ${bulkmodulus} ${cunits}"
print "Shear Modulus 1 = ${shearmodulus1} ${cunits}"
print "Shear Modulus 2 = ${shearmodulus2} ${cunits}"
print "Poisson Ratio = ${poissonratio}"