forked from lgcrego/Dynemol
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Berendsen.f
319 lines (250 loc) · 8.73 KB
/
Berendsen.f
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
module Berendsen_m
use constants_m
use syst ! using all syst
use parameters_m , only: PBC
use MD_read_m , only: MM , atom , molecule, species
use for_force , only: forcefield
use f_inter_m , only: stressr , stresre
use VV_Parent , only: VV
public :: Berendsen , VirialKinetic , barostat
type, extends(VV) :: Berendsen
contains
procedure :: VV1
procedure :: VV2
end type Berendsen
interface Berendsen
module procedure constructor
end interface
! module variables ...
real*8 :: stresvv(3,3)
contains
!
!
!
!===================================
function constructor() result( me )
!===================================
implicit none
type(Berendsen) :: me
!local variable ...
! select atomic or molecular kinetic energy to calculate the temperature ...
me % thermostat_type = NINT( float(maxval(molecule%N_of_atoms)) / float(MM%N_of_molecules) )
end function constructor
!
!
!
!========================
subroutine VV1( me , dt )
!========================
implicit none
class(Berendsen) , intent(inout) :: me
real*8 , intent(in) :: dt
! local variables ...
real*8 :: ai(3)
real*8 :: massa , dt_HALF , dt2_HALF
integer :: i
dt_HALF = dt / two
dt2_HALF = dt_HALF * dt
! VV1 ...
do i = 1 , MM % N_of_atoms
if( atom(i) % flex ) then
massa = mol / atom(i) % mass
ai = atom(i) % ftotal * massa
atom(i) % xyz = atom(i) % xyz + ( atom(i) % vel*dt + dt2_HALF*ai ) * mts_2_Angs
atom(i) % vel = atom(i) % vel + dt_HALF*ai
end if
end do
end subroutine VV1
!
!
!
!=========================
subroutine VV2( me , dt )
!=========================
implicit none
class(Berendsen) , intent(inout) :: me
real*8 , intent(in) :: dt
! local variables ...
integer :: i , j , j1 , j2 , nresid
real*8 :: massa , factor , sumtemp , temp , lambda , dt_half
real*8 :: tmp(3) , V_CM(3) , V_atomic(3)
CALL VirialKinetic( me % thermostat_type , sumtemp )
!######################################################
! Berendsen Thermostat
if( sumtemp == 0.d0 ) then
temp = D_ONE
else
! instantaneous temperature : E_kin/(3/2*NkB) ...
select case (me % thermostat_type)
case (0:1) ! <== molecular ...
temp = sumtemp * iboltz / real( count(molecule%flex) )
case (2:) ! <== atomic ...
temp = sumtemp * iboltz / real( count(atom%flex) )
end select
endif
! Berendsen Thermostat ; turned off for talt == infty ...
If( talt == infty ) then
lambda = D_one
else
lambda = ( dt / (talt*pico_2_sec) ) * ( bath_T / temp - D_ONE )
lambda = SQRT(D_ONE + lambda)
end If
!######################################################
dt_half = dt / two
! VV2 ...
sumtemp = D_zero
select case (me % thermostat_type)
case (0:1) ! <== molecular ...
do i = 1 , MM % N_of_molecules
tmp = D_zero
nresid = molecule(i) % nr
j1 = sum(molecule(1:nresid-1) % N_of_atoms) + 1
j2 = sum(molecule(1:nresid) % N_of_atoms)
do j = j1 , j2
if ( atom(j) % flex ) then
massa = mol / atom(j) % mass
! atom(j) % vel = atom(j) % vel * lambda + ( dt_half * atom(j) % ftotal ) * massa
atom(j) % vel = atom(j) % vel + ( dt_half * atom(j) % ftotal ) * massa
atom(j) % vel = atom(j) % vel * lambda
tmp = tmp + atom(j) % vel / massa
end if
end do
V_CM = tmp / molecule(i) % mass
sumtemp = sumtemp + molecule(i) % mass * sum( V_CM * V_CM )
end do
case (2:) ! <== atomic ...
V_atomic = D_zero
do i = 1 , MM % N_of_atoms
if( atom(i) % flex ) then
massa = mol / atom(i) % mass
! atom(i) % vel = atom(i) % vel * lambda + ( dt_half * atom(i) % ftotal ) * massa
atom(i) % vel = atom(i) % vel + ( dt_half * atom(i) % ftotal ) * massa
atom(i) % vel = atom(i) % vel * lambda
V_atomic = atom(i) % vel
factor = imol * atom(i) % mass
sumtemp = sumtemp + factor * sum( V_atomic * V_atomic )
end if
end do
end select
! instantaneous temperature of the system after contact with thermostat ...
select case (me % thermostat_type)
case (0:1) ! <== molecular ...
me % Temperature = sumtemp * iboltz / real( count(molecule%flex) )
case (2:) ! <atomic ...
me % Temperature = sumtemp * iboltz / real( count(atom%flex) )
end select
! calculation of the kinetic energy ...
me % Kinetic = D_zero
do i = 1 , MM % N_of_atoms
me % Kinetic = me % Kinetic + ( atom(i) % mass ) * sum( atom(i) % vel(:) * atom(i) % vel(:) ) * half
end do
me % Kinetic = me % Kinetic * micro / MM % N_of_Molecules
! calculation of pressure and density ...
CALL barostat( me , dt )
end subroutine VV2
!
!
!
!=====================================================
subroutine VirialKinetic( thermostat_type , sumtemp )
!=====================================================
implicit none
integer , intent(in) :: thermostat_type
real*8 , optional , intent(out) :: sumtemp
! local variables ...
real*8 :: tmp(3) , V_CM(3) , V_atomic(3) , factor , kinetic
integer :: i , j , j1 , j2 , k , nresid
! VV2 and thermostats ...
kinetic = D_zero
stresvv = D_zero
select case ( thermostat_type )
case(0:1) ! <== molecular ...
do i = 1 , MM % N_of_molecules
tmp = D_zero
nresid = molecule(i) % nr
j1 = sum(molecule(1:nresid-1) % N_of_atoms) + 1
j2 = sum(molecule(1:nresid) % N_of_atoms)
do j = j1 , j2
if( atom(j) % flex ) then
tmp = tmp + atom(j)%mass * atom(j)%vel
end if
end do
V_CM = tmp * imol / molecule(i) % mass
forall(k=1:3) stresvv(k,k) = stresvv(k,k) + molecule(i) % mass * V_CM(k) * V_CM(k)
stresvv(1,2) = stresvv(1,2) + molecule(i) % mass * V_CM(1) * V_CM(2)
stresvv(1,3) = stresvv(1,3) + molecule(i) % mass * V_CM(1) * V_CM(3)
stresvv(2,3) = stresvv(2,3) + molecule(i) % mass * V_CM(2) * V_CM(3)
! 2*kinetic energy (molecular) of the system ...
kinetic = kinetic + molecule(i) % mass * sum( V_CM * V_CM )
end do
case (2:) ! <== atomic ...
do i = 1 , MM % N_of_atoms
If( atom(i) % flex ) then
V_atomic = atom(i) % vel
factor = imol * atom(i) % mass
forall( k=1:3) stresvv(k,k) = stresvv(k,k) + factor * V_atomic(k) * V_atomic(k)
stresvv(1,2) = stresvv(1,2) + factor * V_atomic(1) * V_atomic(2)
stresvv(1,3) = stresvv(1,3) + factor * V_atomic(1) * V_atomic(3)
stresvv(2,3) = stresvv(2,3) + factor * V_atomic(2) * V_atomic(3)
! 2*kinetic energy (atomic) of the system ...
kinetic = kinetic + factor * sum( V_atomic * V_atomic )
end if
end do
end select
stresvv(2,1) = stresvv(1,2)
stresvv(3,1) = stresvv(1,3)
stresvv(3,2) = stresvv(2,3)
If( present(sumtemp) ) sumtemp = kinetic
end subroutine VirialKinetic
!
!
!
!==============================
subroutine barostat( me , dt )
!==============================
implicit none
class(Berendsen) , intent(inout) :: me
real*8 , intent(in) :: dt
! local variables ...
integer :: i, j, j1, j2, nresid
real*8 :: mip , volume, massa
real*8 :: Astres(3,3)
massa = sum( species(:) % N_of_molecules * species(:) % mass )
volume = product( MM % box(:) * Angs_2_mts )
me % density = (massa / volume) * milli
if (forcefield == 1) then
else
stresvv = stresvv / (cm_2_Angs * volume)
stressr = stressr / (cm_2_Angs * volume)
stresre = stresre / (cm_2_Angs * volume)
Astres = stresvv + stressr + stresre
endif
! pressure = instantaneous pressure ; press = external pressure ...
if (forcefield == 1) then
else
me % pressure = ( Astres(1,1) + Astres(2,2) + Astres(3,3) ) * third
endif
! Pressurestat ; turned off for talp == infty ...
If( talp == infty ) then
mip = D_one
else
mip = dt * ( 107.0d-6 / (talp * pico_2_sec) ) * ( me % pressure - press )
mip = (D_one + mip)**third
end If
MM % box = MM % box * mip
do i = 1 , MM % N_of_molecules
nresid = molecule(i) % nr
j1 = sum(molecule(1:nresid-1) % N_of_atoms) + 1
j2 = sum(molecule(1:nresid) % N_of_atoms)
do j = j1 , j2
if ( atom(j) % flex ) then
atom(j) % xyz = atom(j) % xyz * mip
atom(j) % xyz = atom(j) % xyz - MM % box * DNINT( atom(j) % xyz * MM % ibox ) * PBC
end if
end do
end do
end subroutine barostat
!
!
!
end module Berendsen_m