********> bugfix.31
Correction Authors: David A. Pearlman and James W. Caldwell
Correction Date: 2/25/92
Revision date: 2/28/92
Revision note: The revision is the additional removal of the DO-loop 898
contained in the context diff, below. Thanks to Julian Tirado-Rives
at Yale for pointing it out.
Programs: MINMD
Severity: serious
Note: This is a 2-part bug: fixing the first part revealed the second.
The first part of the fix is due to Pearlman and the second to Caldwell.
1st Problem: The velocities (and trajectory) calculated for a simulation
where a belly is employed (IBELLY=1) and periodic boundary conditions
are being used will likely be incorrect.
This problem will only occur when the following conditions are met
1) IBELLY=1 (belly simulation)
2) Periodic boundary conditions are being used (periodic box
set up in EDIT).
3) Using the program MINMD (problem does not occur in SANDER/GIBBS).
Affects: The reported and internal velocities will be incorrect. The
trajectory will also be incorrect. The simulation will likely die
within a small number of steps, due to the instabilities imparted
by the incorrect velocities.
Cause: One of the arguments in a call to clear the force array on each
step is incorrect. This argument, which should have been equal
to 3*number_of_atoms is, in fact, larger than that value. The result
is that not only is the force array cleared, but part of the velocity
array is inadvertantly cleared, as well.
2nd Problem: The 1st fix (to lib/force.f) unearthed a second bug, also
dependent on running periodic boundary simulations with a belly.
The symptom was failure of the ensemble to warm up properly because
of the improper setting of reference temperature (kinetic energy)
variables in an attempt to avoid divide by 0.
The fix for this is in the minmd/runmd.f section below.
Fix: Make the following changes.
Key for context diffs (see 0README for more details):
- (delete)
+ (add)
! (change)
--------------------------routine FORCE in src/lib/force.f
*** OLD force.f
--- NEW force.f
***************
*** 110,118
dimension xx(*),ix(*),ih(*)
DIMENSION X(*),F(*),ENE(30),ENER(*),VIR(*)
C
- NR = NPM*NRP+NSM*NRAM
- NR3 = 3*NR
- C
aveper = 0.0d0
aveind = 0.0d0
avetot = 0.0d0
--- 113,118 -----
dimension xx(*),ix(*),ih(*)
DIMENSION X(*),F(*),ENE(30),ENER(*),VIR(*)
C
aveper = 0.0d0
aveind = 0.0d0
avetot = 0.0d0
***************
*** 126,134
do 897 jn = 1,30
ene(jn) = 0.0d0
897 continue
- do 898 jn = 1,nr3
- f(jn) = 0.0d0
- 898 continue
c
DUM = 0.0D0
IMES = 0
--- 126,131 -----
do 897 jn = 1,30
ene(jn) = 0.0d0
897 continue
c
DUM = 0.0D0
IMES = 0
***************
*** 152,158
vt(3) = 0.0d0
vt(4) = 0.0d0
CALL ZEROUT(NREP+NREPC,ENE)
! CALL ZEROUT(NR3,F)
DUM = 0.0D+00
C
C ----- PARTITION THE WORK ARRAY FOR NNBOND -----
--- 149,155 -----
vt(3) = 0.0d0
vt(4) = 0.0d0
CALL ZEROUT(NREP+NREPC,ENE)
! CALL ZEROUT(3*NATOM,F)
DUM = 0.0D+00
C
C ----- PARTITION THE WORK ARRAY FOR NNBOND -----
--------------------------routine RUNMD in src/minmd/runmd.f
*** OLD runmd.f
--- New runmd.f
***************
*** 373,379
else
SCALP = SQRT(1.d0 + DTTP*(EKINP0/EKGP-1.d0))
endif
! scals = sqrt(1.d0 + dtts*(ekins0/ekgs-1.d0))
130 I3 = 0
c write(6,'(t2,''scalp,scals'',2f10.3)')scalp,scals
DO 133 J = 1,NR
--- 375,386 -----
else
SCALP = SQRT(1.d0 + DTTP*(EKINP0/EKGP-1.d0))
endif
! if(ekgs .lt. 1.d-06) then
! scals = 0.0d0
! else
! scals = sqrt(1.d0 + dtts*(ekins0/ekgs-1.d0))
! endif
!
130 I3 = 0
c write(6,'(t2,''scalp,scals'',2f10.3)')scalp,scals
DO 133 J = 1,NR
***************
*** 481,491
IF(NTT.LE.0) GO TO 201
EKGP = ENEW3
EKGS = ENEW4
- c
- c make sure ekgp and ekgs dont get stuck at zero
- c
- IF(EKGP.LT.1.d-4) EKGP = EKINP0
- IF(EKGS.LT.1.d-4) EKGS = EKINS0
C
201 NSTEP = NSTEP+1
T = T + DTX
--- 488,493 -----
IF(NTT.LE.0) GO TO 201
EKGP = ENEW3
EKGS = ENEW4
C
201 NSTEP = NSTEP+1
T = T + DTX
-------------------------------------
Routines affected:
MINMD: Routine FORCE in file .../amber4/src/lib/force.f
Routine RUNMD in file .../amber4/src/minmd/runmd.f