-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathprinte.F90
537 lines (481 loc) · 17.6 KB
/
printe.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
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
#include "copyright.h"
#include "../include/dprec.fh"
#ifndef PBSA
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Emit the final minimization report
!-----------------------------------------------------------------------
! --- REPORT_MIN_RESULTS ---
!-----------------------------------------------------------------------
! Find the maximum component of the gradient (grdmax),
! emit this and other gradient details (printe),
! optionally do pseudocontact shift constraints,
! optionally decompose energies for mm_pbsa,
! and emit nmr related information (nmrptx).
subroutine report_min_results( nstep, gradient_rms, coordinates, &
forces, energies, igraph, xx, ix, ih )
use state
use decomp, only : checkdec, printdec, printpdec
use file_io_dat
implicit none
integer, intent(in) :: nstep
_REAL_, intent(in) :: gradient_rms
_REAL_, intent(in) :: coordinates(*)
_REAL_, intent(in) :: forces(*)
type(state_rec), intent(in) :: energies
character(len=4), intent(in) :: igraph(*) ! atom name map
_REAL_, intent(inout) :: xx(*) ! real dynamic memory
integer, intent(inout) :: ix(*) ! integer dynamic memory
character(len=4), intent(inout) :: ih(*) ! hollerith dynamic memory
# include "box.h"
# include "extra.h"
# include "../include/md.h"
# include "../include/memory.h"
# include "nmr.h"
# include "tgtmd.h"
# include "multitmd.h"
character(len=4) :: atom_name_of_gmax
integer :: atom_number_of_gmax
_REAL_ :: gradient_max
_REAL_ emtmd
if (master) then
write(6, '(/ /20x,a,/)' ) 'FINAL RESULTS'
call grdmax( forces, gradient_max, atom_number_of_gmax )
atom_name_of_gmax = igraph(atom_number_of_gmax)
if (imin /= 5) rewind(MDINFO_UNIT)
call printe( nstep, gradient_rms, gradient_max, energies, &
atom_number_of_gmax, atom_name_of_gmax )
if (idecomp > 0) call checkdec(idecomp)
if (idecomp == 1 .or. idecomp == 2) call printdec(ix)
if (idecomp == 3 .or. idecomp == 4) call printpdec(ix)
if (nmropt > 0) then
if (iredir(7) /= 0) call pcshift(-1,coordinates,forces)
call nmrptx(6)
call nmrptx(MDINFO_UNIT)
call ndvptx(coordinates,forces,ih(m04),ih(m02),ix(i02),nres, &
xx(l95),natom,xx(lwinv),xx(lnmr01),ix(inmr02),6)
end if
if (itgtmd == 2) then
emtmd = 0.0d0
call mtmdcall(emtmd,xx(lmtmd01),ix(imtmd02),coordinates,forces,ih(m04),ih(m02),ix(i02),&
ih(m06),xx(lmass),natom,nres,'PRNT')
end if
if (imin /= 5) call amflsh(MDINFO_UNIT)
end if
return
end subroutine report_min_results
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Emit a minimization progress report
!-----------------------------------------------------------------------
! --- REPORT_MIN_PROGRESS ---
!-----------------------------------------------------------------------
! Find the maximum component of the gradient (grdmax),
! emit this and other gradient details (printe),
! and emit nmr related information (nmrptx).
subroutine report_min_progress( nstep, gradient_rms, forces, energies, &
igraph, charge )
use state
use file_io_dat
use crg_reloc, only : ifcr, crprintcharges, cr_print_charge
implicit none
integer, intent(in) :: nstep
_REAL_, intent(in) :: gradient_rms
_REAL_, intent(in) :: forces(*)
type(state_rec), intent(in) :: energies
character(len=4), intent(in) :: igraph(*) ! atom name map
_REAL_, intent(in) :: charge(*)
# include "extra.h"
# include "../include/md.h"
# include "nmr.h"
character(len=4) :: atom_name_of_gmax
integer :: atom_number_of_gmax
_REAL_ :: gradient_max
if (master) then
call grdmax( forces, gradient_max, atom_number_of_gmax )
atom_name_of_gmax = igraph(atom_number_of_gmax)
if (imin /= 5) rewind(MDINFO_UNIT)
call printe( nstep, gradient_rms, gradient_max, energies, &
atom_number_of_gmax, atom_name_of_gmax )
if (nmropt > 0) then
call nmrptx(6)
call nmrptx(MDINFO_UNIT)
end if
if ( ifcr > 0 .and. crprintcharges > 0 ) then
call cr_print_charge( charge, nstep )
end if
if (imin /= 5) call amflsh(MDINFO_UNIT)
end if
return
end subroutine report_min_progress
#endif /*ifndef PBSA*/
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Compute the maximum gradient component and the corresponding atom
!-----------------------------------------------------------------------
! --- GRDMAX ---
!-----------------------------------------------------------------------
subroutine grdmax( gradient, gradient_max, atom_number_of_gmax )
use constants, only : zero
implicit none
_REAL_, intent(in) :: gradient(*)
! This is actually two-dimensional (3,natoms), but to enable
! vectorization on IA32 SSE platforms they are treated as
! one-dimensional; this may also improve software pipelining !
_REAL_, intent(out) :: gradient_max
integer, intent(out) :: atom_number_of_gmax
# include "../include/memory.h"
integer :: i
_REAL_ :: gi
gradient_max = ZERO
atom_number_of_gmax = 1
do i = 1,3*natom
gi = abs(gradient(i))
if (gi > gradient_max) then
gradient_max = gi
atom_number_of_gmax = i
end if
end do
atom_number_of_gmax = (atom_number_of_gmax - 1)/3 + 1
return
end subroutine grdmax
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ Print status of a minimization calculation.
!-----------------------------------------------------------------------
! --- PRINTE ---
!-----------------------------------------------------------------------
! Output the step number, the gradient rms, the max gradient component,
! and the atom label and atom number of the max gradient component.
subroutine printe( nstep, gradient_rms, gradient_max, ene, &
atom_number_of_gmax, atom_name_of_gmax )
#ifdef APBS
use file_io_dat
#endif
#ifdef RISMSANDER
use sander_rism_interface, only : rismprm, RISM_NONE, RISM_FULL, RISM_INTERP,&
rism_calc_type, rism_thermo_print
#endif
use qmmm_module, only : qmmm_nml
use cns_xref
use state ! Access to energy_rec
use charmm_mod, only : charmm_active
use crg_reloc, only : ifcr
use emap,only : temap,scemap
use ff11_mod, only : cmap_active
use sebomd_module, only : sebomd_obj
implicit none
integer, intent(in) :: nstep
_REAL_, intent(in) :: gradient_rms
_REAL_, intent(in) :: gradient_max
type(state_rec), intent(in) :: ene
integer, intent(in) :: atom_number_of_gmax
character(len=4), intent(in) :: atom_name_of_gmax
# include "../include/md.h"
# include "ew_mpole.h"
# include "ew_cntrl.h"
# include "tgtmd.h"
_REAL_ epot,enonb,enele,ehbond,ebond,eangle,edihed,enb14,eel14,egb,epb
_REAL_ econst,epolar,aveper,aveind,avetot,esurf,edisp,diprms,dipiter
_REAL_ dipole_temp,escf,dvdl,enemap
#ifdef RISMSANDER
_REAL_ :: erism
_REAL_ :: pot_array(potential_energy_rec_len)
#endif /*RISMSANDER*/
_REAL_ :: ect
! ----- Extract Energies. -----
epot = ene%pot%tot
enonb = ene%pot%vdw
enele = ene%pot%elec
ehbond = ene%pot%hbond
egb = ene%pot%gb
epb = ene%pot%pb
ebond = ene%pot%bond
eangle = ene%pot%angle
edihed = ene%pot%dihedral
enb14 = ene%pot%vdw_14
eel14 = ene%pot%elec_14
econst = ene%pot%constraint
epolar = ene%pot%polar
dvdl = ene%pot%dvdl
aveper = ene%aveper
aveind = ene%aveind
avetot = ene%avetot
esurf = ene%pot%surf
diprms = ene%diprms
dipiter = ene%dipiter
dipole_temp = ene%dipole_temp
escf = ene%pot%scf
edisp = ene%pot%disp
enemap = ene%pot%emap
#ifdef RISMSANDER
erism = ene%pot%rism
#endif /*RISMSANDER*/
ect = ene%pot%ct
write(6,9018)
write(6,9028) nstep, epot, gradient_rms, gradient_max, &
atom_name_of_gmax, atom_number_of_gmax
write(6,9038) ebond,eangle,edihed
! CHARMM SPECIFIC ENERGY TERMS
if (charmm_active) write(6,9039) ene%pot%angle_ub, &
ene%pot%imp, &
ene%pot%cmap
#ifdef RISMSANDER
if( igb == 0 .and. ipb == 0 .and. rismprm%rism == 0) then
#else
if( igb == 0 .and. ipb == 0 ) then
#endif
write(6,9048) enonb,enele,ehbond
else if ( igb == 10 .or. ipb /= 0 ) then
write(6,9050) enonb,enele,epb
#ifdef APBS
else if ( igb == 6 .and. mdin_apbs ) then
write(6,9050) enonb,enele,epb
#endif /* APBS */
#ifdef RISMSANDER
else if(rismprm%rism == 1 )then
write(6,9051) enonb,enele,erism
#endif
else
write(6,9049) enonb,enele,egb
end if
write(6,9058) enb14,eel14,econst
if ( ifcr /= 0 ) write(6,9099) ect
! wxw: EMAP energy
if (temap) then
write (6,9062) enemap,scemap
endif
if (qmmm_nml%ifqnt) then
!write the SCF energy
if (qmmm_nml%qmtheory%PM3) then
if (qmmm_nml%qmmm_int == 3) then
write(6,9090) escf ! PM3-MM*
else if (qmmm_nml%qmmm_int == 4) then
write(6,9096) escf ! PM3/MMX2
else
write(6,9080) escf
end if
else if (qmmm_nml%qmtheory%AM1) then
write(6,9081) escf
else if (qmmm_nml%qmtheory%AM1D) then
write(6,9981) escf
else if (qmmm_nml%qmtheory%MNDO) then
write(6,9082) escf
else if (qmmm_nml%qmtheory%MNDOD) then
write(6,9982) escf
else if (qmmm_nml%qmtheory%PDDGPM3) then
write(6,9083) escf
else if (qmmm_nml%qmtheory%PDDGMNDO) then
write(6,9084) escf
else if (qmmm_nml%qmtheory%PM3CARB1) then
write(6,9085) escf
else if (qmmm_nml%qmtheory%DFTB) then
write(6,9086) escf
else if (qmmm_nml%qmtheory%RM1) then
write(6,9087) escf
else if (qmmm_nml%qmtheory%PDDGPM3_08) then
write(6,9088) escf
else if (qmmm_nml%qmtheory%PM6) then
write(6,9089) escf
else if (qmmm_nml%qmtheory%PM3ZNB) then
write(6,9091) escf
else if (qmmm_nml%qmtheory%EXTERN) then
write(6,9092) escf
else if (qmmm_nml%qmtheory%PM3MAIS) then
write(6,9093) escf
else if (qmmm_nml%qmtheory%AM1DHFR) then
write(6,9094) escf
else
write(6,'(" ERROR - UNKNOWN QM THEORY")')
end if
end if
if (sebomd_obj%do_sebomd) then
write(6,9200) sebomd_obj%esebomd
end if
#ifdef PUPIL_SUPPORT
write(6,9900) escf
#endif
if( gbsa > 0 ) write(6,9077) esurf
if (igb == 10 .or. ipb /= 0) write(6,9074) esurf,edisp
call write_cns_xref_min_energies( ene )
#ifdef APBS
if (igb == 6 .and. mdin_apbs ) write(6,9069) esurf
#endif /* APBS */
if (cmap_active .and. ipol > 0 ) then
write(6,9066) epolar, ene%pot%cmap
else
if (cmap_active) write(6,9067) ene%pot%cmap
if (epolar /= 0.0) write(6,9068) epolar
end if
if (econst /= 0.0) write(6,9078) epot-econst
if ( dvdl /= 0.d0) write(6,9100) dvdl
if (induced > 0.and.indmeth < 3) write(6,9190)diprms,dipiter
if (induced > 0.and.indmeth == 3) write(6,9191)diprms, &
dipole_temp
if (itgtmd == 1) then
write (6,'(a,f8.3)') "Current RMSD from reference: ",rmsdvalue
write (6,'(a,f8.3)') "Current target RMSD: ",tgtrmsd
end if
call amflsh(6)
! ----- SEND IT TO THE INFO FILE -----
if (imin /= 5) then
write(7,9018)
write(7,9028) nstep, epot, gradient_rms, gradient_max, &
atom_name_of_gmax, atom_number_of_gmax
write(7,9038) ebond,eangle,edihed
! CHARMM SPECIFIC ENERGY TERMS
if (charmm_active) write(7,9039) ene%pot%angle_ub, &
ene%pot%imp, &
ene%pot%cmap
#ifdef RISMSANDER
if( igb == 0 .and. ipb == 0 .and. rismprm%rism == 0) then
#else
if( igb == 0 .and. ipb == 0 ) then
#endif
write(7,9048) enonb,enele,ehbond
else if ( igb == 10 .or. ipb /= 0 ) then
write(7,9050) enonb,enele,epb
#ifdef APBS
else if ( igb == 6 .and. mdin_apbs ) then
write(7,9050) enonb,enele,epb
#endif /* APBS */
#ifdef RISMSANDER
else if ( rismprm%rism == 1 ) then
write(7,9051) enonb,enele,erism
#endif
else
write(7,9049) enonb,enele,egb
end if
write(7,9058) enb14,eel14,econst
! wxw: EMAP energy
if (temap) then
write (7,9062) enemap,scemap
endif
if (qmmm_nml%ifqnt) then
!write the SCF energy
if (qmmm_nml%qmtheory%PM3) then
if (qmmm_nml%qmmm_int == 3) then
write(7,9090) escf ! PM3-MM*
else if (qmmm_nml%qmmm_int == 4) then
write(7,9096) escf ! PM3/MMX2
else
write(7,9080) escf
end if
else if (qmmm_nml%qmtheory%AM1) then
write(7,9081) escf
else if (qmmm_nml%qmtheory%AM1D) then
write(7,9981) escf
else if (qmmm_nml%qmtheory%MNDO) then
write(7,9082) escf
else if (qmmm_nml%qmtheory%MNDOD) then
write(7,9982) escf
else if (qmmm_nml%qmtheory%PDDGPM3) then
write(7,9083) escf
else if (qmmm_nml%qmtheory%PDDGMNDO) then
write(7,9084) escf
else if (qmmm_nml%qmtheory%PM3CARB1) then
write(7,9085) escf
else if (qmmm_nml%qmtheory%DFTB) then
write(7,9086) escf
else if (qmmm_nml%qmtheory%RM1) then
write(7,9087) escf
else if (qmmm_nml%qmtheory%PDDGPM3_08) then
write(7,9088) escf
else if (qmmm_nml%qmtheory%PM6) then
write(7,9089) escf
else if (qmmm_nml%qmtheory%PM3ZNB) then
write(7,9091) escf
else if (qmmm_nml%qmtheory%EXTERN) then
write(7,9092) escf
else if (qmmm_nml%qmtheory%PM3MAIS) then
write(7,9093) escf
else if (qmmm_nml%qmtheory%AM1DHFR) then
write(7,9094) escf
else
write(7,'(" ERROR - UNKNOWN QM THEORY")')
end if
end if
if (sebomd_obj%do_sebomd) then
write(7,9200) sebomd_obj%esebomd
end if
#ifdef PUPIL_SUPPORT
write(7,9900) escf
#endif
if( gbsa > 0 ) write(7,9077) esurf
if ( igb == 10 .or. ipb /= 0 ) write(7,9074) esurf,edisp
#ifdef APBS
if (igb == 6 .and. mdin_apbs ) write(7,9069) esurf
#endif /* APBS */
! FF11 CMAP SPECIFIC ENERGY TERMS
if (cmap_active .and. epolar /= 0.0 ) then
write(7,9066) epolar, ene%pot%cmap
else
if (cmap_active) write(7,9067) ene%pot%cmap
if (epolar /= 0.0) write(7,9068) epolar
end if
if (econst /= 0.0) write(7,9078) epot-econst
if ( dvdl /= 0.d0) write(7,9100) dvdl
if (induced > 0.and.indmeth < 3) write(7,9190)diprms,dipiter
if (induced > 0.and.indmeth == 3) write(7,9191)diprms, &
dipole_temp
end if
#ifdef RISMSANDER
if(rismprm%rism==1 .and. rismprm%write_thermo==1)then
if(rism_calc_type(nstep) == RISM_FULL)&
call rism_thermo_print(.false.,transfer(ene%pot,pot_array))
end if
#endif /*RISMSANDER*/
9018 format (/ /,3x,'NSTEP',7x,'ENERGY',10x,'RMS',12x,'GMAX',9x, &
'NAME',4x,'NUMBER')
9028 format(1x,i6,2x,3(2x,1pe13.4),5x,a4,2x,i7,/)
9038 format (1x,'BOND = ',f13.4,2x,'ANGLE = ',f13.4,2x, &
'DIHED = ',f13.4)
9039 format (1x, 'UB = ', f13.4, 2x, 'IMP = ', f13.4, 2x, &
'CMAP = ', f13.4)
9048 format (1x,'VDWAALS = ',f13.4,2x,'EEL = ',f13.4,2x, &
'HBOND = ',f13.4)
9049 format (1x,'VDWAALS = ',f13.4,2x,'EEL = ',f13.4,2x, &
'EGB = ',f13.4)
9050 format (1x,'VDWAALS = ',f13.4,2x,'EEL = ',f13.4,2x, &
'EPB = ',f13.4)
#ifdef RISMSANDER
9051 format (1x,'VDWAALS = ',f13.4,2x,'EEL = ',f13.4,2x, &
'ERISM = ',f13.4)
#endif
9058 format (1x,'1-4 VDW = ',f13.4,2x,'1-4 EEL = ',f13.4,2x, &
'RESTRAINT = ',f13.4)
9062 format (1x,'EMAP = ',f14.4, ' EMSCORE = ',f8.4)
9066 format (1x,'EPOLAR = ',f13.4,2x,'CMAP = ',f13.4)
9067 format (1x,'CMAP = ',f13.4)
9068 format (1x,'EPOLAR = ',f13.4)
#ifdef APBS
9069 format (1x,'ENPOLAR = ',f13.4)
#endif /* APBS */
9074 format (1x,'ECAVITY = ',f13.4,2x,'EDISPER = ',f13.4)
9077 format (1x,'ESURF = ',f13.4)
9078 format (1x,'EAMBER = ',f13.4)
9080 format (1x,'PM3ESCF =',f14.4)
9081 format (1x,'AM1ESCF =',f14.4)
9981 format (1x,'AM1DESCF =',f14.4)
9082 format (1x,'MNDOESCF=',f14.4)
9982 format (1x,'MNDODESCF=',f14.4)
9083 format (1x,'PDDGPM3-ESCF=',f14.4)
9084 format (1x,'PDDGMNDO-ESCF=',f14.4)
9085 format (1x,'PM3CARB1-ESCF=',f14.4)
9086 format (1x,'DFTBESCF=',f14.4)
9087 format (1x,'RM1ESCF =',f14.4)
9088 format (1x,'PDDGPM3_08-ESCF=',f14.4)
9089 format (1x,'PM6ESCF =',f14.4)
9090 format (1x,'PM3MMXESCF =',f14.4)
9091 format (1x,'PM3ZNBESCF =',f14.4)
9092 format (1x,'EXTERNESCF =',f14.4)
9093 format (1x,'PM3MAISESCF =',f14.4)
9094 format (1x,'AM1DHFRESCF =',f14.4)
9099 format (1x,'ECRG =',f14.4)
9096 format (1x,'PM3MMX2ESCF =',f14.4)
9190 format(1x,'Dipole convergence: rms = ',e10.3,' iters = ',f6.2)
9191 format(1x,'Dipole convergence: rms = ',e10.3, &
' temperature = ',f6.2)
9100 format (1x,'DV/DL = ',f14.4)
9200 format (1x,'ESEBOMD =',f14.4)
#ifdef PUPIL_SUPPORT
9900 format (1x,'PUPESCF =',f14.4)
#endif
return
end subroutine printe