********> bugfix.65 Correction Author: David A. Pearlman Correction Date: 08/19/94 Programs: Gibbs Severity: moderate Problem: A Gibbs program crash might occur if one attempts to perform a change where ALL the following conditions hold: 1) An atom, "A1", is perturbed that participates in "hydrogen bonding" interactions at either/both of the endpoints (lambda=0/ lambda=1). An interaction is described as "hydrogen bonding" only if explicit 10-12 parameters are provided for that interaction in the parameter file. For example, HO OS 7557. 2385. defines interactions between HO and OS atom types to be hydrogen bonding. A "hydrogen bonding" interaction is modeled by a 10-12 potential, rather than a 6-12. 2) The atom "A1" in question is defined to have a 0.0 interaction (vdW radius/well depth) at one/both of the endpoints. For a hydrogen-bonding endpoint, this would mean the 10-12 coefficients are defined as 0 for interactions with another atom of the system, e.g.: A1 A2 0. 0. For a non-hydrogen bonding endpoint, this would require two conditions to hold: A) one/both of the vdW parameters (r*/epsilon) are defined to be 0 in the non-bonded parameters list for atom A1 B) one/both of the vdW parameters (r*/epsilon) are defined to be 0 in the non-bonded parameters list for atom A2 STDA RE A1 0.00 0.000 A2 0.00 0.000 3) The Thermodynamic Integration (TI) method is being used. Affects: Program would terminate unexpectedly (crash) in routine NNBOND. Cause: Division by the effective interaction distance for an atom pair was being performed to calculate a derivative required for TI. This distance could be 0.0, depending on the criteria described above. The division is not actually necessary, and the equations can be redefined without it. Fix: Make the following changes to routine NNBOND in gibbs ---Sander: Module machinedep.f ----------------------------- *** OLD machinedep.f --- NEW machinedep.f *************** *** 1602,1615 **** EPNBDF = EPNB1 - EPNB0 C RD2 = RLAM*RLAM - RD2P = RLAMP*RLAMP - RD2M = RLAMM*RLAMM RD6 = RD2*RD2*RD2 RD10 = RD6*RD2*RD2 ! RD6P = RD2P*RD2P*RD2P ! RD10P = RD6P*RD2P*RD2P ! RD6M = RD2M*RD2M*RD2M ! RD10M = RD6M*RD2M*RD2M C R10WJ = R6W(JN)*RW(JN)*RW(JN) C --- 1606,1627 ---- EPNBDF = EPNB1 - EPNB0 C RD2 = RLAM*RLAM RD6 = RD2*RD2*RD2 RD10 = RD6*RD2*RD2 ! ! RD2P = RLAMP*RLAMP ! RD3P = RD2P*RLAMP ! RD5P = RD2P*RD3P ! RD6P = RD3P*RD3P ! RD9P = RD6P*RD3P ! RD10P = RD5P*RD5P ! ! RD2M = RLAMM*RLAMM ! RD3M = RD2M*RLAMM ! RD5M = RD2M*RD3M ! RD6M = RD3M*RD3M ! RD9M = RD6M*RD3M ! RD10M = RD5M*RD5M C R10WJ = R6W(JN)*RW(JN)*RW(JN) C *************** *** 1626,1636 **** EHB2 = (6.0D0*EPHB*RD10) * R10WJ EHB1P = EPHBDF*RD10P* ( (5.0D0*RD2P*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2P = (60.0D0*EPHBP*RD10P*RLAMDF/RLAMP) * * ( (RD2P*R12W(JN)) - R10WJ ) EHB1M = EPHBDF*RD10M* ( (5.0D0*RD2M*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2M = (60.0D0*EPHBM*RD10M*RLAMDF/RLAMM) * * ( (RD2M*R12W(JN)) - R10WJ ) C ENB1 = (EPNB*RD6*RD6) * R12W(JN) --- 1638,1648 ---- EHB2 = (6.0D0*EPHB*RD10) * R10WJ EHB1P = EPHBDF*RD10P* ( (5.0D0*RD2P*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2P = (60.0D0*EPHBP*RD9P*RLAMDF) * * ( (RD2P*R12W(JN)) - R10WJ ) EHB1M = EPHBDF*RD10M* ( (5.0D0*RD2M*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2M = (60.0D0*EPHBM*RD9M*RLAMDF) * * ( (RD2M*R12W(JN)) - R10WJ ) C ENB1 = (EPNB*RD6*RD6) * R12W(JN) *************** *** 1637,1647 **** ENB2 = (2.0D0*EPNB*RD6) * R6W(JN) ENB1P = EPNBDF*RD6P* ( (RD6P*R12W(JN)) - * (2.0D0 * R6W(JN)) ) ! ENB2P = (12.0D0*EPNBP*RD6P*RLAMDF/RLAMP) * * ( (RD6P*R12W(JN)) - R6W(JN) ) ENB1M = EPNBDF*RD6M* ( (RD6M*R12W(JN)) - * (2.0D0 * R6W(JN)) ) ! ENB2M = (12.0D0*EPNBM*RD6M*RLAMDF/RLAMM) * * ( (RD6M*R12W(JN)) - R6W(JN) ) C E10 = (EHB1 - EHB2) + (ENB1 - ENB2) --- 1649,1659 ---- ENB2 = (2.0D0*EPNB*RD6) * R6W(JN) ENB1P = EPNBDF*RD6P* ( (RD6P*R12W(JN)) - * (2.0D0 * R6W(JN)) ) ! ENB2P = (12.0D0*EPNBP*RD5P*RLAMDF) * * ( (RD6P*R12W(JN)) - R6W(JN) ) ENB1M = EPNBDF*RD6M* ( (RD6M*R12W(JN)) - * (2.0D0 * R6W(JN)) ) ! ENB2M = (12.0D0*EPNBM*RD5M*RLAMDF) * * ( (RD6M*R12W(JN)) - R6W(JN) ) C E10 = (EHB1 - EHB2) + (ENB1 - ENB2) *************** *** 1672,1682 **** EPHBP = ALMH(1)*EPHB1 + ALMH(2)*EPHB0 EPHBM = ALMH(3)*EPHB1 + ALMH(4)*EPHB0 C ! RD2P = RLAMP*RLAMP ! RD2M = RLAMM*RLAMM ! RD10 = RLAM**10 ! RD10P = RLAMP**10 ! RD10M = RLAMM**10 C R10WJ = R6W(JN)*RW(JN)*RW(JN) C --- 1684,1696 ---- EPHBP = ALMH(1)*EPHB1 + ALMH(2)*EPHB0 EPHBM = ALMH(3)*EPHB1 + ALMH(4)*EPHB0 C ! RD2P = RLAMP*RLAMP ! RD2M = RLAMM*RLAMM ! RD10 = RLAM**10 ! RD9P = RLAMP**9 ! RD9M = RLAMM**9 ! RD10P = RD9P*RLAMP ! RD10M = RD9M*RLAMM C R10WJ = R6W(JN)*RW(JN)*RW(JN) C *************** *** 1684,1696 **** EHB2 = (6.0D0*EPHB*RD10) * R10WJ EHB1P = EPHBDF*RD10P* ( (5.0D0*RD2P*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2P = (60.0D0*EPHBP*RD10P*RLAMDF/RLAMP) * * ( (RD2P*R12W(JN)) - R10WJ ) EHB1M = EPHBDF*RD10M* ( (5.0D0*RD2M*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2M = (60.0D0*EPHBM*RD10M*RLAMDF/RLAMM) * * ( (RD2M*R12W(JN)) - R10WJ ) - C E10 = EHB1 - EHB2 E11 = EHB1P + EHB2P E12 = EHB1M + EHB2M --- 1698,1709 ---- EHB2 = (6.0D0*EPHB*RD10) * R10WJ EHB1P = EPHBDF*RD10P* ( (5.0D0*RD2P*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2P = (60.0D0*EPHBP*RD9P*RLAMDF) * * ( (RD2P*R12W(JN)) - R10WJ ) EHB1M = EPHBDF*RD10M* ( (5.0D0*RD2M*R12W(JN)) - * (6.0D0 * R10WJ) ) ! EHB2M = (60.0D0*EPHBM*RD9M*RLAMDF) * * ( (RD2M*R12W(JN)) - R10WJ ) E10 = EHB1 - EHB2 E11 = EHB1P + EHB2P E12 = EHB1M + EHB2M *************** ------------------------------------------------------------------------------ Temporary workarounds: In some cases, it may be possible to remove a hydrogen-bond interaction definition that assign 0.0 coefficients. Routines affected: GIBBS file ...amber4/src/gibbs/machinedep.f --