Questions and problems?

Gibbs (free energy) questions

Answers on this page are by David Pearlman

Gibbs setup issues

Leap-generated input

Defining states

My biggest concern is to the input of the two states that will be compared by the program. The literature states that the 0 state is the coordinates given by PARM, and the 1 state is contained within the PREP file. Yet the description of the command line for GIBBS has no input PREP file.

The word 'state' refers to the _interpretation_ given to the coordinates, i.e. the nature of the atoms, bonds, etc. There is only one set of coordinates, which is normally a restrt file from sander or gibbs. The PREP 'state' is carried via LINK and EDIT to the PARM topology file (in leap, it goes directly into the topology, and the perturbed state is defined via the xleap 'unit Editor'). From the gibbs manual:

       The  state  defined  in PREP is the =1 state and the state
  given in the PARM is the =0  state.   The  default  state  from
  which  to  start  the  perturbation  is usually =1, because the
  coordinates which are carried from EDIT to PARM to MIN to GIBBS
  will correspond to the PREP state.  

Can I perturb by having two copies of a residue that don't "see" each other, disappearing one while appearing the other?

Yes; it's complicated though..

Shrinking distances w/ dummy atoms

I am perturbing an ethyl group to nothing. The dummy atoms are separated by a min distance of 0.001. The hydrogens are too close to the other hydrogens and there are bad contacts. However I have at frcmod.dat the nonbonded parameters very low for the atoms ...
        	     R*    epsilon
        	DC  0.001   0.0
        	DH  0.001   0.0
so even though they are very close are much smaller in size but they still create bad contacts.

I set the dummy atoms closer than the corresponding real atoms so when they disappear there will not be a big "cap" in the system.

The problem here is that MD trajectories almost always become very unstable if you try to shrink the distances to dummy atoms to extremely small distances. This has to do with increasing vibrational frequencies of short length bonds. In my experience, 0.2 A is typically around the shortest distance you can shrink a bond to without encountering problems.

If you simply insist on shrinking to a shorter distance, you will need to reduce the timestep--possibly to something considerably smaller than 0.5 fsec...Given the tremendously decreased trajectory efficiency that results (i.e. you have to run a lot long to perform the same number of psec of dynamics), you are usually better off with the (possibly) somewhat diminished sampling efficiency associated with shrinking to a larger bond length like 0.2 A.

Interpreting results

Sign of energy; forward/back

> Given the following:
>  (1) I have defined molecule X in the PREP input.
>  (2) I have defined molecule Y in the PARM input.
>  (3) I have equilibrated two separate systems based on this same
>      set of input files.  One was equilibrated at Lambda=1 (X) and the 
>      other was equilibrated at Lamdba=0 (Y).
> If I then perform a Fixed-width window Free Energy Perturbation on each
> equilibrated system separately (one X[lam=1] -> Y[lam=0] perturbation, 
> one Y[lam=0] -> X[lam=1] ), I get "forward" free energy differences which
> are *both* negative.
> Is it true (as it seems to be based on the info on pg 107 of the manual) 
> that the *sign* of both energies will correspond to that of the Y -> X 
> (lam=0 -> lam=1) perturbation?  In other words, is it true that *real* 
> value of delta G for the X -> Y (1->0) is actually the *negative* of the 
> reported "forward" value found in the output?
Yes. You are correct. The "FORWARD" free energy reported is ALWAYS the value corresponding to the change (0->1), regardless of the actual direction in which the simulation was performed. So if you want the free energy for the (1->0) change, take the negative of the energy reported. (NOTE that this applies to Amber 4.0 and later versions. In earlier (3A, 3.0) versions of Gibbs, the meaning of "Forward" for FEP depended on the direction in which the simulation was run.)

Single vs. double-wide trajectories

I decided to try to save computer time by setting the IDWIDE parameter to 1 in order to turn off the double wide sampling ... I ran a test gas phase perturbation. ... the reported final energy differences were not the same for the two calculations. ... Why should these values differ?

Setting IDWIDE = 1 changes the trajectory, and hence the computed free energy (although values calculated from any trajectory will be the same in the limit of infinite sampling, in the finite sampling domain, they will differ). To understand, when IDWIDE=0 (double wide sampling), the trajectory is calculated at the following points:

        0  dl   2dl   3dl   4dl   ...   1


  dG (forward) = dG (0->dl) + dG (dl->2dl) + dG (2dl->3dl) + ...

  dG (reverse) = dG (dl->0) + dG (2dl->dl) + dG (3dl->2dl) + ...

whereas if IDWIDE = 1, then the trajectory is calculated at the following points:
           dl        3dl        5dl      ...


  dG (forward) = 
            -dG (dl->0) + dG (dl->2dl) - dG (3dl->2dl) + dG (3dl->4dl) + ...

  (You don't get a "reverse" free energy).

Note that since are running MD at a different set of "lambda" points, the trajectory will be different, and the net free energy will be different. The difference will be small, of course, if the free energy is converged.

Dihedral angle warnings

	What might be the reason for the warning that the dihedral 
	angles are not between 0 or 180?
          nb-update: total= 186068; regular= 96590; h-b= 77272; pert= 12206
          Warning (PEPHI): phase angle not 0 or 180 for torsion: 
            0- 3- 12- 15 N= 2.0000 Phi= -1.9722
          Warning (PEPHI): phase angle not 0 or 180 for torsion: 
            0- 3- 12- 42 N= 2.0000 Phi= 1.1519

This message is telling you that some of your torsional parameters have
been defined with a phase angle other than the standard values of 0 or 180.
That's phi in the equation:

    Ephi = K (1 + cos (N tau - phi)). 

The only reason you'd typically use a phi other than 0 or 180 is to
impose a particular torsional restraint. Maybe you're doing this. 
(However, I'd suggest that if you are, that you shouldn't really be
doing it in the parm file, but rather with added restraints in the
Sander input). Anyway...

The integers refer to the atom numbers of the affected torsion. However,
you need to convert them to the actual atom numbers using the formula

    atom-number = (I/3) + 1

By figuring out what the atom type parameters are for this torsion, you
should be able to quickly locate the appropriate torisonal parameter in the
parameters file. Be sure to observe that the dihedral parameter input
is formatted, with a format:



(as described in the PARM manual).

If you didn't intend for the phase angles given above, check to see that
your input falls in the correct columns...

Unexpected delta G = 0

I am trying to calculate the free energy difference of CHCl3 mutated to CCl4 inside some cryptates. I just tried the mutation of CHCl3 ALONE to CCl4 in vacuo. Why is the free enrgy I get identically zero?

In CHCl3 (or ChCl4) you have no non-bonded interactions (remember no non-bonds for 1-3 pairs), so there will be no non-bonded contribution. There are no dihedrals. So the only contributions will be from bond and angle changes. To get these calculated correctly and reported, you need to set INTPRT=5 (report intra-pert-group contributions), NTC=NTF=3 (turn SHAKE on for the solute) and ICORC=1 (calculate the "PMF Bond contribution").

Hysteresis & center of mass motion

In perturbing from ACE-GLY-NME dipeptide to ACE-ALA-NME in gas phase, the free energy values show a large hysteresis between the forward and backward runs (double-wide sampling). On the other hand, if I start with the alanine dipeptide as the initial residue and perturb to the glycine dipeptide, there is little hysteresis.

This large hysteresis for the former run disappears if I simply do *not* remove the translational and rotational motion from the system. This problem does not seem to arise if the perturbations are done in a solvated environment, or in a full protein environment (i.e., GLY and ALA are part of a folded protein).

As I have noted in a few papers, one should never remove COM motion when performing free energy calculations. This will definitely have adverse affects on the calculated free energies. I belive this arises because the moment of inertia of the molecule is changing, which will affect the partitioning among particular translational/rotational/vibrational modes. Removing COM motion messes this up.

PMF: the CORC term

> I'm doing a PMF calculation moving an atom in between two dummy
> atoms in the way that the distance atom-dummy1 increases while
> the distance dummy2 decreases. In order to do so I have to use
> the option NCORC=1. I now would like to understand the meaning of
> the CORC value printed in the listing file and its relation to
> the free energy value.
The "CORC" value is, in principal, the contribution from the constrained bonds to the total free energy. With FEP, the components are only approximate in the best of cases, and for various reasons, a lot of non-bonded interactions get lumped into "CORC" in the output. If you're running FEP, you should definitely ignore the CORC value and focus on the total free energy.

With TI, the CORC contribution is better deconvoluted from the rest. But once again, it is not really a good idea to focus on this component alone for most purposes. I'd suggest you still focus primarily on the total free energy.

"Coordinate coupling" = PMF

How can we carry out coordinate coupling with the program GIBBS on AMBER 41 ? This coupling was mentioned in a publication : Rao, B. and al. JACS, vol 111, 9, 1989, 2135.

The procedure that is termed "Coordinate Coupling" in that paper corresponds to the "PMF Bond Contribution," which was detailed in the paper: Pearlman & Kollman (1991) JCP 94, 4532-4545.

The analogue of this contribution, but for TI calculations, is described in Pearlman (1994) JPC 98, 1487-1493.

To include this contribution in your calculations, you set NCORC=1 and NTC=NTF=3 in the Gibbs input. The program takes care of the rest for you automatically.