-
Notifications
You must be signed in to change notification settings - Fork 0
/
MD.f90
458 lines (354 loc) · 12.8 KB
/
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
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
!-------------------------------------------------------------------------------------
! Copyright (c) 2016 Daniel C. Elton
!
! This software is licensed under The MIT License (MIT)
! Permission is hereby granted, free of charge, to any person obtaining a copy of this
! software and associated documentation files (the "Software"), to deal in the Software
! without restriction, including without limitation the rights to use, copy, modify, merge,
! publish, distribute, sublicense, and/or sell copies of the Software, and to permit
! persons to whom the Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
! BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
! DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!-------------------------------------------------------------------------------------
module MD
use dans_timer
use NormalModes, only: EvolveRing
use pot_mod
use dans_timer
use system_mod
use consts
use InputOutput
use force_calc
use main_stuff
contains
!-----------------------------------------------------------------
!-------- Standard full PIMD Velocity-Verlet integration ------
!-----------------------------------------------------------------
subroutine PIMD_VV
use mpi
implicit none
if (pid .eq. 0) then
!Propagate NH chains
if (BEADTHERMOSTAT) call bead_thermostat
if (THERMOSTAT) then
call Nose_Hoover(s, uk, global_chain_length, vxi_global, tau, delt2, 3*Natoms*Nbeads, temp*Nbeads)
PPt = PPt*s
endif
!update momenta a half step w/ old forces
PPt = PPt - MASSCON*dRRt*delt2
!Normal modes stuff
if (Nbeads .gt. 1) then
do i = 1, Nwaters
!Evolve ring a full step
Call EvolveRing(RRt(:,3*i-2,:), PPt(:,3*i-2,:), Nbeads, massO)
Call EvolveRing(RRt(:,3*i-1,:), PPt(:,3*i-1,:), Nbeads, massH)
Call EvolveRing(RRt(:,3*i-0,:), PPt(:,3*i-0,:), Nbeads, massH)
enddo
else
do i = 1,Nwaters
do k = 1,Nbeads
RRt(:,3*i-2,k) = RRt(:,3*i-2,k) + imassO*PPt(:,3*i-2,k)*delt
RRt(:,3*i-1,k) = RRt(:,3*i-1,k) + imassH*PPt(:,3*i-1,k)*delt
RRt(:,3*i-0,k) = RRt(:,3*i-0,k) + imassH*PPt(:,3*i-0,k)*delt
enddo
enddo
endif
!calculate centroid positions
RRc = sum(RRt,3)/Nbeads
!calculate centroid momenta
PPc = sum(PPt,3)/Nbeads
!check PBCs
call PBCs(RRt, RRc)
endif!if (pid .eq. 0)
!call force routine
call full_bead_forces
if (pid .eq. 0) then
!update kinetic energy
call calc_uk
!write stuff out if necessary
call start_timer("WritingOut")
call write_out
call stop_timer("WritingOut")
!update momenta a half step w/ new forces
PPt = PPt - MASSCON*dRRt*delt2
!Propagate NH chains
if (BEADTHERMOSTAT) call bead_thermostat
call calc_uk
if (THERMOSTAT) then
call Nose_Hoover(s, uk, global_chain_length, vxi_global, tau, delt2, 3*Natoms*Nbeads, temp*Nbeads)
PPt = PPt*s
endif
if (BAROSTAT) call Pcouple
endif!(pid .eq. 0)
end subroutine PIMD_VV
!-----------------------------------------------------------------
!-------- " Monomer PIMD" --------------------------------------
! THIS IS DEPRECATED --- EVERYTHING IS DONE WITH CONTRACTED ROUTINE BELOW
!-----------------------------------------------------------------
subroutine monomer_PIMD
use system_mod
use consts
use main_stuff
use InputOutput
use NormalModes
use mpi
use force_calc
implicit none
double precision :: e1, virialmon, virialcmon
double precision, dimension(3,Natoms,Nbeads) :: dRRmon
double precision, dimension(3,Natoms) :: dRRc
double precision, dimension(3,3) :: dr1, r1
double precision, dimension(3) :: roh1, roh2
if (t .eq. 1) dRRmon = 0
Upott = 0
virialmon = 0
virialcmon = 0
if (pid .eq. 0) then
!Propagate NH chains
if (BEADTHERMOSTAT) call bead_thermostat
if (THERMOSTAT) then
call Nose_Hoover(s, uk, global_chain_length, vxi_global, tau, delt2, 3*Natoms*Nbeads, temp*Nbeads)
PPt = PPt*s
endif
!update momenta a half step w/ old forces
PPt = PPt - MASSCON*dRRt*delt2
!Normal modes stuff
if (Nbeads .gt. 1) then
do i = 1, Nwaters
!Evolve ring a full step
call EvolveRing(RRt(:,3*i-2,:), PPt(:,3*i-2,:), Nbeads, massO)
call EvolveRing(RRt(:,3*i-1,:), PPt(:,3*i-1,:), Nbeads, massH)
call EvolveRing(RRt(:,3*i-0,:), PPt(:,3*i-0,:), Nbeads, massH)
enddo
else
do i = 1,Nwaters
do k = 1,Nbeads
RRt(:,3*i-2,k) = RRt(:,3*i-2,k) + imassO*PPt(:,3*i-2,k)*delt
RRt(:,3*i-1,k) = RRt(:,3*i-1,k) + imassH*PPt(:,3*i-1,k)*delt
RRt(:,3*i-0,k) = RRt(:,3*i-0,k) + imassH*PPt(:,3*i-0,k)*delt
enddo
enddo
endif
!calculate centroid positions
RRc = sum(RRt,3)/Nbeads
!calculate centroid momenta
PPc = sum(PPt,3)/Nbeads
!check PBCs
call PBCs(RRt, RRc)
!intermolecular force calculation
call potential(RRc, RRc, Upot, dRRc, virt, virialc, dip_momI, dip_momE, chg, t, BAROSTAT, sys_label)
Upott(:) = Upot !potential energy for the ENTIRE system (all images)
!intramolecular force calculation
do j = 1, Nbeads
do iw = 1, Nwaters
iO=3*iw-2; iH1 = 3*iw-1; iH2=3*iw-0
!PS potential energy surface calculation
r1(1:3, 1:3) = RRt(1:3, (/iO, iH1, iH2/), j)
if (SIESTA_MON_CALC) then
call siesta_monomer(r1, dr1, e1)
else
call pot_nasa(r1, dr1, e1, box, boxi)
endif
dRRmon(1:3, (/iO, iH1, iH2/), j) = dr1
Upott(j) = Upott(j) + e1
!monomer centroid virial
roh1 = RRc(1:3, iH1) - RRc(1:3, iO)
roh1 = roh1 - box*anint(roh1*boxi) !PBC
roh2 = RRc(1:3, iH2) - RRc(1:3, iO)
roh2 = roh2 - box*anint(roh2*boxi) !PBC
virialcmon = virialcmon + dot_product(roh1, dr1(:,2))
virialcmon = virialcmon + dot_product(roh2, dr1(:,3))
!monomer normal virial
roh1 = RRt(1:3, iH1, j) - RRt(1:3, iO, j)
roh1 = roh1 - box*anint(roh1*boxi) !PBC
roh2 = RRt(1:3, iH2, j) - RRt(1:3, iO, j)
roh2 = roh2 - box*anint(roh2*boxi) !PBC
virialmon = virialmon + dot_product(roh1, dr1(:,2))
virialmon = virialmon + dot_product(roh2, dr1(:,3))
enddo
enddo
!update dRRt
do j = 1, Nbeads
dRRt(:,:,j) = dRRc + dRRmon(:,:,j)
enddo
!update Upot
Upot = sum(Upott) !potential energy for the ENTIRE system
virial = virialmon + virt(1,1) + virt(2,2) + virt(3,3) !virial for the ENTIRE system (all images)
virialc = virialcmon/Nbeads + virialc
call calc_monomer_dip_moments(dip_momIt, RRt)
!calculate electronic polarization dipoles using TTM method
if (pot_model .eq. 6) call dip_ttm(RRc, dip_momI, dip_momE, chg, t)
!add electronic polarization dipoles to monomer dipoles
do j = 1, Nbeads
dip_momIt(:,:,j) = dip_momIt(:,:,j) + dip_momE
dip_momEt(:,:,j) = dip_momE
enddo
!update kinetic energy
call calc_uk
!write stuff out if necessary
call start_timer("WritingOut")
call write_out
call stop_timer("WritingOut")
!update momenta a half step w/ new forces
PPt = PPt - MASSCON*dRRt*delt2
!Propagate NH chains
if (BEADTHERMOSTAT) call bead_thermostat
call calc_uk
if (THERMOSTAT) then
call Nose_Hoover(s, uk, global_chain_length, vxi_global, tau, delt2, 3*Natoms*Nbeads, temp*Nbeads)
PPt = PPt*s
endif
if (BAROSTAT) call Pcouple
endif!(pid .eq. 0)
end subroutine monomer_PIMD
!--------------------------------------------------------------------------------------------
!------------- Contracted MD with intermolecular forces on monomer with multiple time step
!-- This subroutine performs evaluates only the intramolecular (fast) forces on the beads and
!-- evaluates the intermolecular (slow) forces on the centroid. It also uses a multiple timestep method
!-- (see Tuckerman, pg 118). The momentum and positions will be updated with the intramolecular forces
!-- every intra_timesteps times. For instance if the 'outer timestep' (normal timestep) is .5 ps and intra_timesteps = 5
!-- then the inner timestep is .1 ps
!--------------------------------------------------------------------------------------------
subroutine contracted_MD
use main_stuff
use NormalModes
use pot_mod
use dans_timer
use system_mod
use consts
use InputOutput
use force_calc
Implicit None
double precision :: e1, virialmon, virialcmon
double precision, dimension(3,Natoms,Nbeads) :: dRRfast
double precision, dimension(3,Natoms) :: dRRc
double precision, dimension(3,3) :: dr1, r1
double precision, dimension(3) :: roh1, roh2
integer :: tintra, iM
if (t .eq. 1) dRRfast =0
Upott = 0
virialmon = 0
virialcmon = 0
if (pid .eq. 0) then
!Propagate NH chains
if (BEADTHERMOSTAT) call bead_thermostat
if (THERMOSTAT) then
call Nose_Hoover(s, uk, global_chain_length, vxi_global, tau, delt2, 3*Natoms*Nbeads, temp*Nbeads)
PPt = PPt*s
endif
!update momenta a half step w/ old forces
PPt = PPt - MASSCON*dRRt*delt2
!calculate centroid positions
RRc = sum(RRt,3)/Nbeads
!calculate centroid momenta
PPc = sum(PPt,3)/Nbeads
!check PBCs
call PBCs(RRt, RRc)
call start_timer("MonomerPIMD")
!--- intramolecular (fast) forces -------------------------------------------------
do tintra = 1, intra_timesteps
!update momenta with fast forces (delt2fast = deltfast/2)
PPt = PPt - MASSCON*dRRfast*delt2fast
!update positions with fast forces
if (Nbeads .gt. 1) then
do i = 1, Nwaters
Call EvolveRing(RRt(:,3*i-2,:), PPt(:,3*i-2,:), Nbeads, massO)
Call EvolveRing(RRt(:,3*i-1,:), PPt(:,3*i-1,:), Nbeads, massH)
Call EvolveRing(RRt(:,3*i-0,:), PPt(:,3*i-0,:), Nbeads, massH)
enddo
else
do i = 1,Nwaters
do k = 1,Nbeads
RRt(:,3*i-2,k) = RRt(:,3*i-2,k) + imassO*PPt(:,3*i-2,k)*deltfast
RRt(:,3*i-1,k) = RRt(:,3*i-1,k) + imassH*PPt(:,3*i-1,k)*deltfast
RRt(:,3*i-0,k) = RRt(:,3*i-0,k) + imassH*PPt(:,3*i-0,k)*deltfast
enddo
enddo
endif
!update fast forces (intramolecular forces)
!masternode calcuates the intramolecular forces, puts them in dRRfast
virialmon = 0
virialcmon = 0
do j = 1, Nbeads
do iw = 1, Nwaters
iO=3*iw-2; iH1 = 3*iw-1; iH2=3*iw-0
r1(1:3, 1:3) = RRt(1:3, (/iO, iH1, iH2/), j)
call pot_nasa(r1, dr1, e1, box, boxi)
dRRfast(1:3, (/iO, iH1, iH2/), j) = dr1
!if last timestep in loop update monomer energy
!and calculate dipole moments using dip. mom. surface
if (tintra .eq. intra_timesteps) then
Upott(j) = Upott(j) + e1
!monomer centroid virial first
roh1 = RRc(1:3, iH1) - RRc(1:3, iO)
roh1 = roh1 - box*anint(roh1*boxi) !PBC
roh2 = RRc(1:3, iH2) - RRc(1:3, iO)
roh2 = roh2 - box*anint(roh2*boxi) !PBC
virialcmon = virialcmon + dot_product(roh1, dr1(:,2))
virialcmon = virialcmon + dot_product(roh2, dr1(:,3))
!monomer normal virial
roh1 = RRt(1:3, iH1, j) - RRt(1:3, iO, j)
roh1 = roh1 - box*anint(roh1*boxi) !PBC
roh2 = RRt(1:3, iH2, j) - RRt(1:3, iO, j)
roh2 = roh2 - box*anint(roh2*boxi) !PBC
virialmon = virialmon + dot_product(roh1, dr1(:,2))
virialmon = virialmon + dot_product(roh2, dr1(:,3))
endif
enddo
enddo
!update momenta with fast forces
PPt = PPt - MASSCON*dRRfast*delt2fast
enddo!tintra = 1..
!--- end intramolecular (fast) forces ----------------------------------------------
!calculate centroid positions
RRc = sum(RRt,3)/Nbeads
!calculate centroid momenta
PPc = sum(PPt,3)/Nbeads
!check PBCs
call PBCs(RRt, RRc)
call stop_timer("MonomerPIMD")
!intermolecular force calculation
call potential(RRc, RRc, Upot, dRRc, virt, virialc, dip_momI, dip_momE, chg, t, BAROSTAT, sys_label)
Upott(:) = Upott(:) + Upot !potential energy for the ENTIRE system (all images)
!update dRRt
do j = 1, Nbeads
dRRt(:,:,j) = dRRc
enddo
call calc_monomer_dip_moments(dip_momIt, RRt)
! !calculate electronic polarization dipoles using TTM method
if (pot_model .eq. 6) call dip_ttm(RRc, dip_momE, chg, t)
!write(*,*) dip_momE
!add electronic polarization dipoles to monomer dipoles
do j = 1, Nbeads
dip_momIt(:,:,j) = dip_momIt(:,:,j) + dip_momE
dip_momEt(:,:,j) = dip_momE
enddo
!update Upot, virial and virialc
Upot = sum(Upott) !potential energy for the ENTIRE system
virial = virialmon + virt(1,1) + virt(2,2) + virt(3,3) !virial for the ENTIRE system (all images)
virialc = virialcmon/Nbeads + virialc
!update kinetic energy
call calc_uk
!write stuff out if necessary
call start_timer("WritingOut")
call write_out
call stop_timer("WritingOut")
!update momenta a half step w/ slow forces
PPt = PPt - MASSCON*dRRt*delt2
!Propagate NH chains
if (BEADTHERMOSTAT) call bead_thermostat
call calc_uk
if (THERMOSTAT) then
call Nose_Hoover(s, uk, global_chain_length, vxi_global, tau, delt2, 3*Natoms*Nbeads, temp*Nbeads)
PPt = PPt*s
endif
if (BAROSTAT) call Pcouple
endif!(pid .eq. 0)
end subroutine contracted_MD
end module MD