********> bugfix.51
Correction Author: David A. Pearlman
Correction Date: 10/28/92
Programs: GIBBS
Severity: Moderate
Problem:
If a PMF calculation is carried out with using *torsional*
constraints on topological dihedrals of the molecule (i.e.
dihedral angles which represent four atoms joined in the
molecular topology passed from parm), the contributions to the
free energy due to the constrained torsion itself will not
be included.
This problem affects only PMF calculations using torsions.
PMF simulations using bonds and valence angles are not affected.
This problem does not affect PMF calculations for a torsion which
is not a toplogical torsion.
Affects:
The reported free energy for the PMF will be off by the enthalpic
contribution from the underlying torsion. The remainder of the
simulation will be correct, and this missing contribution can be
accurately estimated post-run as described below in the "workarounds"
section.
Cause:
An IF statement was incorrect, which resulted in certain flags
not being appropriately set.
Fix: Make the following change to routine RSTIN
---To routine RSTIN in file rstin.f -----------------------------------------
*** OLD rstin.f
--- NEW rstin.f
***************
*** 513,519 ****
C IN THE NORMAL ARRAYS. IF ITOR=2, DON'T ZERO OUT PARAM. (WILL BE
C NEEDED FOR TORSIONAL CONTRIBUTION )
C
! IF (IZE.EQ.1.AND.ITOR.NE.2) THEN
IF (NPHIH.GT.0)
* CALL ZERPRM(NPHIH,IX(I40),IX(I42),IX(I44),IX(I46),
* IX(I48),4,0,PK,IAT,IFH,IZE)
--- 513,519 ----
C IN THE NORMAL ARRAYS. IF ITOR=2, DON'T ZERO OUT PARAM. (WILL BE
C NEEDED FOR TORSIONAL CONTRIBUTION )
C
! IF (IZE.EQ.1.OR.ITOR.EQ.2) THEN
IF (NPHIH.GT.0)
* CALL ZERPRM(NPHIH,IX(I40),IX(I42),IX(I44),IX(I46),
* IX(I48),4,0,PK,IAT,IFH,IZE)
------------------------------------------------------------------------------
Temporary workarounds:
You must fix the code and recompile to avoid the problem.
HOWEVER, if you have already run a simulation affected by this
problem, you can obtain the correct potential by adding to the
values reported by Gibbs the potential energies for rotation
about the constrained torsion. A total energy corresponding to
a synchronous rotation about the constrained torsion by
all torsions centered on the central bond of the torsion should
be used. E.g. the correction for a constrained X-CT-CT-X
(3-fold aliphatic rotation with 9 specific torsion centered on the
CT-CT bond) would be
Ecorr = 1.3 * (1 + cos(3*tau)),
where tau is the value of the torsion angle (radians).
The correction will be approximate, to the extent that the 9 specific
torsions do not, in fact, move exactly in synch, and you are constraining
only one of them. The error due to this approximation is probably very
small.
Routines affected:
GIBBS Routine RSTIN in file ...amber4/src/rstin.f
---