********> 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