********> bugfix.80 - revised Correction Author: David A. Pearlman Correction Date: 12/09/96, revised 2/3/97 Programs: GIBBS Severity: Modest Problem: If one is performing a Free Energy Perturbation (FEP) simulation where all of the following conditions are met, the reported free energies will be incorrect: 1) FEP is being performed, not TI (ITIMTH = 0) 2) Intra-pert group contributions from internals are not determined (INTPRT < 3) 3) "Bond PMF" contributions are being calculated (NCORC = 1) This bug only occurs when this full set of options is chosen. Affects: The reported internal contribution to the free energy (BADH) will likely appear spuriously large. Cause: The routine that determines the "Bond PMF" contributions assumed that intra-perturbed group torsional energies are being calculated, and thus relied on a quantity that should have been previously calculated. However, the routine that should have stored this quantity only stored the appropriate value when INTPRT.GE.3. REVISION NOTE: For those who applied the original bugfix.80, the additional elements of the revision are: all but the last 'hunk' of the connrg.f changes, and the gibb.f change. Fix: Make the following changes to connrg.f ------------------------------------------------------------------------------ *** OLD connrg.f --- NEW connrg.f *************** *** 27,33 **** C************************************************************************ C SUBROUTINE CCNRG(XREF0,XCC,WINV,NBRST,NARST,NPRST,NTHETA, ! * NBONH,NBONA,NBPER,NGPER,MGPER,NPHIH,NPHIA, * NDPER,MDPER,ICBADD,NTB,NPSCAL,NTCEFF,NTCEF2,NTF, * NATOM,NTYPES,INBPER,IELPER,IRINGC,ITIMTH,IBIGM, * IMGSLT,ICORC,INTPRT,NTORCC,ITCUMX,IPERON,IFNRGM, --- 27,33 ---- C************************************************************************ C SUBROUTINE CCNRG(XREF0,XCC,WINV,NBRST,NARST,NPRST,NTHETA, ! * NBONH,NBONA,NBPER,MBPER,NGPER,MGPER,NPHIH,NPHIA, * NDPER,MDPER,ICBADD,NTB,NPSCAL,NTCEFF,NTCEF2,NTF, * NATOM,NTYPES,INBPER,IELPER,IRINGC,ITIMTH,IBIGM, * IMGSLT,ICORC,INTPRT,NTORCC,ITCUMX,IPERON,IFNRGM, *************** *** 154,159 **** --- 154,161 ---- IF (NITP.EQ.0) GO TO 500 EPRT2(21) = 0.0D0 EPRT2(22) = 0.0D0 + EPRT2(24) = 0.0D0 + EPRT2(25) = 0.0D0 EINTRA(1) = 0.0D0 EINTRA(2) = 0.0D0 CALL ZEROUT(3*NATOM,XX(L95)) *************** *** 160,170 **** IF (NTF.LE.2) * CALL PBOND(NBPER,IX(I18+NBONA),IX(I20+NBONA), * IX(I22+NBONA),RK,REQ,NTB,XCC,XX(L95),DUMM1, ! * EPRT2,0,ALM,IELPER,INBPER,DUMM2,IX(I68),IFIN,0, * 0) CALL PANGL(NGPER,IX(I32+NTHETA),IX(I34+NTHETA), * IX(I36+NTHETA),IX(I38+NTHETA),TK,TEQ,NTB,XCC, ! * XX(L95),DUMM1,EPRT2,0,ALM,IELPER,INBPER,DUMM2, * IX(I68),IFIN,0,0) CALL PEPHI(NDPER,NATOM,IX(I50+NPHIA),IX(I52+NPHIA), * IX(I54+NPHIA),IX(I56+NPHIA),IX(I58+NPHIA), --- 162,172 ---- IF (NTF.LE.2) * CALL PBOND(NBPER,IX(I18+NBONA),IX(I20+NBONA), * IX(I22+NBONA),RK,REQ,NTB,XCC,XX(L95),DUMM1, ! * EPRT2,MBPER,ALM,IELPER,INBPER,DUMM2,IX(I68),IFIN,0, * 0) CALL PANGL(NGPER,IX(I32+NTHETA),IX(I34+NTHETA), * IX(I36+NTHETA),IX(I38+NTHETA),TK,TEQ,NTB,XCC, ! * XX(L95),DUMM1,EPRT2,MGPER,ALM,IELPER,INBPER,DUMM2, * IX(I68),IFIN,0,0) CALL PEPHI(NDPER,NATOM,IX(I50+NPHIA),IX(I52+NPHIA), * IX(I54+NPHIA),IX(I56+NPHIA),IX(I58+NPHIA), *************** *** 178,188 **** C for the internals. Place these totals in EPSTG(26)/EPSTG(27) C IF (IFR.EQ.1) THEN ! EPSTG(10) = EPRT2(21) + EPRT2(31) + EPRT2(41) EPSTG(20) = EPSTG(10)-EPSTG(5) EPSTG(26) = EPSTG(20) ELSE IF (IFR.EQ.2) THEN ! EPSTG(15) = EPRT2(22) + EPRT2(32) + EPRT2(42) EPSTG(25) = EPSTG(15)-EPSTG(5) EPSTG(27) = EPSTG(25) END IF --- 180,198 ---- C for the internals. Place these totals in EPSTG(26)/EPSTG(27) C IF (IFR.EQ.1) THEN ! IF (INTPRT.GE.3) THEN ! EPSTG(10) = EPRT2(21) + EPRT2(31) + EPRT2(41) ! ELSE ! EPSTG(10) = EPRT2(24) + EPRT2(34) + EPRT2(52) ! END IF EPSTG(20) = EPSTG(10)-EPSTG(5) EPSTG(26) = EPSTG(20) ELSE IF (IFR.EQ.2) THEN ! IF (INTPRT.GE.3) THEN ! EPSTG(15) = EPRT2(22) + EPRT2(32) + EPRT2(42) ! ELSE ! EPSTG(15) = EPRT2(25) + EPRT2(35) + EPRT2(55) ! END IF EPSTG(25) = EPSTG(15)-EPSTG(5) EPSTG(27) = EPSTG(25) END IF ------------------------------------------------------------------------------ Fix: Make the following changes to gibb.f ------------------------------------------------------------------------------ *** OLD gibb.f --- NEW gibb.f *************** *** 2369,2375 **** CALL SECOND(TIME0) IF (IDIFRG.EQ.0) THEN CALL CCNRG(XX(L26),XX(L25),XX(L20),NBRST,NARST,NPRST,NTHETA, ! * NBONH,NBONA,NBPER,NGPER,MGPER,NPHIH,NPHIA, * NDPER,MDPER,ICBADD,NTB,NPSCAL,NTCEFF,NTCEF2, * NTF,NATOM,NTYPES,INBPER,IELPER,IRINGC,ITIMTH, * IBIGM,IMGSLT,ICORC,INTPRT,NTORCC,ITCUMX,IPERON, --- 2369,2375 ---- CALL SECOND(TIME0) IF (IDIFRG.EQ.0) THEN CALL CCNRG(XX(L26),XX(L25),XX(L20),NBRST,NARST,NPRST,NTHETA, ! * NBONH,NBONA,NBPER,MBPER,NGPER,MGPER,NPHIH,NPHIA, * NDPER,MDPER,ICBADD,NTB,NPSCAL,NTCEFF,NTCEF2, * NTF,NATOM,NTYPES,INBPER,IELPER,IRINGC,ITIMTH, * IBIGM,IMGSLT,ICORC,INTPRT,NTORCC,ITCUMX,IPERON, ------------------------------------------------------------------------------ Temporary workarounds: You can avoid this problem by simply setting INTPRT.GE.3. In general, it is *strongly* recommended that for all simulations one determine *all* intra-perturbed group contributions by setting INTPRT=5 (see, e.g. Pearlman & Connelly (1995) JMB 248, 696). These contributions can frequently be substantial, and are not particularly problematic or costly to determine. Routines affected: GIBBS Routine CCNRG in file amber41/src/gibbs/connrg.f GIBBS Routine FORCE in file .../amber4/src/gibbs/gibb.f --