********> bugfix.67
Correction Author: David A. Pearlman
Correction Date: 10/20/94

Programs: Gibbs

Severity: moderate

Problem: 

       When attempting to run a TI simulation (IDIFRG=1), the
       dihedral contribution to the free energy could be incorrectly
       calculated if any torsions represented by a Fourier series
       of more than one term change with lambda. That is, if any
       of the torsions that change with lambda are represented by
       a series of dihedral definitions in the PARM file with 
       PN(I) < 0. (See the PARM manual for more details)
  
       This bug does not affect any FEP (IDIFRG=0) runs.
       This bug will also not occur if none of the torsions
       are parameterized with muti-term Fourier series.

Affects: 

       The internal energy reported (BADH) might be in error.

Cause: 

       Algorithmic errors in the routine that calculates the 
       torsional energy.

Fix:  Make the following changes to routine PEPHI in gibb.f


--------------------------------------------------------------------
*** OLD gibb.f
--- NEW gibb.f
***************
*** 3460,3465 ****
--- 3460,3466 ----
        DIMENSION EEPGN(3),ENBPN(3),XXIJ(3),XXKJ(3),XXKL(3),XXIL(3)
        DIMENSION ESINS(2)
        DIMENSION TAUKEP(1),ALMCC(1),IPERC(1)
+       DIMENSION 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)
***************
*** 3488,3507 ****
    820 CONTINUE
        ESINS(1) = 0.0D0
        ESINS(2) = 0.0D0
  C
- C     ----- LOOP OVER TWO SETS -----
- C
-       DO 600 KPERT = 1,2
-       INDD = (KPERT-1)*NPHI
-       IF(IELPER.GT.0 .OR. INBPER.EQ.2) INDD = 0
-       EDIHP(KPERT) = 0.0
- C
        DO 200 IPHI = 1,NPHI
        I3 = IP(IPHI)
        J3 = JP(IPHI)
        K3 = KP(IPHI)
        L3 = LP(IPHI)
        IC = ICP(IPHI+INDD)
        KDIV = 1
  C
  C     ----- L3.LT.0 FOR IMPROPER TORSION ANGLES
--- 3489,3511 ----
    820 CONTINUE
        ESINS(1) = 0.0D0
        ESINS(2) = 0.0D0
+       EDIHP(1) = 0.0D0
+       EDIHP(2) = 0.0D0
  C
        DO 200 IPHI = 1,NPHI
        I3 = IP(IPHI)
        J3 = JP(IPHI)
        K3 = KP(IPHI)
        L3 = LP(IPHI)
+ C
+ C     ----- LOOP OVER TWO SETS -----
+ C
+       DO 600 KPERT = 1,2
+       INDD = (KPERT-1)*NPHI
+       IF(IELPER.GT.0 .OR. INBPER.EQ.2) INDD = 0
        IC = ICP(IPHI+INDD)
+       NDUP(KPERT) = 1
+       IF (KPERT.EQ.1) ICKEEP(NDUP(KPERT)) = IC
        KDIV = 1
  C
  C     ----- L3.LT.0 FOR IMPROPER TORSION ANGLES
***************
*** 3661,3677 ****
  C negligable, and saves us a lot of extra calculations here.
  C
        IF (IDIFRG.EQ.1 .AND. KPERT.EQ.2) THEN
!          ICOLD = ICP(IPHI)
!          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
!          ESINS(1)  = ESINS(1) + RKU*SARGUM * (AP*DN-DPHI)
! C
!          IF (IPHI.GT.MPHI) ESINS(2) = ESINS(2)+RKU*SARGUM * (AP*DN-DPHI)
        END IF
  C
  C Calculate the coordinate coupling contribution, if necessary
--- 3665,3686 ----
  C negligable, and saves us a lot of extra calculations here.
  C
        IF (IDIFRG.EQ.1 .AND. KPERT.EQ.2) THEN
!          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 is zero for unmatched torsions.
  C
!             ESINAD = RKU*SARGUM * (AP*DN-DPHI)
!          ELSE
!             ESINAD = 0.0D0
!          END IF
!          ESINS(1)  = ESINS(1) + ESINAD
!          IF (IPHI.GT.MPHI) ESINS(2) = ESINS(2)+ESINAD
        END IF
  C
  C Calculate the coordinate coupling contribution, if necessary
***************
*** 3746,3757 ****
  C
    301 IF(PN(IC).GT.0) GO TO 121
        IC = IC+1
        GO TO 120
    121 CONTINUE
  C
  C     ----- COMPUTE NONBONDED INTERACTIONS -----
  C
!       IF(KDIV.EQ.0 .OR. KPERT.EQ.2) GO TO 200
  C
        RIL2 = 0.0
        DO 510 M = 1,3
--- 3755,3769 ----
  C
    301 IF(PN(IC).GT.0) GO TO 121
        IC = IC+1
+       NDUP(KPERT) = NDUP(KPERT) + 1
+       IF (KPERT.EQ.1) ICKEEP(NDUP(KPERT)) = IC
        GO TO 120
    121 CONTINUE
+ 
  C
  C     ----- COMPUTE NONBONDED INTERACTIONS -----
  C
!       IF (KDIV.EQ.0 .OR. KPERT.EQ.2) GO TO 600
  C
        RIL2 = 0.0
        DO 510 M = 1,3
***************
*** 4073,4080 ****
        F(L3+2) = F(L3+2)+YA
        F(L3+3) = F(L3+3)+ZA
  C
-   200 CONTINUE
    600 CONTINUE
  C
  C     ----- CALCULATE THE PERTURBATION ENERGY APPROPRIATELY -----
  C
--- 4085,4092 ----
        F(L3+2) = F(L3+2)+YA
        F(L3+3) = F(L3+3)+ZA
  C
    600 CONTINUE
+   200 CONTINUE
  C
  C     ----- CALCULATE THE PERTURBATION ENERGY APPROPRIATELY -----
  C
--------------------------------------------------------------------

Temporary workarounds: 

     Use FEP. If you wish to use TI, and are encountering this error,
     you must install the fix and recompile.

Routines affected:

     GIBBS file ...amber4/src/gibbs/gibb.f        


--