********> bugfix.69 Author: David A. Pearlman Date: 03/26/96 Programs: Gibbs Severity: Serious Problem: When running a free energy calculation using Thermodynamic Integration (IDIFRG=1), the reported free energy contribution from torsions (included in the "BADH" component) may be wrong if the periodicity (N) or phase (phi) of any torsion changes with lambda. This problem can only occur when N or phi changes with lambda, and only with Thermodynamic Integration. The problem does NOT occur if only the force constant (K) for a torsion changes, or if FEP (IDIFRG=0) has been used. Affects: The reported total free energy and the contribution from internals (BADH) may be incorrect. All other components will be correct in all cases. (Correctness may generally be test by comparing the results for TI and FEP over a small number of windows run in vacuo). Workarounds: None Fix: Make the following changes to routine PEPHI ---To routine PEPHI in file gibb.f ----------------------------------------- *** OLD gibb.f --- NEW gibb.f *************** *** 4044,4053 **** DIMENSION CONTOR(4,1),ITPT(1) DIMENSION DC(6),T(6),EDIHP(3),EELPT(3),EEPGP(3),EEPGE(3) DIMENSION EEPGN(3),ENBPN(3),XXIJ(3),XXKJ(3),XXKL(3),XXIL(3) ! DIMENSION ESINS(2),DALMF(2),DALMR(2),RMLT(2) DIMENSION TAUKEP(1),ALMCC(1),IPERC(1) DIMENSION DUMA(4),IDUMA(4) ! DIMENSION EKEEP(20),ESINKP(20),ICKEEP(20),NDUP(2) C EQUIVALENCE (XXIJ(1),XIJ),(XXIJ(2),YIJ),(XXIJ(3),ZIJ) EQUIVALENCE (XXKJ(1),XKJ),(XXKJ(2),YKJ),(XXKJ(3),ZKJ) --- 4044,4053 ---- DIMENSION CONTOR(4,1),ITPT(1) DIMENSION DC(6),T(6),EDIHP(3),EELPT(3),EEPGP(3),EEPGE(3) DIMENSION EEPGN(3),ENBPN(3),XXIJ(3),XXKJ(3),XXKL(3),XXIL(3) ! DIMENSION DALMF(2),DALMR(2),RMLT(2) DIMENSION TAUKEP(1),ALMCC(1),IPERC(1) DIMENSION DUMA(4),IDUMA(4) ! DIMENSION EKEEP(20),ICKEEP(20),NDUP(2) C EQUIVALENCE (XXIJ(1),XIJ),(XXIJ(2),YIJ),(XXIJ(3),ZIJ) EQUIVALENCE (XXKJ(1),XKJ),(XXKJ(2),YKJ),(XXKJ(3),ZKJ) *************** *** 4093,4100 **** EEPGN(I) = 0.0D0 ENBPN(I) = 0.0D0 820 CONTINUE - ESINS(1) = 0.0D0 - ESINS(2) = 0.0D0 EDIHP(1) = 0.0D0 EDIHP(2) = 0.0D0 C --- 4093,4098 ---- *************** *** 4261,4277 **** EDIHP(KPERT) = EDIHP(KPERT)+E IF(IPHI.GT.MPHI) EEPGP(KPERT) = EEPGP(KPERT)+E C - C Store k(lam)*sin(n*tau-phi) * {tau*((n(1)-n(0)) - (phase(1)-phase(0))), - C for use in differential calculation of energy. ESINS(1) is - C the sum for all pert torsions. ESINS(2) is the same thing, summed only - C over torsions with 1/more atoms not in pert group. We only need to - C do this evaluation on the second iteration of the KPERT flag. - C - C Note that we do not bother to evaluate K(lam), tau, and phi at ALMH, - C which in the case of slow growth (only), will actually be lamda+-d_lam/2. - C We only evaluate them at lambda itself. The difference is most surely - C negligable, and saves us a lot of extra calculations here. - C C Because we can have a Fourier series of more than one term C defined for a specific torsion (in which case PN(I) < 0 for all C but one of the terms), we need to keep track of the number --- 4259,4264 ---- *************** *** 4280,4317 **** C torsion in the lambda=1 (KPERT=1) state; NDUP(2) is the number C of duplicates for this torsion in the lambda=0 (KPERT=2) state. C ! C For torsions that have a counterpart in both the lambda = 0 and ! C lambda = 1 states, we need to calculate ! C ESINAD = {tau*((n(1)-n(0)) - (phase(1)-phase(0))). ! C If there is no counterpart, then this term effectively drops out, and we ! C set it to 0 here. ! C ! C We store ESINAD in ESINKP(NDUP(2)). We store E in EKEEP(NDUP(KPERT)). ! C We need these to calculate the free energy component contributions, C if requested (below). C ! IF (IDIFRG.EQ.1 .AND. KPERT.EQ.2) THEN ! c ICOLD = ICP(IPHI) ! IF (NDUP(2).LE.NDUP(1)) THEN ! ICOLD = ICKEEP(NDUP(2)) ! DN = ABS(PN(ICOLD)) - ABS(PN(IC)) ! DPHI = PHASE(ICOLD) - PHASE(IC) ! RKU = ALM(1)*PK(ICOLD) + ALM(2)*PK(IC) ! IF (IOLDFM.EQ.0) SARGUM = SIN(ABS(PN(IC))*AP-PHASE(IC)) C ! C Define for/rev values of K (RKU), N (RN) and PHI (RPHI) C ! ESINAD = RKU*SARGUM * (AP*DN-DPHI) ! ELSE ! ESINAD = 0.0D0 ! END IF C ! ESINS(1) = ESINS(1) + ESINAD ! IF (IPHI.GT.MPHI) ESINS(2) = ESINS(2)+ESINAD ! ESINKP(NDUP(2)) = ESINAD ! END IF ! EKEEP(NDUP(KPERT)) = E C IF (ITORON.GT.0) GO TO 301 C C If FEP is being peformed with torsional constraints, --- 4267,4299 ---- C torsion in the lambda=1 (KPERT=1) state; NDUP(2) is the number C of duplicates for this torsion in the lambda=0 (KPERT=2) state. C ! C We store E in EKEEP(NDUP(KPERT)). ! C We need this to calculate the free energy component contributions, C if requested (below). C ! EKEEP(NDUP(KPERT)) = E C ! C --------------------------------------------------------------------------- C ! C Bugfix 3/26/96: Was incorrectly calculating the TI derivative for ! C case where n (periodicity) or phi (phase changing). Because the mixed ! C lambda torsional contribution is calculated as ! C lambda* K1 (1 + cos (n1*tau-phi1)) + (1-lambda) K0 (1+cos (n0*tau-phi0)) ! C the only lambda dependence appears in the pre-force constant terms. ! C Therefore, there's no reason to add additional contributions when ! C n or phi changes... C ! C Previously, was INCORRECTLY calculating the derivative for changes in ! C n or phi as ! C k(lam)*sin(n*tau-phi) * {tau*((n(1)-n(0)) - (phase(1)-phase(0))) C + C A block of code that was calculating values for ESINAD, ESINS() and + C ESINKP() here has been deleted, and references to these variables + C (which are now always 0) have been removed elsewhere in the code. + C + C -- dap (3/26/96) + C --------------------------------------------------------------------------- + C IF (ITORON.GT.0) GO TO 301 C C If FEP is being peformed with torsional constraints, *************** *** 4818,4827 **** IF (IDIFRG.EQ.1) THEN IF (KPERT.EQ.1) THEN ECNTF1 = -EKEEP(IMM) * DLF * FACT1 ! ECNTF1 = -EKEEP(IMM) * DLR * FACT1 ELSE ! ECNTF1 = (+EKEEP(IMM)+ESINKP(IMM)) * DLF * FACT1 ! ECNTF1 = (+EKEEP(IMM)+ESINKP(IMM)) * DLR * FACT1 END IF END IF IF (IAPER(IIAT).GT.0) THEN --- 4800,4809 ---- IF (IDIFRG.EQ.1) THEN IF (KPERT.EQ.1) THEN ECNTF1 = -EKEEP(IMM) * DLF * FACT1 ! ECNTR1 = -EKEEP(IMM) * DLR * FACT1 ELSE ! ECNTF1 = (+EKEEP(IMM)) * DLF * FACT1 ! ECNTR1 = (+EKEEP(IMM)) * DLR * FACT1 END IF END IF IF (IAPER(IIAT).GT.0) THEN *************** *** 4963,4970 **** C DEDLT is for all torsions. DEDLT2 is for torsions not entirely in pert grp. C EP = (EDIHP(1)*ALM(1)+EDIHP(2)*ALM(2)) ! DEDLT = EDIHP(1) - EDIHP(2) - ESINS(1) ! DEDLT2 = EEPGP(1) - EEPGP(2) - ESINS(2) C EPERT(40) = EP EPERTD(41) = DEDLT*DLF --- 4945,4952 ---- C DEDLT is for all torsions. DEDLT2 is for torsions not entirely in pert grp. C EP = (EDIHP(1)*ALM(1)+EDIHP(2)*ALM(2)) ! DEDLT = EDIHP(1) - EDIHP(2) ! DEDLT2 = EEPGP(1) - EEPGP(2) C EPERT(40) = EP EPERTD(41) = DEDLT*DLF -----------------------------------------------------------------------------