forked from TinkerTools/tinker
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmdstat.f
320 lines (320 loc) · 10.7 KB
/
mdstat.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
320
c
c
c ###################################################
c ## COPYRIGHT (C) 1990 by Jay William Ponder ##
c ## All Rights Reserved ##
c ###################################################
c
c #################################################################
c ## ##
c ## subroutine mdstat -- compute averages over a trajectory ##
c ## ##
c #################################################################
c
c
c "mdstat" is called at each molecular dynamics time step to
c form statistics on various average values and fluctuations,
c and to periodically save the state of the trajectory
c
c
subroutine mdstat (istep,dt,etot,epot,ekin,temp,pres)
use atoms
use bath
use bound
use boxes
use inform
use iounit
use limits
use mdstuf
use molcul
use units
use usage
use warp
implicit none
integer istep
integer modstep
real*8 dt,temp,pres
real*8 etot,epot,ekin
real*8 pico,dens
real*8 fluctuate,fluctuate2
real*8 potfluct,potfluct2
real*8 kinfluct,kinfluct2
real*8 tfluct,pfluct,dfluct
real*8 tfluct2,pfluct2,dfluct2
real*8 etot_sum,etot2_sum
real*8 etot_ave,etot2_ave
real*8 epot_sum,epot2_sum
real*8 epot_ave,epot2_ave
real*8 ekin_sum,ekin2_sum
real*8 ekin_ave,ekin2_ave
real*8 temp_sum,temp2_sum
real*8 temp_ave,temp2_ave
real*8 pres_sum,pres2_sum
real*8 pres_ave,pres2_ave
real*8 dens_sum,dens2_sum
real*8 dens_ave,dens2_ave
save etot_sum,etot2_sum
save epot_sum,epot2_sum
save ekin_sum,ekin2_sum
save temp_sum,temp2_sum
save pres_sum,pres2_sum
save dens_sum,dens2_sum
c
c
c set number of steps for block averages of properties
c
modstep = mod(istep,iprint)
c
c zero out summation variables for new averaging period
c
if (modstep.eq.1 .or. iprint.eq.1) then
etot_sum = 0.0d0
etot2_sum = 0.0d0
epot_sum = 0.0d0
epot2_sum = 0.0d0
ekin_sum = 0.0d0
ekin2_sum = 0.0d0
temp_sum = 0.0d0
temp2_sum = 0.0d0
pres_sum = 0.0d0
pres2_sum = 0.0d0
dens_sum = 0.0d0
dens2_sum = 0.0d0
end if
c
c print energy, temperature and pressure for current step
c
if (verbose) then
if (modstep .eq. 1) then
if (use_bounds .and. integrate.ne.'STOCHASTIC') then
write (iout,10)
10 format (/,4x,'MD Step',7x,'E Total',4x,'E Potential',
& 6x,'E Kinetic',7x,'Temp',7x,'Pres',/)
else
write (iout,20)
20 format (/,4x,'MD Step',7x,'E Total',4x,'E Potential',
& 6x,'E Kinetic',7x,'Temp',/)
end if
end if
if (use_bounds .and. integrate.ne.'STOCHASTIC') then
write (iout,30) istep,etot,epot,ekin,temp,pres
30 format (i10,3f15.4,2f11.2)
else
write (iout,40) istep,etot,epot,ekin,temp
40 format (i10,3f15.4,f11.2)
end if
flush (iout)
end if
c
c print header for the averages over a group of recent steps
c
if (modstep .eq. 0) then
pico = dble(istep) * dt
write (iout,50) iprint,istep
50 format (/,' Average Values for the Last',i6,' Out of',
& i9,' Dynamics Steps')
if (digits .ge. 8) then
write (iout,60) pico
60 format (/,' Simulation Time',5x,f19.8,' Picosecond')
else if (digits .ge. 6) then
write (iout,70) pico
70 format (/,' Simulation Time',5x,f17.6,' Picosecond')
else
write (iout,80) pico
80 format (/,' Simulation Time',5x,f15.4,' Picosecond')
end if
end if
c
c compute total energy and fluctuation for recent steps
c
etot_sum = etot_sum + etot
etot2_sum = etot2_sum + etot**2
if (modstep .eq. 0) then
etot_ave = etot_sum / dble(iprint)
etot2_ave = etot2_sum / dble(iprint)
fluctuate2 = etot2_ave - etot_ave**2
if (fluctuate2 .gt. 0.0d0) then
fluctuate = sqrt(fluctuate2)
else
fluctuate = 0.0d0
end if
if (digits .ge. 8) then
write (iout,90) etot_ave,fluctuate
90 format (' Total Energy',8x,f19.8,' Kcal/mole',3x,
& '(+/-',f13.8,')')
else if (digits .ge. 6) then
write (iout,100) etot_ave,fluctuate
100 format (' Total Energy',8x,f17.6,' Kcal/mole',3x,
& '(+/-',f11.6,')')
else
write (iout,110) etot_ave,fluctuate
110 format (' Total Energy',8x,f15.4,' Kcal/mole',3x,
& '(+/-',f9.4,')')
end if
end if
c
c compute average potential energy and its fluctuation
c
epot_sum = epot_sum + epot
epot2_sum = epot2_sum + epot**2
if (modstep .eq. 0) then
epot_ave = epot_sum / dble(iprint)
epot2_ave = epot2_sum / dble(iprint)
potfluct2 = epot2_ave - epot_ave**2
if (potfluct2 .gt. 0.0d0) then
potfluct = sqrt(potfluct2)
else
potfluct = 0.0d0
end if
if (digits .ge. 8) then
write (iout,120) epot_ave,potfluct
120 format (' Potential Energy',4x,f19.8,' Kcal/mole',3x,
& '(+/-',f13.8,')')
else if (digits .ge. 6) then
write (iout,130) epot_ave,potfluct
130 format (' Potential Energy',4x,f17.6,' Kcal/mole',3x,
& '(+/-',f11.6,')')
else
write (iout,140) epot_ave,potfluct
140 format (' Potential Energy',4x,f15.4,' Kcal/mole',3x,
& '(+/-',f9.4,')')
end if
end if
c
c compute average kinetic energy and its fluctuation
c
ekin_sum = ekin_sum + ekin
ekin2_sum = ekin2_sum + ekin**2
if (modstep .eq. 0) then
ekin_ave = ekin_sum / dble(iprint)
ekin2_ave = ekin2_sum / dble(iprint)
kinfluct2 = ekin2_ave - ekin_ave**2
if (kinfluct2 .gt. 0.0d0) then
kinfluct = sqrt(kinfluct2)
else
kinfluct = 0.0d0
end if
if (digits .ge. 8) then
write (iout,150) ekin_ave,kinfluct
150 format (' Kinetic Energy',6x,f19.8,' Kcal/mole',3x,
& '(+/-',f13.8,')')
else if (digits .ge. 6) then
write (iout,160) ekin_ave,kinfluct
160 format (' Kinetic Energy',6x,f17.6,' Kcal/mole',3x,
& '(+/-',f11.6,')')
else
write (iout,170) ekin_ave,kinfluct
170 format (' Kinetic Energy',6x,f15.4,' Kcal/mole',3x,
& '(+/-',f9.4,')')
end if
end if
c
c compute the average temperature and its fluctuation
c
temp_sum = temp_sum + temp
temp2_sum = temp2_sum + temp**2
if (modstep .eq. 0) then
temp_ave = temp_sum / dble(iprint)
temp2_ave = temp2_sum / dble(iprint)
tfluct2 = temp2_ave - temp_ave**2
if (tfluct2 .gt. 0.0d0) then
tfluct = sqrt(tfluct2)
else
tfluct = 0.0d0
end if
if (digits .ge. 8) then
write (iout,210) temp_ave,tfluct
210 format (' Temperature',9x,f19.6,' Kelvin',6x,
& '(+/-',f13.6,')')
else if (digits .ge. 6) then
write (iout,220) temp_ave,tfluct
220 format (' Temperature',9x,f17.4,' Kelvin',6x,
& '(+/-',f11.4,')')
else
write (iout,230) temp_ave,tfluct
230 format (' Temperature',9x,f15.2,' Kelvin',6x,
& '(+/-',f9.2,')')
end if
end if
c
c compute the average pressure and its fluctuation
c
if (use_bounds) then
pres_sum = pres_sum + pres
pres2_sum = pres2_sum + pres**2
if (modstep .eq. 0) then
pres_ave = pres_sum / dble(iprint)
pres2_ave = pres2_sum / dble(iprint)
pfluct2 = pres2_ave - pres_ave**2
if (pfluct2 .gt. 0.0d0) then
pfluct = sqrt(pfluct2)
else
pfluct = 0.0d0
end if
if (digits .ge. 8) then
write (iout,240) pres_ave,pfluct
240 format (' Pressure',12x,f19.6,' Atmosphere',2x,
& '(+/-',f13.6,')')
else if (digits .ge. 6) then
write (iout,250) pres_ave,pfluct
250 format (' Pressure',12x,f17.4,' Atmosphere',2x,
& '(+/-',f11.4,')')
else
write (iout,260) pres_ave,pfluct
260 format (' Pressure',12x,f15.2,' Atmosphere',2x,
& '(+/-',f9.2,')')
end if
end if
c
c compute the average density and its fluctuation
c
dens = (1.0d24/volbox) * (totmass/avogadro)
dens_sum = dens_sum + dens
dens2_sum = dens2_sum + dens**2
if (modstep .eq. 0) then
dens_ave = dens_sum / dble(iprint)
dens2_ave = dens2_sum / dble(iprint)
dfluct2 = dens2_ave - dens_ave**2
if (dfluct2 .gt. 0.0d0) then
dfluct = sqrt(dfluct2)
else
dfluct = 0.0d0
end if
if (digits .ge. 8) then
write (iout,270) dens_ave,dfluct
270 format (' Density',13x,f19.8,' Grams/cc',4x,
& '(+/-',f13.8,')')
else if (digits .ge. 6) then
write (iout,280) dens_ave,dfluct
280 format (' Density',13x,f17.6,' Grams/cc',4x,
& '(+/-',f11.6,')')
else
write (iout,290) dens_ave,dfluct
290 format (' Density',13x,f15.4,' Grams/cc',4x,
& '(+/-',f9.4,')')
end if
end if
end if
c
c declare deformation value for potential energy smoothing
c
if (use_smooth) then
if (modstep .eq. 0) then
if (digits .ge. 8) then
write (iout,300) deform
300 format (' Deformation',9x,f19.8,' Sqr Angs')
else if (digits .ge. 6) then
write (iout,310) deform
310 format (' Deformation',9x,f17.6,' Sqr Angs')
else
write (iout,320) deform
320 format (' Deformation',9x,f15.4,' Sqr Angs')
end if
end if
end if
c
c ensure any output is written to the storage device
c
if (modstep .eq. 0) flush (iout)
return
end