forked from libAtoms/pymatnest
-
Notifications
You must be signed in to change notification settings - Fork 0
/
fortran_MC_MD.F90
317 lines (255 loc) · 9.13 KB
/
fortran_MC_MD.F90
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
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
function fortran_seed_size()
integer :: fortran_seed_size
call random_seed(size=fortran_seed_size)
end function fortran_seed_size
subroutine fortran_set_seed(n_seed, seed)
integer :: seed(n_seed)
call random_seed(put=seed)
end subroutine fortran_set_seed
subroutine fortran_MC_atom_velo(N, vel, mass, n_steps, step_size, nD, KEmax, final_KE, n_try, n_accept)
implicit none
integer :: N
double precision :: vel(3,N), mass(N)
integer :: n_steps, nD
double precision :: step_size, KEmax, final_KE
integer :: n_try, n_accept
integer :: d_i
double precision :: d_r, KE, dKE, d_vel(3)
integer :: i_step, i_at, t_i
integer :: order(N)
n_try = N*n_steps
n_accept = 0
KE = 0.5*sum(spread(mass,1,3)*vel**2)
do i_step=1, n_steps
do i_at=1, N
order(i_at) = i_at
end do
do i_at=1, N-1
call random_number(d_r); d_i = floor(d_r*(N-i_at+1))+i_at
if (d_i /= i_at) then
t_i = order(i_at)
order(i_at) = order(d_i)
order(d_i) = t_i
endif
end do
do i_at=1, N
d_i = order(i_at)
call random_number(d_vel)
d_vel = 2.0*step_size*(d_vel-0.5)
if (nD==2) d_vel(3) = 0.0
dKE = 0.5*mass(d_i)*(sum((vel(:,d_i)+d_vel(:))**2) - sum(vel(:,d_i)**2))
if (KE + dKE < KEmax) then
vel(1:3,d_i) = vel(1:3,d_i) + d_vel(1:3)
KE = KE + dKE
n_accept = n_accept + 1
endif
end do
end do
final_KE = KE
end subroutine fortran_MC_atom_velo
subroutine fortran_MC_atom(N, Z, pos, vel, mass, n_extra_data, extra_data, cell, n_steps, &
step_size_pos, step_size_vel, Emax, nD, KEmax, final_E, n_try, n_accept_pos, n_accept_vel)
implicit none
integer :: N
integer :: Z(N)
double precision :: pos(3,N), vel(3,N), mass(N), cell(3,3)
integer :: n_extra_data
double precision :: extra_data(n_extra_data,N)
integer :: n_steps
double precision :: step_size_pos, step_size_vel, Emax, KEmax, final_E
integer :: n_try, n_accept_pos, n_accept_vel, nD
logical :: do_vel
integer :: d_i
double precision :: d_r, E, dE, KE, dKE, d_pos(3), d_vel(3)
double precision, external :: ll_eval_energy
integer, external :: ll_move_atom_1
integer :: i_step, i_at, t_i
integer :: order(N)
double precision :: vel_pos_rv
do_vel = (step_size_vel /= 0.0)
n_try = N*n_steps
n_accept_pos = 0
n_accept_vel = 0
E = ll_eval_energy(N, Z, pos, n_extra_data, extra_data, cell)
if (do_vel) then
KE = 0.5*sum(spread(mass,1,3)*vel**2)
E = E + KE
endif
do i_step=1, n_steps
do i_at=1, N
order(i_at) = i_at
end do
do i_at=1, N-1
call random_number(d_r); d_i = floor(d_r*(N-i_at+1))+i_at
if (d_i /= i_at) then
t_i = order(i_at)
order(i_at) = order(d_i)
order(d_i) = t_i
endif
end do
do i_at=1, N
d_i = order(i_at)
if (do_vel) then
call random_number(d_vel)
d_vel = 2.0*step_size_vel*(d_vel-0.5)
if (nD==2) d_vel(3)=0.0
call random_number(vel_pos_rv)
endif
if (do_vel .and. vel_pos_rv < 0.5) then
dKE = 0.5*mass(d_i)*(sum((vel(:,d_i)+d_vel(:))**2) - sum(vel(:,d_i)**2))
dE = dKE
if (E + dE < Emax .and. (KEmax <= 0.0 .or. KE + dKE < KEmax)) then
vel(1:3,d_i) = vel(1:3,d_i) + d_vel(1:3)
E = E + dE
KE = KE + dKE
n_accept_vel = n_accept_vel + 1
endif
endif
call random_number(d_pos)
d_pos = 2.0*step_size_pos*(d_pos-0.5)
if (nD==2) d_pos(3)=0.0
n_accept_pos = n_accept_pos + ll_move_atom_1(N, Z, pos, n_extra_data, extra_data, cell, d_i, d_pos, Emax-E, dE)
E = E + dE
if (do_vel .and. vel_pos_rv >= 0.5) then
dKE = 0.5*mass(d_i)*(sum((vel(:,d_i)+d_vel(:))**2) - sum(vel(:,d_i)**2))
dE = dKE
if (E + dE < Emax .and. (KEmax <= 0.0 .or. KE + dKE < KEmax)) then
vel(1:3,d_i) = vel(1:3,d_i) + d_vel(1:3)
E = E + dE
KE = KE + dKE
n_accept_vel = n_accept_vel + 1
endif
endif
end do
end do
final_E = E
end subroutine fortran_MC_atom
subroutine rotate_dir(N, dir, max_ang)
implicit none
integer :: N
double precision :: dir(3, N), max_ang
double precision :: d_r, ang, v1, v2, c, s
integer :: ii, i1, i2, c1, c2
do ii=1, 3*N
call random_number(d_r); i1 = floor(d_r*N)+1
call random_number(d_r); i2 = floor(d_r*N)+2
call random_number(d_r); c1 = floor(d_r*3)+1
call random_number(d_r); c2 = floor(d_r*3)+2
call random_number(d_r); ang = max_ang*2.0*(d_r-0.5)
c = cos(ang)
s = sin(ang)
v1 = c*dir(i1,c1) + s*dir(i2,c2)
v2 = -s*dir(i1,c1) + c*dir(i2,c2)
dir(i1,c1) = v1
dir(i2,c2) = v2
end do
end subroutine rotate_dir
subroutine fortran_GMC_atom(N, Z, pos, mass, n_extra_data, extra_data, cell, n_steps, &
Emax, final_E, n_try, n_accept, d_pos, no_reverse, pert_ang, debug)
implicit none
integer :: N
integer :: Z(N)
double precision :: pos(3,N), mass(N), cell(3,3)
integer :: n_extra_data
double precision :: extra_data(n_extra_data,N)
integer :: n_steps
double precision :: Emax, final_E
integer :: n_try, n_accept
double precision :: d_pos(3,N)
integer :: no_reverse
double precision :: pert_ang
integer :: debug
double precision :: E, Fhat(3,N), d_pos_reflect(3,N), &
last_good_pos(3,N), last_good_d_pos(3,N), last_good_E
double precision, external :: ll_eval_energy, ll_eval_forces
integer :: n_reflect, n_reverse
integer :: i_step
n_reflect = 0
n_reverse = 0
if (no_reverse /= 0) then
last_good_pos = pos
last_good_d_pos = d_pos
last_good_E = ll_eval_energy(N, Z, pos, n_extra_data, extra_data, cell)
end if
do i_step=1, n_steps
if (no_reverse == 0) then
last_good_pos = pos
last_good_d_pos = d_pos
endif
call rotate_dir(N, d_pos, pert_ang)
pos = pos + d_pos
E = ll_eval_forces(N, Z, pos, n_extra_data, extra_data, cell, Fhat)
if (isnan(E) .or. E >= Emax) then ! reflect or reverse
! E = ll_eval_forces(N, Z, pos, n_extra_data, extra_data, cell, Fhat)
Fhat = Fhat / sqrt(sum(Fhat*Fhat))
d_pos = d_pos - 2.0*Fhat*sum(Fhat*d_pos)
n_reflect = n_reflect + 1
if (no_reverse == 0) then ! step and reverse if needed
pos = pos + d_pos
E = ll_eval_energy(N, Z, pos, n_extra_data, extra_data, cell)
if (isnan(E) .or. E >= Emax) then ! reflect or reverse
pos = last_good_pos
d_pos = - last_good_d_pos
! undo reflection that just happened, replace it with a reversal
n_reflect = n_reflect - 1
n_reverse = n_reverse + 1
endif
endif
endif
end do
if (no_reverse /= 0) then
n_try = 1
! accept/reject
if (E < Emax) then ! accept
final_E = E
n_accept = 1
else ! reject (E=NaN should also end up here) and reverse dir
pos = last_good_pos
d_pos = -last_good_d_pos
final_E = last_good_E
n_accept = 0
endif
else
if (n_reverse > 0) E = ll_eval_energy(N, Z, pos, n_extra_data, extra_data, cell)
! we'll use the ration n_reflect/(n_reflect+n_reverse) as the tuning parameter
n_try = n_reflect + n_reverse
n_accept = n_reflect
final_E = E
endif
end subroutine fortran_GMC_atom
subroutine fortran_MD_atom_NVE(N, Z, pos, vel, mass, n_extra_data, extra_data, cell, n_steps, timestep, final_E, debug)
implicit none
integer :: N
integer :: Z(N)
double precision :: pos(3,N), vel(3,N), mass(N), cell(3,3)
integer :: n_extra_data
double precision :: extra_data(n_extra_data,N)
integer :: n_steps
double precision :: timestep, final_E
integer :: debug
integer i_step, i
double precision :: forces(3,N), acc(3,N), PE
double precision, external :: ll_eval_forces, ll_eval_energy
! initialize accelerations
PE = ll_eval_forces(N, Z, pos, n_extra_data, extra_data, cell, forces)
do i=1, 3
acc(i,:) = forces(i,:) / mass(:)
end do
if (debug > 0) print *, "initial PE KE E ", PE, 0.5*sum(spread(mass,1,3)*vel**2), &
PE+0.5*sum(spread(mass,1,3)*vel**2)
do i_step=1, n_steps
! Verlet part 1
vel = vel + 0.5*timestep*acc
pos = pos + timestep*vel
! new accelerations at t+dt
PE = ll_eval_forces(N, Z, pos, n_extra_data, extra_data, cell, forces)
do i=1, 3
acc(i,:) = forces(i,:) / mass(:)
end do
! Verlet part 2
vel = vel + 0.5*timestep*acc
if (debug > 0) print *, "step PE KE E ", i_step, PE, 0.5*sum(spread(mass,1,3)*vel**2), &
PE+0.5*sum(spread(mass,1,3)*vel**2)
end do
final_E = PE + 0.5*sum(spread(mass,1,3)*vel**2)
end subroutine fortran_MD_atom_NVE