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