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.

 
  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


Notes

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.

Description of the sander and gibbs parameters

To reverse the calculation:
  1. Change the sign of isldyn
  2. Change almda (to 0.0 or 1.0)


SHAKE

Old version (with parm)

  1. Start with an equilibrated protein file (sander mdrest file).
  2.  
  3. Rerun parm to create a prmtop file for a perturbation run.
  4. You need the edtbin file from edit and a parm.in file (see below for examples).
     
  5. Set up the gibbs.in file (see below).
  6. Note that it is not possible to run PMF with shake, therefore you may need to re-equilibrate the coordinates without shake.
     
  7. Run gibbs (cgibbs give an input line)
  8. 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