-
Notifications
You must be signed in to change notification settings - Fork 0
/
in.load
207 lines (163 loc) · 7.02 KB
/
in.load
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
##################################################################################################################
# setup MD calculation to calibrate dislocation velocity as a function of applied stress at finite temperature #
# SG - August 14th, 2013 #
# SG - tested on August 15th, 2013 #
# reference: Groh et al, IJP 2009 #
# TD - Revised for Nickel. March 12, 2015 #
# BH - Revised for a general material. January 13, 2017 #
# BH - Revised for different post-processing method. February 15, 2017 # #
# #
##################################################################################################################
##################################################################################################################
# Edge dislocation
# 1. static to obtain the core structure
# 2. equilibrate temperature at 300K (nve + rescaling)
# 3. apply a constant force to calculate the dislocation velocity. Force is applied on a rigid body
##################################################################################################################
##################################################################################################################
# User set variables
# initTemp : desired temperature
# sigma : shear stress in bar units
# material : material name
# atom_file : name of the file containing the atom positions
#
# Miscellaneous control variables
# equilTime : Number of increment to equilibrate the temperature
# runTime : Number of increment to calibrate the temperature
# energyConv : Energy conversion term from eV to J to convert atomic force to stress
#
# output:
# dump.minimize: atoms before and after minimization
# dump.equilibration: atoms during temperature equilibration
# dump.shear: wrapped atoms in the simulation box (perfect to look at the velocity)
# dump.shear.unwrap: unwrapped atoms (to check that everything is going well)
# output file (log.lammps): monitor the temperature during the MD run
#
##################################################################################################################
# Variable definitions
variable initTemp equal 0. # desired temperature
variable sigma equal 15000. # applied stress in bar
variable material string Ta # material symbol
variable atom_file string data.110_cg # the configuration was generated by SG with the preprocessor dislocation.f90
variable equilTime equal 50000 # number of increment to equilibrate the temperature
variable runTime equal 100000 # number of increment to calibrate the velocity
variable energyConv equal 1602191.7 # conversion factor
dimension 3
boundary p s p
units metal
atom_style atomic
read_data ${atom_file}
# store initial position of bottom and top planes along y
variable tmp0 equal "ylo+15"
variable ylo0 equal ${tmp0}
variable tmp1 equal "yhi-15"
variable yhi0 equal ${tmp1}
# variable for dumping
variable ymid1 equal "0.5*ylo + 0.5*yhi - 30.0"
variable ymiddlenegative equal ${ymid1}
variable ymid2 equal "0.5*(ylo + yhi) + 30.0"
variable ymiddlepositive equal ${ymid2}
# define potential
pair_style eam/fs
pair_coeff * * WTa.eam.fs Ta Ta Ta
neighbor 2.0 bin
neigh_modify delay 5
# definition of the upper and lower blocks
region upper block INF INF ${yhi0} INF INF INF units box
region lower block INF INF INF ${ylo0} INF INF units box
# definition of the group
group upper region upper
group lower region lower
group mobile subtract all upper lower
# fix top and botton group
fix 1 lower setforce 0.0 0.0 NULL
fix 2 upper setforce 0.0 0.0 NULL
# compute specific quantities
compute pot_energy all pe/atom
compute stress all stress/atom NULL
dump 1 all custom 1000 dump.minimize id x y z c_pot_energy
# create the dislocation
timestep 0.010
minimize 0.0 1.0e-8 10000 100000
###################################################
# step 2: equilibrate the temperature #
###################################################
timestep 0.001
undump 1
reset_timestep 0
unfix 1
unfix 2
# define temperature
compute temp1 mobile temp
# initialize the velocities
velocity mobile create ${initTemp} 16723 units box
velocity mobile zero linear
velocity mobile zero angular
# equilibrate the temperature
fix 1 mobile nve
fix 2 mobile temp/rescale 1 ${initTemp} ${initTemp} 1.0 0.5
fix_modify 2 temp temp1
# boundary condition
fix 3 lower setforce 0.0 0.0 NULL
fix 4 upper setforce 0.0 0.0 NULL
thermo 100
thermo_modify temp temp1
thermo_style custom step pe ke temp pxx pyy pzz pxy pyz pxz
dump 1 all custom 50000 dump.equilibration id x y z c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
run ${equilTime}
########################################################################################################
# step 3: MD at finite temperatute and with constant force #
# force (eV/angstrom) = stress (bar) * S (angstrom*angstrom) / (number of atoms in upper) / 1602191.7) #
########################################################################################################
unfix 1
unfix 2
unfix 3
unfix 4
undump 1
uncompute temp1
reset_timestep 0
timestep 0.002
# define the force to apply
variable nupper equal count(upper)
variable nuper1 equal ${nupper}
print "number of atoms in upper == ${nupper}"
variable tmp2 equal "lx"
variable tmp3 equal "lz"
variable tmp4 equal v_tmp2*v_tmp3/v_nupper*v_sigma/v_energyConv
variable appforce equal ${tmp4}
# define temperature
compute temp1 mobile temp
# define velocity on boundary
velocity upper set 0. 0. 0. units box
velocity lower set 0. 0. 0. units box
# define boundary conditions
fix 1 upper setforce NULL 0. NULL
fix 2 lower setforce 0. 0. NULL
fix 3 upper aveforce ${appforce} 0. 0.
fix 4 upper rigid group 1 upper
# no temperature control
fix 5 mobile nve
# Calculate the shear stress xy direction
compute 1 mobile stress/atom NULL
compute 2 mobile reduce sum c_1[4] # XY component
variable shearstress equal "(c_2)/((1*lx*lz*(ly-30)) * 10)" # 10 is to convert bar to MPa; d=1=dimension so shearstress = sum of stress/dV
variable steps equal "dt * step"
# Calculate the shear strain xy direction
compute 3 upper displace/atom
compute 4 upper msd
variable tmp1 equal "c_4[1]" # msd dx
variable L2 equal v_tmp1
thermo 100
thermo_style custom step temp c_4[1] v_L2 v_tmp1
thermo_modify temp temp1
run 0
print "Initial Length x, L2: ${L2}"
variable strainxy2 equal "-((lx) - v_L2)/ly"
variable tmp equal "lx"
variable L1 equal ${tmp}
print "Initial Length x, L1: ${L1}"
variable strainxy equal "-((lx) - v_L1)/ly"
dump 1 all custom 5000 dump.shear id mass x y z
#dump 2 all custom 5000 dump.shear.unwrap id mass xu yu zu
fix press all print 1000 "${steps} ${strainxy} ${strainxy2} ${shearstress}" file values.txt screen no # ps ; no unit ; MPa
run ${runTime}