Free-energy calculations with Gibbs in Amber

1. Start from an equilibrated protein structure
2. Rerun tleap and set the peturbation data
   This is done with the commands
   set x.1.ATOM pertCharge newcharge  (note the capital C - necessary) change the charge of this atom
   set x.1.ATOM pertType newname  (note the capital T  - necessary) change the type
   set x.1.ATOM pertName newname  (note the capital N - necessary) change the name
3. Check that you have got the right data with
   desc x.1.ATOM
4. check x
5. saveamberparmpert x prmtop prmcrd
6. quit
7. set up the perturbation run using the template files below
8. Run gibbs (cgibbs give an input line) (some perturbations can also be run with sander).

gibbs -i gibbs.in -o gibbs.out -r prest -ms pstat -p prmtop -a comp -c input_coord
-i is the name of the input file
-o is the name of the output file
-r is the name of the restart file; it will contain the final coordinates
-ms is the name of a file with a free energy statistics
-p is the name of the topology file (normally prmtop, but pparm as default)
-a is the name of a file with free energy contributions (if iperat is set)
-c is the name of the file with input coordinates. The first time it should be prmcrd; if you have already run gibbs, you should instead use the restart file (prest) from the previous calculation.


Note that gibbs only allows the old prmtop file format:
Convert the prmtop file to old format:
new2oldparm <prmtop >prmtop.old

• Free energy perturbation formula:

Delta G = -RT ln < exp(-(V(lambda_i+1)-V(lambda_i))/RT) >_lambda_i
 
• Thermodynamic integration formula:

DeltaG = Integral_from_0_to_1 <dV/dlambda>_lambda dlambda
 
• State (lambda =) 1 is the one defined in prep (link and edit, i.e. in the unperturbed prmtop file), whereas state (lambda =) 0 is the one defined as a perturbation in parm.

 
• The forward energy always corresponds to the process represented by increasing lambda from 0 to 1. Thus, it is normally the negative of the process run in reality (since state 0 is the perturbation).

 
• Enthalpies and entropies are typically more than ten times less accurate than free energies.

GRP1
ATOM   49 1034
ATOM 1271 1535
...
ATOM 5950 9999
END
END

RES   -1 308
END
END
 4866  4864  4863  4845     0     1     2    1.00000   0.00000
   0.00000   0.00000   0.00000 999.00000    0    0
 4822  4842  4844  4845     0     1     2    1.00000   0.00000
   0.00000   0.00000   0.00000 999.00000    0    0

• isldyn (p. 209) determines the method and the direction of the perturbation. +-2 gives the dynamically modified window method. -2 gives a simulation from lambda = 1 to 0; +2 in the reverse direction.
• nstmeq and nstmul (p. 210) gives the number of steps for equilibration and sampling
• amxmov (p. 211) is the target energy for each lambda step (in kcal/mole)
• amxrst (p. 212) is the energy above which the calculation is discarded and a smaller step is taken. Should typically be 5*amxmov (default)
• almdl0 (p. 212) is the starting lambda step
• dlmax (p. 212) is the maximum lambda step
• idifrg (p. 209) gives the method of the calculation. 1 gives thermodynamic integration
• isande=1 (p. 205) says that enthalpy and entropy should be calculated.
• almda (p. 208) is the stating value of lambda
• idsx0 (p. 202) controls the behaviour of vanishing atoms
• ielper (p. 201) is an option to decouple electrostatics and van der Waals interaction
• ncorc (p. 198) determines if the energy due to constraints should be accumulated. Thus ncorc=1 necessarily for PMF calculations.
• intprt (p. 203) determines if the energy of the perturbed group should be included
• ndmpmc (p. 210) determine the frequency of sampling of perturbation energies
• itimth (p. 199) is the method to calculate constraint free energy contributions with NCOR=1. itimth=0 gives the potential forces method that is best. itimth=1 gives the constraint forces which must be used if there are any constraints forming a ring.
• iperat (p. 205) determine if contributions from various groups should be sampled. It is connected to lines -5 to -7 in the file.
• intr (p. 207) is a flag that there is a PMF calculation. It is connected to the last four lines in the file.

These lines gives the constraint:
atom1 atom2 atom3 atom4 umbrella overwrite_options_in_force_field restraint/constraint lambda_start lambda_end
keq_start eq_start keq_end eq_end options_for_1-and 2-dimensional scans

• irest=1 (p. 190) indicate a restart, where accumulated energies and lambda are read from the prest file.
• ntx=5 (p. 192) indicates that both coordinates and velocities should be read from the prest file. ntx=1 only coordinates.
• init (p. 197) is the way to obtain velocities. init=4 should be combined with ntx=5, whereas init=3 should be used with ntx=1
• nrun (p. 193) is the number of MD runs. The restart file is written only at the end of each run, so it may be vise to run several runs.
• nstlim (p. 197) is the number of MD steps per nrun
• dt (p. 197) is the time step in ps. Use 0.001-0.002 with SHAKE and 0.0005-0.001 without.
• temp0 (p. 194) is the desired temperature.
• ntt=1 (p. 194) is the default Berendsen temperature coupling algorithm.
• taup is the temperature relaxation time. Should be 0.1-0.4, default 0.1.
• ntc (p. 198) is the SHAKE flag. ntc=1 no shake. ntc=3 use shake.
• ntf (p. 200) determines which forces to calculate. Should be set to the same value as ntc.
• ishkfl (p. 199) is the behaviour if shake fails. ishkfl=1 is OK.
• nsnb (p. 201) is the number of steps between each non-bonded list update. Should be about 25 fs. Thus, use nsnb=50 without shake and nsnb=15 with
• cut (p. 204) is the cut off for non-bonded interactions
• cut2nd (p. 204) is the cut-off for non-bonded interactions calculated only once every nsnb steps
  
• cutprt (p. 204) is the cut-off for non-bonded interactions between the solute and the solvent
• idiel=1 (p. 201) gives a normal constant dielectric constant (epsilon).
• dielc (p. 204) is the value of the dielectric constant.
• scnb (p. 204) is the scale factor for 1-4 van der Waals interactions. scnb=2.0 according to the force field.
• scee (p. 204) is the scale factor for 1-4 electric interactions. scnb=1.2 according to the 1994 force field. (2.0 in the old 1991 force field)
• ntpr (p. 204) is the frequency of printing energies
• ntwx, ntwv, ntwe (p. 204) are the frequencies of writing coordinates, velocities, and energies, respectively, to sampling files. -1 indicates that no such files are written.
• ntb=0 (p. 193) no periodic boundaries
• ntp=0  (p. 195) no constant pressure
• taup (p. 196) time constant for pressure relaxation.
• ibelly=1 (p. 190) indicate that only atoms given in the belly group are allowed to move.
• ivcap (p. 213) controls the cap flag. Default is 0, the cap is read from prmtop. ivcap=1 indicates that a new cap atom pointer should be read from natcap.
• natcap (p. 213) is the new cap atom pointer.

1. Change the sign of isldyn
2. Change almda (to 0.0 or 1.0)

• Without shake: ntc=1, ntf=1, nsnb=50, dt=0.0005
• With shake:    ntc=3, ntf=3, nsnb=15, dt=0.002

Old version (with parm)

1. Start with an equilibrated protein file (sander mdrest file).
 
2. Rerun parm to create a prmtop file for a perturbation run.
You need the edtbin file from edit and a parm.in file (see below for examples).
 
3. Set up the gibbs.in file (see below).
Note that it is not possible to run PMF with shake, therefore you may need to re-equilibrate the coordinates without shake.
 
4. Run gibbs (cgibbs give an input line)
gibbs -i gibbs.in -o gibbs.out -r prest -ms pstat -p prmtop -a comp -c input_coord
-i is the name of the input file
-o is the name of the output file
-r is the name of the restart file; it will contain the final coordinates
-ms is the name of a file with a free energy statistics
-p is the name of the topology file (normally prmtop, but pparm as default)
-a is the name of a file with free energy contributions (if iperat is set)
-c is the name of the file with input coordinates. The first time it should be prmcrd; if you have already run gibbs, you should instead use the restart file (prest) from the previous calculation.

There is one group for each perturbed residue (give the residue name and number) and then one row for each atom in the residue, containing the number, the name, the new type, and the new charge.

An explicit example:
Title
BIN FOR MOD4    PERT
    0    0    0
    1    0    0
PERTURBATION
IP
RES  IP     1
 1  IP    Li  1.0000
END
END