********> bugfix.39
Correction Author: David A. Pearlman
Correction Date: 5/09/92  

Programs: SANDER

Severity: moderate

Problem: 

       The "weight change" instructions to modify the force constants for the
       torsional (TORSION) and improper (IMPROP) terms do not function
       correctly. The ALL and INTERN weight change instructions will also
       work improperly as they apply to the torsion angles.

       The standard torsional term is of the form
             K(1) + K(2)*COS(N*TAU-PHI)

       where K(1) = K(2) (the force constant).

       Because of the bug, when an attempt was made to change the force
       constant for torsions, only K(1) was being modified. K(2) remained
       as in the original topology file from PARM.

       This problem does *not* affect any added torsional restraints. The
       force constants for those terms are properly modified by REST
       commands.

Affects: 

       Because only the constant term in the torsional energy expression
       is changing, the torsional derivatives are not affected, so the
       trajectory remains as if the torsional force constants had not 
       been changed. The reported torsional energy changes, but the 
       change only reflects a change in the relative 0-point energy for 
       the torsional terms.

Cause: 

       The nmr-refinement-related code, including that for the "weight change"
       instructions, was originally written to be included in an older
       version of Amber (3.0). In the transition to Amber 3.0A, the
       torsional routine was vectorized, and a couple of new arrays were
       added. These arrays were initialized once, at the start of the program,
       and their values depended on the force constants for the torsional
       terms. When the nmr-refinement code was merged with the updated
       code, this change was inadvertantly overlooked.

       For the fix, we update the appropriate arrays any time the torsional
       force constants are changed.

Fix:  Make the following changes. See 0README for a description of 
      context diffs (used for modwt.f).

---To routine MODWT in src/sander/modwt.f-------------------------------------


*** OLD modwt.f
--- NEW modwt.f	
***************
*** 153,158 ****
--- 153,159 ----
        END IF
        IVDW = 0
        IHBP = 0
+       ITORFC = 0
  C
  C Initialize the change indicator. Any energy term not modified by one
  C of the change instructions should be set to the default weight of 1.0
***************
*** 251,256 ****
--- 252,258 ----
                ICHANG(3) = 1
                IF (ABS(WT).LE.SMALL) ICHOLD(3) = 2
                IF (FIXED) ICHOLD(3) = 2
+               ITORFC = 1
            END IF
  C
  C TYPE = VDW/ATTRACT/NB/ALL:
***************
*** 343,348 ****
--- 345,351 ----
                ICHANG(9) = 1
                IF (ABS(WT).LE.SMALL) ICHOLD(9) = 2
                IF (FIXED) ICHOLD(9) = 2
+               ITORFC = 1
            END IF
  C
  C TYPE = REST/RESTS:
***************
*** 488,492 ****
--- 491,499 ----
  C
        IF (IHBP.EQ.1) CALL DECNVH(AG,BG,NHB,RADHB,WELHB)
  C
+ C SETGMS sets a couple of arrays which depend on the torsional term
+ C force constants, if those force constants were changed here.
+ C
+       IF (ITORFC.EQ.1) CALL SETGMS(NUMPHI+NIMPRP)
        RETURN
        END
------------------------------------------------------------------------------
---Add the following routine to the end of file src/lib/set.f-----------------

      SUBROUTINE SETGMS(NUMPHI)
C
C Subroutine SET GaMS.
C
C This routine sets the values of the GAMC() and GAMS() arrays, which
C are used in vectorized torsional energy routines. This routine is only
C called when the torsional force constants are changed in routine MODWT.
C Otherwise, these arrays are set only once, by a call to DIHPAR from RDPARM,
C at the start of the program. This is an abbreviated version of DIHPAR which
C passes most arguments by common, not call-list.
C
C Author: David A. Pearlman
C Date: 5/92
C
#ifdef DPREC 
      implicit double precision (a-h,o-z)
#endif
C
      COMMON/PARMS/RK(5000),REQ(5000),TK(900),TEQ(900),PK(900),PN(900),
     $             PHASE(900),CN1(1830),CN2(1830),SOLTY(60),
     $             GAMC(900),GAMS(900),IPN(900),FMN(900)
C
      DATA ZERO,ONE,TENM3,TENM10/0.0d+00,1.0d+00,1.0d-03,1.0d-06/
      DATA FOUR,PI/4.0d+00,3.141592653589793d+00/
C
      PIM = FOUR*ATAN(ONE)
      DO 100 I = 1,NUMPHI
        DUM = PHASE(I)
        IF(ABS(DUM-PI).LE.TENM3) DUM = SIGN(PIM,DUM)
        DUMC = COS(DUM)
        DUMS = SIN(DUM)
        IF(ABS(DUMC).LE.TENM10) DUMC = ZERO
        IF(ABS(DUMS).LE.TENM10) DUMS = ZERO
        GAMC(I) = DUMC*PK(I)
        GAMS(I) = DUMS*PK(I)
  100 CONTINUE
      RETURN
      END
------------------------------------------------------------------------------
------------------------------------------------------------------------------


Temporary workarounds: 

     If you wish to modify the force constants for the torsional/improper
     terms, there are no temporary workarounds. You must install this fix.

Routines affected:

     SANDER  Routine MODWT in  .../amber4/src/sander/modwt.f
     SANDER  File .../amber4/src/lib/set.f