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
|
# Numerical difference calculation
# of error in forces, virial stress, and Born matrix
# adjustable parameters
variable nsteps index 500 # length of run
variable nthermo index 10 # thermo output interval
variable ndump index 500 # dump output interval
variable nlat index 3 # size of box
variable fdelta index 1.0e-4 # displacement size
variable vdelta index 1.0e-6 # strain size for numdiff/virial
variable bdelta index 1.0e-8 # strain size for numdiff Born matrix
variable temp index 10.0 # temperature
variable nugget equal 1.0e-6 # regularization for relerr
units metal
atom_style atomic
atom_modify map yes
lattice fcc 5.358000
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
create_box 1 box
create_atoms 1 box
mass 1 39.903
velocity all create ${temp} 2357 mom yes dist gaussian
pair_style lj/cut 5.0
pair_coeff 1 1 0.0102701 3.42
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
timestep 0.001
fix nve all nve
# define numerical force calculation
fix numforce all numdiff ${nthermo} ${fdelta}
variable ferrx atom f_numforce[1]-fx
variable ferry atom f_numforce[2]-fy
variable ferrz atom f_numforce[3]-fz
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
compute faverrsq all reduce ave v_ferrsq
variable fsq atom fx^2+fy^2+fz^2
compute favsq all reduce ave v_fsq
variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget}))
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
# define numerical virial stress tensor calculation
compute myvirial all pressure NULL virial
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
variable errxx equal f_numvirial[1]-c_myvirial[1]
variable erryy equal f_numvirial[2]-c_myvirial[2]
variable errzz equal f_numvirial[3]-c_myvirial[3]
variable erryz equal f_numvirial[4]-c_myvirial[6]
variable errxz equal f_numvirial[5]-c_myvirial[5]
variable errxy equal f_numvirial[6]-c_myvirial[4]
variable verrsq equal "v_errxx^2 + &
v_erryy^2 + &
v_errzz^2 + &
v_erryz^2 + &
v_errxz^2 + &
v_errxy^2"
variable vsq equal "c_myvirial[1]^2 + &
c_myvirial[3]^2 + &
c_myvirial[3]^2 + &
c_myvirial[4]^2 + &
c_myvirial[5]^2 + &
c_myvirial[6]^2"
variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget}))
# define numerical Born matrix calculation
compute bornnum all born/matrix numdiff ${bdelta} myvirial
compute born all born/matrix
variable berr vector c_bornnum-c_born
variable berrsq equal "v_berr[1]^2 + &
v_berr[2]^2 + &
v_berr[3]^2 + &
v_berr[4]^2 + &
v_berr[5]^2 + &
v_berr[6]^2 + &
v_berr[7]^2 + &
v_berr[8]^2 + &
v_berr[9]^2 + &
v_berr[10]^2 + &
v_berr[11]^2 + &
v_berr[12]^2 + &
v_berr[13]^2 + &
v_berr[14]^2 + &
v_berr[15]^2 + &
v_berr[16]^2 + &
v_berr[17]^2 + &
v_berr[18]^2 + &
v_berr[19]^2 + &
v_berr[20]^2 + &
v_berr[21]^2"
variable bsq equal "c_born[1]^2 + &
c_born[2]^2 + &
c_born[3]^2 + &
c_born[4]^2 + &
c_born[5]^2 + &
c_born[6]^2 + &
c_born[7]^2 + &
c_born[8]^2 + &
c_born[9]^2 + &
c_born[10]^2 + &
c_born[11]^2 + &
c_born[12]^2 + &
c_born[13]^2 + &
c_born[14]^2 + &
c_born[15]^2 + &
c_born[16]^2 + &
c_born[17]^2 + &
c_born[18]^2 + &
c_born[19]^2 + &
c_born[20]^2 + &
c_born[21]^2"
variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget}))
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr
thermo ${nthermo}
run ${nsteps}
|
