********> bugfix.68 Correction Author: David A. Pearlman Correction Date: 01/11/95 Programs: Gibbs Severity: serious Problem: Attempting to calculate a Potential of Mean Force curve using the FEP method in a periodic box could give erroneous results. The manifestation of this bug is a PMF curve that appears extremely noisy and/or non-physical. This bug ONLY affects the free energy contributions of *ADDED* constraints (INTR >0 and ITOR=2 for the added constraint), and only when FEP is used. It does NOT affect the calculation of the "PMF bond contribution" in free energy calculations (NCORC=1), nor does it affect ANY simulation run using Thermodynamic Integration (IDIFRG = 1). This bug is only manifest when 1) two DISTINCT molecules are "joined" through an added constraint and 2) at some point during the simulation the constraint crosses the periodic boundary, i.e. those two molecules are not related by the unitary symmetry operator. Affects: The free energy PMF profile generated can be wrong. Cause: Determining the constraint contribution when using FEP requires that the "lambda+-d_lambda" states be generated by rigid body translation/rotation about constrained internals. The routine that carries out this translation or rotation assumes that the atoms involved in the internal are related by a unitary symmetry operator. If, instead, you have a situation such as ------------------------ | | (A) | A B | | | | | | | ------------------------ which actually corresponds to the unitary symmetry operator case ------------------------ | | (B) B | A | | | | | | | ------------------------ the routine responsible for the rigid body translations/rotations would fail to work properly. Fix: Make the following changes to giba.f, which include adding a new subroutine UNIIMG which takes care of making sure any molecules connected by added constraints are related by the unitary symmetry operator. -------------------------------------------------------------------- *** OLD giba.f --- NEW giba.f *************** *** 432,437 **** --- 432,443 ---- C CALL GTLAM2(IPERON,IPRNOW,ALMCC,IDWIDE,IBNDLM, * X(L30),X(L40),NRP,RDCORC,MOVEOK) + C + CALL UNIIMG( NBONH ,NBONA ,NBPER ,NBRST , + * NARST ,NPRST ,IX(ICR28) ,IX(ICR20) ,NSPM , + * IX(I70) ,X(L95) ,BOX ,NATOM ,IDIFRG , + * X(L30) ,0) + END IF C C Call RUNMD NRUN times. Each time, NSTLIM steps will be performed. *************** *** 4208,4213 **** --- 4214,4228 ---- IF(NTB.EQ.0) GOTO 610 IF (ITRSLU.EQ.1) CALL SHIAG(NSPM,1,0,NSP,X, nr) IF (ITRSLU.NE.1) CALL SHIAG(NSPM,NSPSOL,NSPSTR,NSP,X, nr) + C + C UNIIMG makes sure that any molecules connected by constraints are + C connected by unit symmetry when FEP used. + C + CALL UNIIMG( NBONH ,NBONA ,NBPER ,NBRST , + * NARST ,NPRST ,IX(ICR28) ,IX(ICR20) ,NSPM , + * IX(I70) ,XX(L95) ,BOX ,NATOM ,IDIFRG , + * X ,1) + 610 IF(NTP.LE.0) GOTO 625 EKCMT(1) = 0.0 EKCMT(2) = 0.0 *************** *** 4367,4373 **** --- 4382,4394 ---- 100 IF(NTB .ne. 0) then IF (ITRSLU.EQ.1) CALL SHIAG(NSPM,1,0,NSP,X, nr) IF (ITRSLU.NE.1) CALL SHIAG(NSPM,NSPSOL,NSPSTR,NSP,X, nr) + CALL UNIIMG( NBONH ,NBONA ,NBPER ,NBRST , + * NARST ,NPRST ,IX(ICR28) ,IX(ICR20) ,NSPM , + * IX(I70) ,XX(L95) ,BOX ,NATOM ,IDIFRG , + * X ,1) endif + + C C If we're doing microwindows, see if we've collected enough data C for this window. If so, set WIDONE = .true. *************** *** 7983,7986 **** --- 8004,8242 ---- END IF C RETURN + END + c---------------------------------------------------------------------- + SUBROUTINE UNIIMG( NBONH ,NBONA ,NBPER ,NBRST , + * NARST ,NPRST ,ITORTY ,NATRST ,NSPM , + * NSP ,ITAG ,BOX ,NATOM ,IDIFRG , + * X ,IFLAG) + C + C Subroutine UNIfy IMaGe + C + C This routine makes sure that any molecules joined by added + C constraints are described by a single symmetry operator. + C That is, no periodic imaging is required to calculate + C the true internal parameters between the molecules. This is + C required to properly apply the rigid body translations and + C rotations required if calculating free energy contributions when + C employing FEP (not needed with TI). + C + C Once the solute molecules are translated, the central box will + C be repositioned so that it is centered on the tagged molecules + C (unless IFLAG=1). + C + C Note that this routine assumes that any individual molecule exists + C as a single unified image, and it assumes a rectangular box (NTM=0). + C + C Author: David A. Pearlman + C Date: 1/95 + C + C + C Argument list: + C ------------- + C NBONH: Number of bonds containing at least one hydrogen, and not being + C perturbed, in the original topology list (from PARM). + C NBONA: Number of bonds containing no hydrogens, and not being perturbed, + C in the original topology list (from PARM). + C NBPER: Number of perturbed bonds in the original list (from PARM). + C NBRST: Total number of bond restraints defined on input to this run + C (See routine RSTIN) + C NARST: NBRST + Total number of angle restraints defined on input to this run + C (See routine RSTIN) + C NPRST: NARST + Total number of torsional restraints defined on input to this + C run (See routine RSTIN). + C ITORTY(I): Flag for each internal, defining the type of + C con/res-traint. Only internals for which ITORTY(I)=2 are to be + C constrained using this routine. + C NATRST(4,I): The atoms defining added constraint I. + C NSPM: Number of molecules in the system + C NSP(I): Number of atoms in molecule I + C ITAG(*): Scratch array used to keep track of which residues are involved + C in constraints. (Length >= NRES) + C BOX(3): Dimensions of the box in the three directions + C NATOM: Number of atoms in the system + C IDIFRG: =0, then FEP being used + C IDIFRG: =1, then TI being used. + C X(3,*): Atomic coordinate array + C IFLAG: 0, recenter box on any molecules related by constraints before return + C 1, just ensure all molecules decribed by same symmetry operator, + C do no recenter box. + C + C + C + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + C + DIMENSION ITORTY(*) + DIMENSION NATRST(4,*) + DIMENSION IM(4),BOX(3) + DIMENSION NSP(*),ITAG(*),X(3,*) + C + C Don't need to do anything if NPRST = 0 (no added constraints) or + C if IDIFRG = 1 (TI, not FEP) + C + IF (NPRST.LE.0) RETURN + IF (IDIFRG.EQ.1) RETURN + C + DO 5 I = 1,NSPM + ITAG(I) = 0 + 5 CONTINUE + C + C Loop over all constrained and perturbed bonds, and over + C all added constraints. + C + ITOTBN = NBONH+NBONA+NBPER + ISTRT = ITOTBN+1 + C + DO 10 ICONS=ISTRT,NPRST+ITOTBN + C + C Skip this interation for certain bonds, depending on the value of NTC: + C + IP = ICONS - ITOTBN + C + C For added constraints, skip if it's actually a restraint (ITORTY(IP).NE.2) + C + IF (ITORTY(IP) .NE. 2) GO TO 10 + C + C Determine the type of constraint and set the atom pointers IM + C + IM(1) = NATRST(1,IP)/3 + 1 + IM(2) = NATRST(2,IP)/3 + 1 + IM(3) = 0 + IM(4) = 0 + IF (IP.GT.NARST) THEN + IM(3) = NATRST(3,IP)/3 + 1 + IM(4) = NATRST(4,IP)/3 + 1 + ELSE IF (IP.GT.NBRST) THEN + IM(3) = NATRST(3,IP)/3 + 1 + END IF + C + C Figure out which molecule each atom of the constraint belongs to + C and tag that molecule as being in a constraint + C + INDX = 0 + DO 100 IMOL = 1,NSPM + IF (IM(1).GT.INDX .AND. IM(1).LE.INDX+NSP(IMOL)) THEN + ITAG(IMOL)= 1 + IM(1) = 0 + END IF + IF (IM(2).GT.INDX .AND. IM(2).LE.INDX+NSP(IMOL)) THEN + ITAG(IMOL)= 1 + IM(2) = 0 + END IF + IF (IM(3).GT.INDX .AND. IM(3).LE.INDX+NSP(IMOL)) THEN + ITAG(IMOL)=1 + IM(3) = 0 + END IF + IF (IM(4).GT.INDX .AND. IM(4).LE.INDX+NSP(IMOL)) THEN + ITAG(IMOL)=1 + IM(4) = 0 + END IF + IF (IM(1).EQ.0 .AND. IM(2).EQ.0 .AND. + * IM(3).EQ.0 .AND. IM(4).EQ.0) GO TO 10 + INDX = INDX + NSP(IMOL) + 100 CONTINUE + 10 CONTINUE + C + C Figure out the centers of geometry of all the tagged molecules, and + C translate them all to be in same symmetry frame as the first tagged + C molecule. + C + IHIT = 0 + INDX = 0 + DO 110 IMOL = 1,NSPM + IF (ITAG(IMOL).EQ.1) THEN + IHIT = IHIT + 1 + XCOM = 0.0D0 + YCOM = 0.0D0 + ZCOM = 0.0D0 + DO 120 JJ = INDX+1,INDX+NSP(IMOL) + XCOM = XCOM + X(1,JJ) + YCOM = YCOM + X(2,JJ) + ZCOM = ZCOM + X(3,JJ) + 120 CONTINUE + XCOM = XCOM/FLOAT(NSP(IMOL)) + YCOM = YCOM/FLOAT(NSP(IMOL)) + ZCOM = ZCOM/FLOAT(NSP(IMOL)) + C + C If IHIT=1 (first tagged molecule), use this as the "mother" position. + C Otherwise, translate the current molecule to be as close as possible + C to the mother: + C + IF (IHIT.EQ.1) THEN + XCOM0 = XCOM + YCOM0 = YCOM + ZCOM0 = ZCOM + ELSE + XDIF = XCOM-XCOM0 + YDIF = YCOM-YCOM0 + ZDIF = ZCOM-ZCOM0 + BOXTX = 0.0D0 + BOXTY = 0.0D0 + BOXTZ = 0.0D0 + 121 IF (XDIF+BOXTX.GT.BOX(1)/2.0D0) THEN + BOXTX = -BOX(1) + GO TO 121 + ELSE IF (XDIF+BOXTX.LT.-BOX(1)/2.0D0) THEN + BOXTX = +BOX(1) + GO TO 121 + END IF + 122 IF (YDIF+BOXTY.GT.BOX(2)/2.0D0) THEN + BOXTY = -BOX(2) + GO TO 122 + ELSE IF (YDIF+BOXTY.LT.-BOX(2)/2.0D0) THEN + BOXTY = +BOX(2) + GO TO 122 + END IF + 123 IF (ZDIF+BOXTZ.GT.BOX(3)/2.0D0) THEN + BOXTZ = -BOX(3) + GO TO 123 + ELSE IF (ZDIF+BOXTZ.LT.-BOX(3)/2.0D0) THEN + BOXTZ = +BOX(3) + GO TO 123 + END IF + C + C Translate the second molecule so both com's are in the same box + C + DO 130 JJ = INDX+1,INDX+NSP(IMOL) + X(1,JJ) = X(1,JJ) + BOXTX + X(2,JJ) = X(2,JJ) + BOXTY + X(3,JJ) = X(3,JJ) + BOXTZ + 130 CONTINUE + END IF + INDX = INDX + NSP(IMOL) + END IF + 110 CONTINUE + C + C If any molecules were translated, now recenter the box on all tagged + C molecules (only if IFLAG .NE. 1): + C + IF (IHIT.GT.1 .AND. IFLAG.NE.1) THEN + XCOM = 0.0D0 + YCOM = 0.0D0 + ZCOM = 0.0D0 + NUSE = 0 + INDX = 0 + DO 140 IMOL = 1,NSPM + IF (ITAG(IMOL).GT.0) THEN + DO 150 IA = INDX+1,INDX+NSP(IMOL) + XCOM = XCOM + X(1,IA) + YCOM = YCOM + X(2,IA) + ZCOM = ZCOM + X(3,IA) + 150 CONTINUE + NUSE = NUSE + NSP(IMOL) + END IF + INDX = INDX + NSP(IMOL) + 140 CONTINUE + XTRAN = XCOM/FLOAT(NUSE) - BOX(1)/2.0D0 + YTRAN = YCOM/FLOAT(NUSE) - BOX(2)/2.0D0 + ZTRAN = ZCOM/FLOAT(NUSE) - BOX(3)/2.0D0 + DO 160 I = 1,NATOM + X(1,I) = X(1,I) - XTRAN + X(2,I) = X(2,I) - YTRAN + X(3,I) = X(3,I) - ZTRAN + 160 CONTINUE + END IF + C + RETURN + C END ------------------------------------------------------------------------------ Temporary workarounds: The routine that wraps molecules back into the central box in Gibbs does so on a molecule basis. You can ensure that all the molecules joined by constraints in your system are wrapped back into the central box as one unit by defining these molecules as joined by non-covalent linkages in the LINK program (using a '***' spacer residue). See the LINK documentation for more details. If you are attempting to use this workaround without installing the above fix, ensure that the molecules of interest are properly related by the unitary symmetry operator at the beginning of the MIN/MD/GIBBS run (if they are OK at the beginning of MINMD, they should remain so at the end of MINMD, and at the end of GIBBS). Routines affected: GIBBS file ...amber4/src/gibbs/giba.f --