Free-energy calculations with Gibbs in Amber
Note that gibbs is no longer included in Amber9.
It was included in Amber8, but there are problems to compile it.
It is compiled and running in Amber7.
How to start up a gibbs calculations
Free energy perturbation (FEP) or potential of mean force (PMF)
calculations.
- Start from an equilibrated protein structure
- 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
- Check that you have got the right data with
desc x.1.ATOM
- check x
- saveamberparmpert x prmtop prmcrd
- quit
- set up the perturbation run using the template files below
- 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
Notes
- 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.
This is a simple gibbs free energy perturbation input file (gibbs.in)
Title
&cntrl
isldyn=-2,nstmeq=40,nstmul=40,amxmov=0.01,
amxrst=0.05,almdl0=0.0010,dlmax=0.010,
idifrg=1,isande=1,almda=1.0,
idsx0=-1,ielper=0,ncorc=0,intprt=0,
irest=0,ntx=5,init=4,
nrun=100,nstlim=10000,dt=0.002,
temp0=300.0,ntt=1,tautp=0.2,
ntc=3,ntf=3,ishkfl=1,
nsnb=15,cut=15.0,cut2nd=20.0,cutprt=1000.0,
idiel=1,dielc=1.0,
scnb=2.0,scee=1.2,
ntpr=1000,ntwx=-1,ntwv=-1,ntwe=-1,ndmpmc=1,
ntb=0,ntp=0,taup=0.2,
ibelly=0,itimth=0,iperat=0,
intr=0
&end
This is a gibbs.in file for a torsion PMF
Title
&cntrl
isldyn=-2,nstmeq=40,nstmul=40,amxmov=0.01,
almdl0=0.0010,dlmax=0.010,
idifrg=1,isande=1,almda=1.0,
idsx0=-1,ielper=0,ncorc=1,intprt=0,
irest=0,ntx=5,init=4,
nrun=100,nstlim=10000,dt=0.0005,
temp0=300.0,ntt=1,tautp=0.2,
ntc=1,ntf=1,ishkfl=1,
nsnb=50,cut=15.0,cut2nd=20.0,cutprt=1000.0,
idiel=1,dielc=1.0,
scnb=2.0,scee=1.2,
ntpr=1000,ntwx=-1,ntwv=-1,ntwe=-1,ndmpmc=1,
ntb=0,ntp=0,taup=0.2,
ibelly=1,ivcap=1,natcap=4923,itimth=1,iperat=4,
intr=2
&end
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
Description of the gibbs-specific parameters.
- 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
Description of the sander and gibbs parameters
- 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.
To reverse the calculation:
- Change the sign of isldyn
- Change almda (to 0.0 or 1.0)
SHAKE
- 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)
- Start with an equilibrated protein file (sander mdrest file).
- 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).
- 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.
- 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.
This is a minimal parm.in file (no perturbed group)
Title
BIN FOR MOD4 PERT
0 0 0
1 0 0
PERTURBATION
None
END
END
This is a parm.in file with a perturbed group
Ferrochelatase binding MMP
BIN FOR MOD4 PERT
0 0 0
1 0 0
PERTURBATION
MMP
RES MMP 307
1 CHA CD -0.1542
2 HHA HC 0.1157
3 C1A CI -0.0470
4 NA NM -0.0089
5 CCA CT -0.2514
6 HCA1 HC 0.0988
7 HCA2 HC 0.0988
8 HCA3 HC 0.0988
...
Next residue
RES NEW 308
1 CHA CD -0.1542
...
END
END
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