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

---