How to start an Amber equilibration
 Form to fill in

New version (with LEaP)

Currently, we prefer to use a combination of Maestro and changepdb, as described on the Maestro home page.
Follow these instructions and then go to point 13 below.

  1. Follow the instructions on the maestro home page and then go to point 13 below.

  2. Start with a pdb file
  3.  
  4. Read all the initial rows in the PDB file (before the coordinates).
  5. Read also the article describing the structure.
    Find out the following
    1. Is the molecule a monomer of a polymer? Are all subunits included in the PDB file?
    2. Are there any Cys-Cys cross links?
    3. Are there any missing parts of the structure (check remark 465, missing residues, but also the start and the end of the file).
    4. Are there any unnormal molecules in the structure (check HET and HETNAM)?
    5. At what pH was the structure collected.
     
  6. If there are unnusual residues in the structure, create for these molecules
    (first check if it is already in /home/bio/Data/Amber):
  7. 1.  prep.in files  
    2.  parameter files  
    Copy the final files to /home/bio/Data/Amber
     
  8. Characterise the protein using changepdb :
  9. 1. Number of atoms and residues
    2. cc: calculate the total charge
    3. cap centre: gives the centre of the protein and the radius.
    4. occ: check which residues have occupation number less than 1.00
    Optionally move the center of mass of the protein to the CAP centre or to another centre (command m), recommended if you will solvate in a water cap.  Do not forget to write the new coordinates to a file with command w.
     
  10. If only a part of the protein shall be included in the simulation. Remove the other parts using changepdb, command cut.

  11. If alternative conformations are present, remove the residue conformation (found in point 4 above) with the lowest occupancy and note which residues are disordered. This can be done automatically with changepdb, command occ.

  12. Ions from the buffer (e.g. SO4) may also be removed.

  13. It is advisable to check the crystal structure with the programs
    WHAT CHECK (pats of WHAT IF) or REDUCE

    whatcheck
    Enter name of pdb file (not processed)

    Output is in file pdbout.txt
    First look for errors (normally few) and warnings (often quite a few)
    grep '# Error' pdbout.txt
    grep '# Warning' pdbout.txt

    Then check in particularly the following two sections in the file pdbout.txt:
    His, Asn, Gln, side chain flips
    Histidine type assignments

    reduce -build -DB /away/bio/Reduce/reduce.3.14.080821.src/reduce_het_dict.txt file.pdb >fileH.pdb 2>reduce.out
    This gives a fully protonated file with optimised positions and possibly flipped rings.

  14. If there are any His residues in the structure, you must determine the protonation status of these:
  15. HID has a proton on ND1
    HIE has a progon on NE2
    HIP has protons on both ND1 and NE2 (i.e. it is positively charged).

    First run the file through PROPKA (http://propka.ki.ku.dk/) and check the output. It will indicate which residues (also Asp, Glu, Lys, Arg, Tyr, Cys) to check, but it often makes mistakes, especially close to unusual residues.

    Run changepdb, command hh.
    This program gives you four types of information for each His residue:
    1. Solvent exposure.
    2. Number of charged residues in the neighbourhood, together with the sum of their charge.
    3. Possible H-bond acceptors and donators close to the HD1 and NE2 atoms.
    4. A series of files his1, his2, ..., containing all atoms within 8 of each His residue. Examine these with rasmol or SwissPDB
    The pKa of His is about 6.5, therefore at pH 7, sligthly less than half of the His residues should be HIP.
    If a nitrogen atom is an acceptor, the other nitrogen must be protonated.
    If the residue is exposed to water it is probably HIP (especially if there is a strong H-bond to one of the atoms), otherewise it is probably neutral, unless it interacts with Asp, Glu, or another negatively charged residue.
    HIP residues are characterised by strong and short hydrogen bonds, often to negatively charged residues.
    If the total charge of the protein is very negative, this may be an indication that many His residues may be HIP. You should also consider the possibility that the crystallographers have mixed up the ND1 and CD2, and the NE2 and CE1 atoms (can only be decided from the H-bond pattern, not from the electron density).
    On the basis of this result, change the name of all HIS residues in the pdb file (changepdb command his).
     
  16. Some further information and more check is obtained by running changepdb command bc (buried charges), without using any charges. It will check if there are any charged groups that are buried in the protein. Normally, they are no or very few, unless they form ionic pairs. If there are any, they can be predicted to have unnormal pKa values.

  17. Change other names in the pdb file, if necessary
  18. Set up the Amber calculation (construct the topology and coordinate files) using the program leap (sample file).
    Run tleap (or xleap)
    1. logfile name_of_logfile     (set the name of a log file)
    2. set default pdbwritecharges on       (ensures that charges are written to pdb files)
    3. x=loadpdb pdb_file      (reads the pdb file)
    4. Add extra bonds (e.g. CYS-CYS crosslinks, e.g.
        bond x.22.SG x.65.SG  (note that the atom names must be capital; note also that leap renumbers the residues, especially if there are several subunits; do savepdb x pdbout and look in the pdbout file to get the appropriate residue numbers).
    5. Solvate the protein, if desired, e.g.
        solvateOct x TIP3PBOX 10
    for an truncated periodic octahedron (best); use 10-25 outside the protein.
        solvatecap x TIP3PBOX {0.0 0.0 0.0} 30  
       (adds a water cap around the origin of radius 30 ; capital letters of TIP3PBOX are nessesary; if you did not move the centre of mass to the cap centre, we need to replace 0.0 0.0 0.0 with the actual cap centre, calculated with changepdb, command cap).
    Use a radius of 6 plus the radius of the protein, obtained by changepdb, command cap.
    6. saveamberparm x prmtop prmcrd  (saves the topology and coordinates)
    7. quit   (stops the program)
    Check in each step that you do not get any error messages and if you do, solve them!
    Additional leap commands can be found below.

    How to avoid charged termini.

    Optionally, you may want to add counter ions.
    This is done with the leap command
    addIons x Na+ 203 Cl- 185
    (done after solvateOct; slow).
    Sometimes, it gives strange results (counter ions in strings). If so, try instead (also much faster):
    addIonsRand x Na+ 203 Cl- 185 4
    However, since it replaces any water molecule with a counter ion, the ions may end up inside the protein.
    This can be avoided by running leap both without and with counter ions and saving pdbout files of both calculations. Then run changepdb, command fci (fix counter ions) on the pdbout file with counter ions and read in the pdbout file without counter ions. The program the overwrites the prmcrd file with counter ions that are more than 4 from the protein (and also from each other).
    You need to read in ion parameters, e.g. with
    loadAmberParams frcmod.ionsjc_tip3p

    If you have any non-standard residues, be careful. Amber does not give you any notice if any of the improper dihedral terms are missing - they are only skipped. This is particularly dangerous for -NH2 groups with GAFF, which are often described by nh-hn bonds, although the x-x-nh-hn dihedral is missing and therefore skipped (x-x-n-hn exists, but sometimes it is too small to keep the group planar if the N charge is too negative). Optimise the molecule in solvent and check the structure. I have sometimes used:
    ca-hn-nh-hn        10.0          180.          2.           ! increased from 1.1

    When leap is done, build a pdb file by changeparm (a pdb file built from telap with savepdb does not always have the same order of the atoms as in the prmtop and prmcrd file!):
    mkpdbout
  19.  
  20. This should already be (better) done with Maestro:
    Optionally, you can improve the positions of the hydrogen atoms with pol_h and gwh.
    Is so, you must put crystal water molecules in a separate file called watpdb, before running leap, and then run leap again with the final result file from gwh.
    Thus:
    1. Move crystal waters to a separate file wat.pdb

    2. Run tleap as above, but do not set pdbwritecharges, do not solvate and do not check:
    x=loadpdb pdb_file
    bond x.22.SG x.65.SG
    savepdb x pdbout
    saveamberparm x prmtop prmcrd
    quit

    3. Run pol_h (it gives the position of the polar hydrogens, but not on water molecules):
  21. pol_h <pdbout >ppdbout
    This usually takes ~1 minute.
    This is the protonated file for the protein.
    Add a proper title to this file (a row starting with REMARK)

    4. Convert the prmtop file to old format:
    new2oldparm <prmtop >prmtop.old

    5. Run gwh (it calculates improved positions of the water hydrogens)
    gwh -w wat.pdb -elstat -p prmtop.old <ppdbout >wpdb
    It usually takes about 5 minutes, and it may crash.
    Add this file to the protonated pdb file:
    cat wpdb >> ppdbout
    remove then END record before the water molecules

    6. After all this run tleap on this file (ppdbout) as discriped in the previous point.

  22. Look at the pdbout file with rasmol (or another program to see that the solvent cap has been placed correct and that the structure looks reasonable.

  23. Characterise the final coordinates (pdbout) with changepdb :
  24. 1. Number of atoms and residues
    2. CCharge: calculate the total charge, now including the actual charges
    3. CAp centre: gives the centre of the protein and the radius.
    4. SD (subunit distance) from the CAP centre.
    5. Short contacts. Make sure that there are no pathological short contacts. A better check is made by changeparm (below).
     
  25. You should also check the coordinates with program changeparm:
    1. Read in the prmtop file
    2. Select Energy and read in the prmcrd file (of type mdrest). Check for large energies and unnormal bonds.

  26. Optionally change the charge of some atoms.
  27. Change to the proper charge in file pdbout.
    Then use changeparm, option Modify (and Write) to include the new charges into prmtop.
    Check that the charges add up to an integer by changepdb, otpion CC.
    Normally we use this algorithm to get a correct total charge (it is done automatically by changepdb if you have a pdbout file, a turbomole (mulliken), gaussian (gauss.out), or resp (resp.out) file with the fitted charges, and a file giving the correspondance between the quantum and MM atoms and a list of the junction atoms (comqum.dat):
    1. For residues, for which all atoms are treated by quantum mechanics, all quantum chemical charges are used directly.
    2. For the other residues, the total charge of the residue in the quantum system should also be the total charge of the amino acid in the protein (except if the residue has a charge outside the quantum system).
    3. The charge on the junction atom (the one that is converted to a hydrogen atom in the quantum system) is modified so that the total charge of the residue becomes correct.
    4. If there are more than one junction atoms in the same residue, the change in charge (from the quantum charge) should be the same for all junction atoms.
    This is the algorithm used by fixcharge in ComQum (perhaps this program can be used).
     
  28. Normally, not needed any longer (if you use belly list):
    If all atoms except hydrogens and waters should be fixed, use changepdb, command fix, to make a belly group for the fixed atoms (this is necessary since edit may reorder the atoms). You also need to know the number of crystal water molecules, since only the non-crystal WAT molecules should be allowed to move freely.
  29. The new list is by default written to file temp, and this group should replace the belly group in all sander.in files.
       
  30. If you use a spherical system (not recommended), you should run the calculations with a somewhat smaller cap radius than that set in the leap.in file.
    This is because the harmonic restrains allow the waters to move somewhat outside the cap radius and because the waters close to the protein are not well packed.
    Therefore, you may get a too low density or a system with vacuum spaces on the surface.
    For a 70 cap, the optimum value was 65 . The restraint energy was then 181 kcal/mol (291 kcal for 64 ; 66 gave empty space in the coordinates, but restraint energies of 84 kcal/mol).
    A guess is that -4 is appropriate for 60-69 ; -3 for 50-59 , -2 for 30-49 , and -1 for smaller systems.

    The cap radius is changed by changeparm, command c.

  31. Octahedral system (recommended):
    Copy the sander.in files from here.

    Spherical system:Construct the sander.in files. Start with the sander.in1 and sander.in2 templates below.
  32. You may want to change the CAP list to include also the crystal water molecules. This is done by:
    ivcap=1,matcap=the_last_non-water_atom
    At least the following two sander jobs should be run:
    1.   120 000 * 2 ps simulated annealing with shake (no change)
    2.   10 000 steps molecular mechanics with shake (change nrun=0, imin=1; remove ",vlimit=20.0,nmropt=1" and all the rows starting with &wt or &rst)
    If the job crashes due to shake:
    COORDINATE RESETTING CANNOT BE ACCOMPLISHED,
    DEVIATION IS TOO LARGE
    you may need to run one or two of these:
    sander.in00. 100 steps molecular mechanics without shake (as 2. and ntc=1,ntf=1, nsnb=50,ntpr=10)
    sander.in0.   1000 steps molecular mechanics with shake (as 2 and ntpr=100)
    sander is started by
    sander -i sander.in -o sander.out -r mdrest -e mden -c prmcrd_or_mdres_n-1
    For example:
    sander -i sander.in1 -o sander.out1 -r mdrest1 -e mden1 -c prmcrd

    sander -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c mdrest1
    sander -i sander.in00 -o sander.out00 -r mdrest00 -e mden00 -c prmcrd
    Run first csander integer,
    where integer is the number of the sander.in1 file (here 1). This will remove files which are not allowed to be present (alternatively, you can start sander with the flag -O.
    -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
    -e is a file with energy statistics
    -c is the name of the file with input coordinates. The first time it should be prmcrd; if you have already run sander, you should instead use the restart file (mdres) from the previous calculation.
    -x is the name of a file with coordinates sampled at each ntvx time step.

    If you use a spherical system with belly, you cannot use pmemd or GPU version.
     
  33. When all the sander jobs are completed, build a pdb file from the final coordinates:
  34. ambanalys
    Read pdbout and mdrest2, and then build a pdb-file.
    Optionally, you may first want to cut the residue numbers from pdb.in into pdbout1 (with emacs) so that they are the same as in the pdb file (in pdbout they are consecutively numbered from 1).
    You can also look at the energies form the output file and get a energy file that can be directly plotted by Gnuplot.
    You can also get the average coordinates and get a r.m.s. analysis from this program (to be compared with B-factors in the crystal coordinates file).
    You can also collect a concentrated result file from the sander.out* files by puramb
     
  35. Characterise the pdb file (changepdb):
  36. 1. Number of atoms and residues
    2. CCharge: calculate the total charge, now including the actual charges
    3. CAp centre: gives the centre of the protein and the radius.
    4. SD (subunit distance) from the CAP centre.
    5. CM
    6. Short contacts. Make sure that there are no pathological short contacts.

How to avoid charged ends of the protein (AMT and OXT)


  1. cp $AMBERHOME/dat/leap/cmd/leaprc.ff03 .
  2. Remove all the rows under addPdbResMap, or only the few rows actually used
Alternatively, if you do not want any charged ends at all, issue the leap command:
clearPdbResMap

If you want only selected termini:

  1. Run command changepdb, command ws on your pdb file (it works only locally, not on aurora or abisko)
    It gives you a line like: seq = {ARG GLY ALA SER GLU ARG TYR LEU VAL PHE PRO ALA LYS ALA VAL ARG PRO VAL THR CYS MET TSA HOH HOH HOH HOH }
  2. Insert this line in your leap.in file, before you read the pdb file
    Also change it so that you add the letter N before a residue that should have a N-terminal and C before a residue that you want to be C-terminal.
    For example:
    seq = {NARG GLY ALA SER GLU ARG TYR LEU VAL PHE PRO ALA LYS ALA VAL ARG PRO VAL THR CYS CMET TSA HOH HOH HOH HOH }
    where I made the first residue N-terminal and the last one C-terminal
  3. Change the
    x=loadPdb file
    command in your leap.in file to
    x=loadPdbUsingSeq file seq


Sample sander.in files for octahedral systems (recommended)

Replace 1111 with the number of the first non-crystal water molecule.

Sometimes you get gas bubbles in you final structure. Then, try to use a larger box of water molecules (larger value for solvateOct).
restraint_wt=10000 gives typical movement of atoms by 0.003 A on average and 0.006 A maximum. Deviations i metal-ligand distances of up to 0.01 A. This is typically acceptable.
restraint_wt=1000 gives typical movement of atoms by 0.03 A on average and 0.3 A maximum. Deviations in metal-ligand distances of up to 0.13 A. This is  too large to be acceptable.

sander.in00
Minimisation without shake.
 &cntrl
  irest=0,ntx=1,
  nstlim=0,dt=0.0005,
  imin=1,maxcyc=1000,drms=0.001,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=1,ntf=1,
  nsnb=40,cut=8.0,dielc=1.0,
  ntpr=100,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
  ntb=1,ntp=0,taup=0.2,
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ntr=1,restraint_wt=10000,
  restraintmask=":1-1111 & !@H="
 &end

sander.in0
10 ps MD simulation without shake
&cntrl
  irest=0,ntx=1,
  nstlim=20000,dt=0.0005,
  imin=0,maxcyc=100,drms=0.001,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=1,ntf=1,
  nsnb=40,cut=8.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
  ntb=1,ntp=0,taup=0.2,
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ntr=1,restraint_wt=10000,
  restraintmask=":1-1111 & !@H="
 &end

sander.in1
1 ns equilibration with shake and constant volume
 &cntrl
  irest=1,ntx=5,
  nstlim=500000,dt=0.002,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=2,ntf=2,
  nsnb=15,cut=8.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
  ntb=1,ntp=0,taup=0.2,
  ipol=0,igb=0,
  ntr=1,restraint_wt=10000,
  restraintmask=":1-1111 & !@H="
 &end

sander.in2
10 ns simulated annealing at constant pressure
 &cntrl
  irest=1,ntx=5,
  nstlim=5000000,dt=0.002,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=2,ntf=2,
  nsnb=15,cut=8.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
  ntb=2,ntp=1,pres0=1.0,taup=1.0,
  ipol=0,igb=0,
  nmropt=1,
  ntr=1,restraint_wt=10000,
  restraintmask=":1-1111 & !@H="
 &end

 &wt type='TEMP0',istep1=     0,istep2=1600000,value1=370.00,value2=370.00  &end

 &wt type='TEMP0',istep1=1600001,istep2=5000000,value1=370.00,value2=  0.00  &end

 &wt type='TAUTP',istep1=     0,istep2=2000000,value1=  0.20,value2=  0.20  &end

 &wt type='TAUTP',istep1=200001,istep2=3200000,value1=  1.00,value2=  1.00  &end

 &wt type='TAUTP',istep1=3200001,istep2=4000000,value1=  0.50,value2=  0.50  &end

 &wt type='TAUTP',istep1=4000001,istep2=5000000,value1=  0.05,value2=  0.05  &end

 &wt type='END' &end

 &rst iat=0 &end

sander.in 3
Minimisation with shake
 &cntrl
  irest=0,ntx=1,
  nstlim=0,dt=0.002,
  imin=1,maxcyc=10000,drms=0.001,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=2,ntf=2,
  nsnb=15,cut=8.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,ioutfm=0,ntxo=1,
  ntb=2,ntp=1,pres0=1.0,taup=1.0,
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ntr=1,restraint_wt=10000,
  restraintmask=":1-1111 & !@H="
 &end


Sample file sander.in1

Test to use ig=-1, ntt=3, igb=10

500 000 x 2 fs MD simulated annealing with SHAKE, only H atoms and WAT free
Title
 &cntrl
  irest=0,ntx=1,
  nstlim=500000,dt=0.002,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=2,ntf=2,
  nsnb=15,cut=15.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2,
  ipol=0,igb=0,
  vlimit=20.0,nmropt=1,
 
ibelly=1,bellymask=":WAT | @H="
 &end

 &wt type='TEMP0',istep1=     0,istep2=160000,value1=370.00,value2=370.00  &end

 &wt type='TEMP0',istep1=160001,istep2=500000,value1=370.00,value2=  0.00  &end

 &wt type='TAUTP',istep1=     0,istep2=200000,value1=  0.20,value2=  0.20  &end

 &wt type='TAUTP',istep1=200001,istep2=320000,value1=  1.00,value2=  1.00  &end

 &wt type='TAUTP',istep1=320001,istep2=400000,value1=  0.50,value2=  0.50  &end

 &wt type='TAUTP',istep1=400001,istep2=500000,value1=  0.05,value2=  0.05  &end

 &wt type='END' &end

 &rst iat=0 &end

 
If you instead use an explicit belly list, it should come directly after &rst.

It is probably better to use periodic boundary conditions instead. Set:
cut=8.0
ntb=1


Sample file sander.in2

10 000 steps MM with SHAKE
Title
 &cntrl
  irest=0,ntx=1,
  nstlim=0,dt=0.002,
  imin=1,maxcyc=10000,drms=0.001,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=2,ntf=2,
  nsnb=15,cut=15.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ibelly=1,bellymask=":WAT | @H="
 &end


Sample file sander.in00

1000 steps MM without SHAKE
Title
 &cntrl
  irest=0,ntx=1,
  nstlim=0,dt=0.0005,
  imin=1,maxcyc=1000,drms=0.001,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=1,ntf=1,
  nsnb=40,cut=15.0,dielc=1.0,
  ntpr=100,ntwx=0,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2,
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ibelly=1,bellymask=":WAT | @H="
 &end


Sample file sander.in0

20000 steps MD (10 ps) without SHAKE
Title
&cntrl
  irest=0,ntx=1,
  nstlim=20000,dt=0.0005,
  imin=0,maxcyc=100,drms=0.001,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=1,ntf=1,
  nsnb=40,cut=15.0,dielc=1.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2,
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ibelly=1,bellymask=":WAT | @H="
 &end

Sample file sander.in3

300000 x 2 fs MD equilibration, without belly
PBC and Langevin thermostat

Title
 &cntrl
  irest=0,ntx=5,
  nstlim=300000,dt=0.002,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=2,ntf=2,
  cut=8.0,
  ntpr=1000,ntwx=0,ntwv=0,ntwe=0,iwrap=1,
  ntb=2,ntp=1,pres0=1.0,taup=1.0,
  ipol=0,igb=0,
  &end

Sample sander.in files for a long MD simulation with PBE (e.g. entropy)

UR 8/10-15

sander.in1
Restrained minimization
  &cntrl
   irest=0,ntx=1,
   imin=1,maxcyc=1000,drms=0.0001,ntmin=2,
   ntc=2,ntf=2,
   cut=8.0,
   ntpr=100,ntwx=0,ntwv=0,ntwe=0,
   ipol=0,igb=0,
   ntr=1,restraint_wt=100,
   restraintmask="!:WA= & !@H="
  &end

sander.in2
Restrained NVT equilibration
  &cntrl
   irest=0,ntx=1,
   nstlim=10000,dt=0.002,
   temp0=300.0,ntt=3,gamma_ln=2.0,
   ntc=2,ntf=2,
   cut=8.0,
   ntpr=500,ntwx=0,ntwv=0,ntwe=0,
   ntb=1,
   ipol=0,igb=0,
   ntr=1,restraint_wt=50,
   restraintmask="!:WA= & !@H="
  &end

sander.in3
Restrained NPT equilibration
  &cntrl
   irest=1,ntx=5,
   nstlim=10000,dt=0.002,
   temp0=300.0,ntt=3,gamma_ln=2.0,
   ntc=2,ntf=2,
   cut=8.0,
   ntpr=500,ntwx=0,ntwv=0,ntwe=0,
   ntb=2,ntp=1,pres0=1.0,taup=1.0,
   ipol=0,igb=0,
   ntr=1,restraint_wt=50,
   restraintmask="!:WA= & !@H="
  &end

sander.in4
Equilibration, 1 ns NPT
  &cntrl
   irest=1,ntx=5,
   nstlim=500000,dt=0.002,
   temp0=300.0,ntt=3,gamma_ln=2.0,
   ntc=2,ntf=2,
   cut=8.0,
   ntpr=500,ntwx=0,ntwv=0,ntwe=0,iwrap=1,
   ntb=2,ntp=1,pres0=1.0,taup=1.0,
   ipol=0,igb=0,
   ntr=0
  &end

sander.in5
Production 20 ns
  &cntrl
   irest=1,ntx=5,
   nstlim=10000000,dt=0.002,
   temp0=300.0,ntt=3,gamma_ln=2.0,
   ntc=2,ntf=2,
   cut=8.0,
   ntpr=500,ntwx=500,ntwv=0,ntwe=0,iwrap=1,
   ntb=2,ntp=1,pres0=1.0,taup=1.0,
   ipol=0,igb=0,
   ntr=0
  &end


Old file, without PBC and with Berendsen thermostat (obsolete)

Title
 &cntrl
  irest=0,ntx=5,
  nstlim=300000,dt=0.002,
  imin=0,maxcyc=10000,drms=0.001,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=2,ntf=2,
  nsnb=25,cut=15.0,dielc=1.0,
  ntpr=1000,ntwx=1000,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2,
  ipol=0,igb=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ibelly=0
 &end

Sample file sander.in4 with PBC  and Langevin thermostat

2 ns MD sampling every 5 ps without belly
Title
 &cntrl
  irest=0,ntx=5,
  nstlim=1000000,dt=0.002,
  temp0=300.0,ntt=3,gamma_ln=2.0,
  ntc=2,ntf=2,
  cut=8.0,
  ntpr=2500,ntwx=2500,ntwv=0,ntwe=0,iwrap=1,
  ntb=2,ntp=1,pres0=1.0,taup=1.0,
  ipol=0,igb=0,
  &end



Sample file for sander with restraints

With belly, no PBC and Berendsen, i.e. oldfashioned.
Title
 &cntrl
  irest=1,ntx=5,
  nstlim=2000000,dt=0.0005,
  imin=0,maxcyc=10000,drms=0.001,
  temp0=300.0,ntt=1,tautp=0.2,
  ntc=1,ntf=1,
  nsnb=50,cut=10.0,dielc=1.0,
  ntpr=2000,ntwx=0,ntwv=0,ntwe=0,
  ntb=0,ntp=0,taup=0.2,
  ipol=0,
  ncyc=10,ntmin=1,dx0=0.01,
  ibelly=1,ivcap=1,matcap=9902,
  nmropt=1
 &end

 &wt type='END' &end

 &rst
  iat(1)=396,iat(2)=9005,iat(3)=0,iat(4)=0,
  r1=1.5,r2=2.5,r3=3.5,r4=4.5,
  rk2=50.0,rk3=50.0
 &end

 &rst iat=0 &end

GRP1
ATOM   323   465
ATOM   498   537
END
END

Alternatively (this is the only thing that works recently):
...
&end

  &wt type='END' &end
  DISANG=rst

rst file:
&rst iat=17892, 865, r1=0., r2=2.16, r3=2.16, r4=100., rk2=100., rk3 =100.,/
&rst iat=17892, 7802, r1=0., r2=2.20, r3=2.20, r4=100., rk2=100., rk3 =100.,/



Sander parameters (with page references for the Amber 8 manual)

Commands removed in Amber 7.0 (pages in the Amber 6.0 manual), but still present in Gibbs

Sample of leap commands

source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
loadAmberPrep hpo.in

loadAmberParams o2.dat
x=loadpdb 1icf4
bond x.156.SG x.206.SG
solvatecap x TIP3PBOX {0.0 0.0 0.0} 33
check x

saveamberparm x prmtop prmcrd
quit

A sodium counter ion has should have both atom and residue names Na+ (and atom type IP).

To avoid errors like:
+Currently only Sp3-Sp3/Sp3-Sp2/Sp2-Sp2 are supported
+---Tried to superimpose torsions for: *-C1-C6-*
+--- With Sp0 - Sp0
+--- Sp0 probably means a new atom type is involved
+--- which needs to be added via addAtomTypes
Insert in leap.in (before the in/dat files are read in):
addAtomTypes {
  {"c1" "C" "sp2"}
  {"c2" "C" "sp2"} }

To avoid promlems like:
Bond: maximum coordination exceeded on .R<CUA 465>.A<CU 1>
-- setting atoms pert=true overrides default limits

First try to delete any CONECT entries in the pdb file.
If this does not work, use in leap.in:
set x.465.CU pert true

Sometimes a rectangular box is preferred
The command is then
solvateBox x TIP3Pbox 10
It can then be favourable to rotate the system before, to align the axes of inertia with the coordinate axes.
This can be done inside tleap with the transform command (no, this does not work!):
transform x {
{1 0 0 }
{0 cos q -sin q}
{0 sin q cos q} }
Where q is the rotation angle (but you need to give the proper numbers)

Do it instead in the pdb file before reading it into tleap, i.e. with changepdb, command tr
with a file like this:
0.866000 0.500000 0.000000 0.00000 -0.500000 0.866000 0.000000 0.00000 0.000000 0.000000 1.000000 0.00000


Some useful Leap commands


How to restrain a metal site
TI 1 min, lambda=0.500
  &cntrl
   irest=0,ntx=1,
   imin=1,maxcyc=500,drms=0.0001,ntmin=2,
   ntc=2,ntf=1,
   cut=8.0,
   ntpr=100,ntwx=0,ntwv=0,ntwe=0,
   ntr=1,restraint_wt=100,
   restraintmask="!:WA= & !@H=",
   icfe=1,clambda=0.500,ifsc=1,
   timask1=":C47",scmask1=":C47@H9",
   timask2=":C05",scmask2=":C05@Br9",
   nmropt=1
  &end

  &wt type='END' &end
  DISANG=rst

The file rst (but can also be inserted in the sander.in file)
&rst
    iat=4427,1612,
    r1=0.0,r2=2.4,r3=2.4,r4=100.,
    rk2=50.0,rk3=50.0
   &end

&rst
    iat=4427,1653,
    r1=0.0,r2=2.4,r3=2.4,r4=100.,
    rk2=50.0,rk3=50.0
   &end

&rst
    iat=4427,1699,
    r1=0.0,r2=2.4,r3=2.4,r4=100.,
    rk2=50.0,rk3=50.0
   &end

&rst
    iat=4427,1755,
    r1=0.0,r2=2.4,r3=2.4,r4=100.,
    rk2=50.0,rk3=50.0
   &end

&rst
    iat=4427,5102,
    r1=0.0,r2=2.4,r3=2.4,r4=100.,
    rk2=50.0,rk3=50.0
   &end

&rst
    iat=4427,5166,
    r1=0.0,r2=2.4,r3=2.4,r4=100.,
    rk2=50.0,rk3=50.0
   &end

  &wt type='END', &end
  &rst iat=0, &end


A comment on NPT calculations with restraints
It seems that if you run NPT simulations with restraints, the coordinates of the reference structure is also scaled as the NPT box is scaled to keep the pressure constant. This means that if you restart the same calculation with the the same restraint file (e.g. -ref prmcrd), you will get a discontinuity at the restart, because then the restraint file is suddenly no longer scaled at the restart, so in practice a different restrain structure is used. 

Old versions (prep - link - edit - parm)

  1. Start with a pdb file

  2.  
  3. Read all the initial rows in the PDB file (before the coordinates).

  4. Read also the article describing the structure.
    Find out the following
    1. Is the molecule a monomer of a polymer? Are all subunits included in the PDB file?
    2. Are there any Cys-Cys cross links?
    3. Are there any missing parts of the structure (check the start and the end of the file).
    4. Are there any unnormal molecules in the structure?
    5. At what pH was the structure collected.
     
  5. If there are unnusual residues in the structure, create

  6. 1.  protonate entries (in $OML_AMBER/dat/Mumod/PROTON_INFO.Ulf)
    2.  prep.in files   (in $OML_AMBER/dat/Mumod directory)
    3.  parameter files   (in $OML_AMBER/dat/Mumod/modparm.dat)
     for these molecules.
     
     
  7. Characterise the protein using changepdb :

  8. 1. Number of atoms and residues
    2. cc: calculate the total charge
    3. cap centre: gives the centre of the protein and the radius.
    4. occ: check which residues have occupation number less than 1.00
    Optionally move the center of mass of the protein to the CAP centre or to another centre (command m).
     
  9. Remove the residue conformation (found in point 4 above) with the lowest occupancy and note which residues are disordered.

  10.  
  11. Ions from the buffer (e.g. SO4) may also be removed

  12.  
  13. If there are any His residues in the structure, you must determine the protonation status of these. Run  changepdb, command hh

  14. and give the name of the protein pdb file on the next line (no promt).
    Analyse the resulting file to decide which nitrogens on the His residues are protonated.
    HID has a proton on ND1
    HIE has a progon on NE2
    HIP has protons on both ND1 and NE2 (i.e. it is positively charged)
    The pKa of His is about 6.5, therefore at pH 7, sligthly less than half of the His residues should be HIP.
    The file shows possible acceptors and donators to the His nitrogens.
    If a nitrogen is an acceptor, the other nitrogen must be protonated.
    You may also look at the closest neighbourhood of the His residues with rasmol (with changepdb, command cut, you can cut our only the surroundings of each His residue).
    If the residue is exposed to water, it is probably HIP, otherwise it is probably neutral, unless it interacts with Asp, Glu, or another negatively charged residue.
    HIP residues are characterised by strong and short hydrogen bonds, often to negatively charged residues.
    If the total charge of the protein is very negative, this may be an indication that many His residues may be HIP.
    On the basis of this result, change the name of all HIS residues in the pdb file (changepdb command his).
     
  15. Change other names in the pdb file, if necessary
  16. Move crystal water molecules to the file watpdb

  17.  
  18. Protonate the structure (this is done in the following seven steps).

  19.  
  20. Run protonate (it gives the coordinates of all hydrogens but those of polar hydrogens are not optimal):

  21. $OML_AMBER/exe/protonate -k -d $OML_AMBER/dat/Mumod/PROTON_INFO.Ulf <pdb_file >ppdb_file
     
  22. Run pdbtoamber

  23. 0. The log file is a file with the choices you have done. The default name (logfil) is OK.
    1. Do not check for short contacts.
    2. Give a title if it is missing
    3. Select Amber amino acids
    4. Select both
    5. Select one atom as the cap centre (any one).
    6. Select the whole protein in the all-atom system (give a cut-off of 1000)
    7. No uni-atom system (cut-off=0), no fixed system (cut-off=0), no point-charges, no cap, no perturbation, no simulated annealing, no belly, and no analysis files.

    This gives the files link.in, edit.in, pdb.in, parm.in, and sander.in
     

  24. Ensure that the link.in file is complete.

  25. 1. All non-standard residues are listed at the start of the file.
     

    2. Crosslinks (e.g. Cys-Cys bonds) are included as appropriate.
        This is done by change the first zero on the row before the amino acids to 1 and then adding a description of the crosslink before the end of the file in this format: '   30   51SG  SG      1', where the first two numbers are the amino acid numbers, followed by the two atoms involved in the crosslink, and finally the molecule from which the the first atom come (it may be another than the current).
    See the link example below.
     

    3. Protein ends are charged, if applicable (if you have a whole protein, the first residue is positively charged (NH3+) and the last residue is negatively charged (COO-).
    You indicate this by setting in the row before the amino acids the first letter should be P (protein) and the last digit should be 1. Thus: 'P   0    0    1    3    1' or 'P   1    0    1    3    1' depending on whether you have crosslinks or not.
      It may be necessary also to divide the system into several subsystems (protein + ligands).
     

    4. Ligands are separated by ***.
     

    Then run link
    clink
    This gives files link.out and lnkbin
    Check link.out that:
    1. the job was sucessful
    2. the charge of each molecule is correct.
     

  26. Ensure that the edit.in and pdb.in files are complete.

  27. This is normally the case at this point.
    Then run edit:
    ced
    This gives edit.out, edtbin, and pdbout
    Check edit.out that:
    1. There were no missing or extra atoms in the structure (beginning of the file, between the two words XRAY)).
    2. The five longest bonds are correct (end of the file)
     
  28. Ensure that the parm.in file is correct

  29. This is normally the case.
    Then run parm
    cparm
    This gives files parm.out, prmcrd, and prmtop.
    Check parm.out that the job was sucessful; the last line should be:
    PARM-I-NOPROB, successful completion of parm program.
     
  30. Run pol_h (it gives the position of the polar hydrogens, but not on water molecules):

  31. $OML_AMBER/exe/pol_h <pdbout >hppdb_file
    This usually takes about 5 minutes.
    This is the final protonated file for the protein.
    Add a proper title to this file (a row starting with REMARK ).
     
  32. Run gwh (it calculates the positions of the water hydrogens):

  33. $OML_AMBER/exe/gwh <the_protonated_file >wpdb
    It expects to find the coordinates of the oxygen atoms in file watpdb.
    It usually takes about 5 minutes, and it may crash (then use the addwat option of changepdb on file watpdb and write them to file addfil).
    grep -v TER wpdb | grep -v END >addfil
    This is the final protonated file for the water atoms.
    It is read in by edit below, and this program expects the file name to be addfil, otherwise it must be explicitly defined.
     
  34. If there are symmetry-related subunits, create these using the symmetry information (make a file with the rotation matrix and the translation vector) and the program changepdb, command DIMer.

  35.  
  36. Optionally run changepdb to

  37. 1. Centre the coordinates around a proper atom (move)
    2. Calcuate the cap centre.
    This must be done on both the protein and the water atoms. Therefore you first have to insert the water molecules (from file addfil) in the protein coordinate file, then run changepdb, and then move the resulting coordinates of the water molecules back to  addfil.
     
  38. Make copies of the link.in, edit.in, and parm.in files if you have done any changes to them. They will be destroyed in the next step.

  39.  
  40. Define the system to equilibrate, including the cap centre and the fixed-atom system.

  41. pdbtoamber
    1. Do not check for short contacts.
    2. Change the title if you want.
    3. Select Amber amino acids.
    4. Select both.
    5. Select the proper cap centre atom,

    6. either the centrum of the protein if you can afford to include the full protein, or
      the central atom of the interesting system (active site, site of perturbation, quantum system, etc.). The cap centre is the centre of a sphere of water molecules solvating the protein. It can be changed to the origin in the edit.in file.
    7. Select the radius within which all atoms should be included

    8. preferably, the whole protein (radius = 1000 ),
      otherwise a radius of  25-30 ).
    9. Set radius to 0 for the united-atom system
    10. Select the fixed-atom system (hydrogens for an equilibration of a crystal structure).
    11. No point charges.
    12. Set radius for the CAP (typically 25-30 ).
    13. No perturbation
    14. Select belly.
    15. Select simulated annealing.
    16. No analysis files.
    The files link.in, edit.in,pdb.in, parm.in, and sander.in are made.
       
  42. Ensure that the link.in file is complete.

  43. This was probably done already above.
    1. All non-standard residues are listed.
    2. Crosslinks are included as appropriate.
    3. Protein ends are charged, if applicable.
    4. Ligands are separated by ***.
    Then run link
    clink
    This gives files link.out and lnkbin
    Check link.out that:
    1. the job was sucessful
    2. the charge of each residue is correct.
     
  44. Ensure that the edit.in and pdb.in files are complete.

  45. In particular, ensure that the SOL statement is appropriate. You may want to change the last number in the third row to -1 so that the origin is the cap centre. Then remove also the cap centre: SOL
    HW   OW       4
       0.41700   1.80000   1.80000
        0    0    0   -1
      24.00000   0.00000   0.00000   0.00000 Add the ADD section just before the SOL section: ADD
       0.41700HW  OW This ensures that the water molecules are read in and are treated by the fast-water algorithms in the simulations.
    Then run edit:
    ced
    This gives edit.out, edtbin, and pdbout
    Check edit.out that:
    1. There were no missing or extra atoms in the structure (beginning of the file).
    2. The five longest bonds are correct (end of the file)
     
  46. Ensure that the parm.in file is correct

  47. This is normally the case.
    Then run parm
    cparm
    This gives files parm.out, prmcrd, and prmtop.
    Check parm.out that the job was sucessful.
     
  48. Characterise the final coordinates (pdbout) with changepdb :

  49. 1. Number of atoms and residues
    2. CCharge: calculate the total charge, now including the actual charges
    3. CAp centre: gives the centre of the protein and the radius.
    4. SD (subunit distance) from the CAP centre.
    5. Short contacts. Make sure that there are no pathological short contacts. A better check is made by changeparm (below).
     
  50. You should also check the coordinates with program changeparm.
    1. Read in the prmtop file
    2. Select Energy and read in the prmcrd file (of type mdrest). Check for large energies and unnormal bonds.
  51. Optionally change the charge of some atoms.

  52. Change to the proper charge in file pdbout.
    Then use changeparm, option Modify (and Write) to include the new charges into prmtop.
    Check that the charges add up to an integer by changepdb, otpion CC.
    Normally we use this algorithm to get a correct total charge:
    1. For residues, for which all atoms are treated by quantum mechanics, all quantum chemical charges are used directly.
    2. For the other residues, the total charge of the residue in the quantum system should also be the total charge of the amino acid in the protein (except if the residue has a charge outside the quantum system).
    3. The charge on the junction atom (the one that is converted to a hydrogen atom in the quantum system) is modified so that the total charge of the residue becomes correct.
    4. If there are more than one junction atoms in the same residue, the change in charge (from the quantum charge) should be the same for all junction atoms.
    This is the algorithm used by fixcharge in ComQum (perhaps this program can be used).
     
  53. If all atoms except hydrogens and waters should be fixed, use changepdb, command Fix, to make a new belly group for the fixed atoms (this is necessary since edit may reorder the atoms). You also need to know the number of crystal water molecules, since only the non-crystal WAT molecules should be allowed to move freely.

  54. The new list is by default written to file temp, and this group should replacy the belly group in all sander.in files.
  55. Modify the sander.in files.

  56. You may want to change the CAP list to include also the crystal water molecules. This is done by:
    ivcap=1,matcap=the_last_non-water_atom
    At least the following two sander jobs should be run:
    1.   30 000 * 1.5 ps simulated annealing with shake (no change)
    2.   10 000 steps molecular mechanics with shake (change nrun=0, imin=1; remove ",vlimit=20.0,nmropt=1" and all the rows starting with &wt or &rst)
    If the job crashes due to shake:
    COORDINATE RESETTING CANNOT BE ACCOMPLISHED,
    DEVIATION IS TOO LARGE
    you may need to run one or two of these:
    00. 100 steps molecular mechanics without shake (as 2. and ntc=1,ntf=1, nsnb=50,ntpr=10)
    0.   1000 steps molecular mechanics with shake (as 2 and ntpr=100)
    sander is started by
    $OML_AMBER/exe/sander -i sander.in -o sander.out -r mdrest -e mden -c prmcrd_or_mdres_n-1
    sander -i sander.in1 -o sander.out1 -r mdrest1 -e mden1 -c prmcrd
    sander -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c mdrest1
    sander -i sander.in00 -o sander.out00 -r mdrest00 -e mden00 -c prmcrd
    Run first csander integer,
    where integer is the number of the sander.in1 file (here 1). This will remove files which are not allowed to be present (alternatively, you can start sander with the flag -O.
    -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
    -e is a file with energy statistics
    -c is the name of the file with input coordinates. The first time it should be prmcrd; if you have already run sander, you should instead use the restart file (mdres) from the previous calculation.
    -x is the name of a file with coordinates sampled at each ntvx time step.
     
  57. When all the sander jobs are completed, build a pdb file from the final coordinates:

  58. ambanalys
    Read pdbout and mdrest2, and then build a pdb-file.
    Optionally, you may first want to cut the residue numbers from pdb.in into pdbout1 (with emacs) so that they are the same as in the pdb file (in pdbout they are consecutively numbered from 1).
    You can also look at the energies form the output file and get a energy file that can be directly plotted by Gnuplot.
    You can also get the average coordinates and get a r.m.s. analysis from this program (to be compared with B-factors in the crystal coordinates file).
    You can also collect a concentrated result file from the sander.out* files by puramb
     
  59. Characterise the pdb file (changepdb):

  60. 1. Number of atoms and residues
    2. CCharge: calculate the total charge, now including the actual charges
    3. CAp centre: gives the centre of the protein and the radius.
    4. SD (subunit distance) from the CAP centre.
    5. CM
    6. Short contacts. Make sure that there are no pathological short contacts.


Example files

link.in
Title

CU4      1/usr/local/lib/Amber/dat/Mumod/cu4.dat

DU
    1    0    0    0
COMT
P   0    0    1    3    1
ALA 2*** 2LYS 2*** 2ILE 2GLN 2GLY 2ASN 2ASP 2GLN 2MET 2GLN 2PHE 2ASN 2THR 2ASN 2
*** 2LEU 2SER 2HIP 2PRO 2GLY 2ASN 2LEU 2PRO 2LYS 2ASN 2VAL 2MET 2GAA 2HA6 2ASN 2

Copper
O   1    0    1    3    0
CU4 2

   37    1ND1 CU      1
   84    1SG  CU      1
   87    1ND1 CU      1
   92    1SD  CU      1

QUIT


edit.in
Title
   0    0    0    0
XYZ
OMIT
XRAY
    0    0    0    0
ADD
   0.41700HW  OW
SOL
HW   OW       4
   0.41700   1.80000   1.80000
    0    0    0    1
  25.00000   0.00000   0.00000   0.00000
CAP1
ATOM 3383 3383
END
END
XRAY
    2    0    0    0
QUIT
 

parm.in
Title
BIN FOR MOD4
    0    0    0
    1    0    0

How to build and use a solvation box file from a MD simulation
  1. Run the MD simulation of the pure solvent.
  2. Build a pdb file from the last mdrest file with changeparm command p
  3. Build an off file from the pdb file.
  4. Change the box size back to that in the mdrest file (last line) in the off file
  5. Use the off file.

Leap input to make a silly solvation of the pure solvent:
addPath /home/ulf/Iria/amber/Leap
source leaprc.gaff2
loadAmberPrep acn.in
solv=loadpdb acn.pdb
x=loadpdb acn.pdb
solvateBox x solv 25
saveamberparm x prmtop prmcrd
quit

Leap input file to make the off box:
addPath /home/ulf/Iria/amber/Leap
source leaprc.gaff2
loadAmberPrep dcm.in
dcm_box=loadpdb pdb4
setBox dcm_box centers
saveOff dcm_box acn_box.off
quit

Leap input file to use the box file:
addPath /home/ulf/Iria/amber/Leap
source leaprc.gaff2
loadAmberPrep acn.in
loadoff acn_box.off
x=loadpdb acn.pdb
solvateBox x acn_box 10
saveamberparm x prmtop prmcrd
quit


Some simple ptraj scripts
cpptraj prmtop ptraj.in

Wrap water molecules:
trajin mdcrd5
trajout mdcrd5_wrap trajectory
image origin center familiar com :411
go


Format an mdrest file
trajin mdrest5
trajout mdrest5_form restart
go