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