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..
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.
> 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.)
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
and
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 ...
and
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.
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:
IPT , JPT , KPT , LPT , IDIVF , PK , PHASE , PN
FORMAT(A2,1X,A2,1X,A2,1X,A2,I4,3F15.2)
(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...
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").
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.
> 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.
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.