New version (with LEaP)

1. Follow the instructions on the maestro home page and then go to point 13 below.
   
   
2. Start with a pdb file
   It is normally wise to delete all ANISOU and CONECT lines in the file.
   grep -v ANISOU file | grep -v CONECT > newfile
   
 
3. Read all the initial rows in the PDB file (before the coordinates).
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.
 
4. If there are unnusual residues in the structure, create for these molecules
   (first check if it is already in /home/bio/Data/Amber):
   
1.  prep.in files  
2.  parameter files  
Copy the final files to /home/bio/Data/Amber
 
5. Characterise the protein using changepdb :
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.
 
6. If only a part of the protein shall be included in the simulation. Remove the other parts using changepdb, command cut.
   
   
7. 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.
   


8. Ions from the buffer (e.g. SO4) may also be removed.
   
   
9. 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.
   
   
10. If there are any His residues in the structure, you must determine the protonation status of these:
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).
 
11. 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.
    
    
12. Change other names in the pdb file, if necessary
• Metals
• CYS->CYM, if it is a metal ligand
• CYS->CYX, if it is a Cys-Cys crosslink
 
13. 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.
    solvateBox x TIP3PBOX 10
    for a rectangular system (currently preferred); use 10-25 Å outside the protein.     solvateOct x TIP3PBOX 10
    for an truncated periodic octahedron (good, but often problematic with Amber); 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!
    The most common errors are caused by missing TER commands between molecules that should not be connected.
    Sometimes CONECT lines at the end of the pdb file may cause similar problems.
    Always delete all CONECT lines.
    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
    
 

14. 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):
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.


15. 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.
    
    
16. Characterise the final coordinates (pdbout) with changepdb :
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).
 
17. 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.
   
   
18. Optionally change the charge of some atoms.
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).
 
19. 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.
The new list is by default written to file temp, and this group should replace the belly group in all sander.in files.

 
20. 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.
    
    
21. Rectangular system (recommened):
    Copy the sander.in files from here.
    
    Octahedral system (good, but often problematic with Amber):
    Copy the sander.in files from here.
    
    Spherical system:Construct the sander.in files. Start with the sander.in1 and sander.in2 templates below.
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.
 
22. When all the sander jobs are completed, build a pdb file from the final coordinates:
changeparm
Read pdbtop 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.
 
23. Characterise the pdb file (changepdb):
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.

By the end of the simulation, the protein most likely will have moved outside the simulation box, which can create artefacts in the visualization of the resulting trajectory. Furthermore, it will rotate during the simulation, and visualization of a restricted region, such as the active site, can become awkward. 

To solve this, you have to recentre the protein in the simulation box and eliminate the rotations. Here, is a simple procedure that does this by aligning the protein to a reference structure. 

First, create a file script_recenter.in with the following instructions:

and run cpptraj with

ccptraj -i script_recenter.in

(Marcos 29/10-23)

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
   

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 rectangular systems (recommended)

Sometimes you get gas bubbles in you final structure. If so:

1. Then, build a pdb file of the final structure (pdb.in).
2. Read it into leap and fill the bubbles with a command like
   solvateBox x TIP3PBox 0 0.8
3. mkpdbout
4. Use changepdb, command rwb, to remove water molecules outside the original box (read in the pdb.in file to illustrate the original box). Write it out as pdb.in2
5. Run leap again with pdb.in2 and with
   setBox x centers
   instead of solvateBox
6. Replace the last line in prmcrd with the corresponding line from the mdrest file of the previous calculation (giving gas bubbles).
7. Rerun pmemd with these files.
   Normally, this is enough to get rid of the gas bubbles.
   

 
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.

Sample sander.in files for octahedral systems (recommended)

Neraly always 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.

Sample file sander.in1

 &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

Sample file sander.in00

Sample file sander.in0

Sample file sander.in3

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

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

Sample file sander.in4 with PBC  and Langevin thermostat

Sample file for sander with restraints

• irest=1 (p. 95) indicate a restart.
• ntx=5 (p. 95) indicates that both coordinates and velocities should be read from the prest file. ntx=1 only coordinates.
• nstlim (p. 102) is the number of MD steps per nrun
• dt (p. 103) is the time step in ps. Use 0.001-0.002 with SHAKE and 0.0005-0.001 without.
• imin=1 (p. 95) tells the program to do a minimisation
• maxcyc (p. 102) is the number of cycles in the minimisation
• drms (p. 102) is the convergence criterion ofor the energy gradient; default 0.0001 kcal/mole/Å
• temp0 (p. 103) is the desired temperature.
• ntt=1 (p. 103) is the defaut Berendsen temperature coupling algorithm.
• tautp (p.104) is the temperature relaxation time. Should be 0.1-0.4, default 0.1.
• ntc (p. 105) is the SHAKE flag. ntc=1 no shake. ntc=2 use shake only for H atoms. ntc=3 use shake for all bonds.
• ntf (p. 98) determines which forces to calculate. Should be set to the same value as ntc.
• nsnb (p. 99) 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 shake
• cut (p. 98) is the cut off for non-bonded interactions
• dielc (p. 98) is the value of the dielectric constant.
• ntpr (p. 96) is the frequency of printing energies
• ntwr (p. 96) is the frequency of printing the restart file
  
• ntwx, ntwv, ntwe (p. 96-97) are the frequencies of writing coordinates, velocities, and energies, respectively, to sampling files. 0 indicates that no such files are written.
• ntb=0 (p. 98) no periodic boundaries. ntb=1 constant volume. ntb=2 constant pressure.
  
• ntp=0  (p. 105) no constant pressure. 1= with constant pressure.
• pres0 (p. 105) the desired pressure (in bar).
  
• taup (p. 105) time constant for pressure relaxation.
• ipol=0 (p. 99) polarizabilities anre not included
• scnb (p. 98) is the scale factor for 1-4 van der Waals interactions. scnb=2.0 according to the force field. Not allowed in Amber >12
  
• scee (p. 99) is the scale factor for 1-4 electric interactions. scee=1.2 according to the 1994 force field. (2.0 in the old 1991 force field). Not allowed in Amber >12
• ncyc (p. 102) is the number of initial steepest decent cycles in the energy minimisation, default 10
• ntmin (p. 102) is the algorithm for the energy minimisation: 0 = conjugate gradient, 2 = steepest decent; 1 = a combination (c.f. ncyc).
• dx0 (p. 102) is the initial step length in energy minimisations
• dxm (p. 102) is the maximum step length in energy minimisations; No longer allowed
  
• ibelly=1 (p. 101) indicate that only atoms given in the belly group are allowed to move.
  Specified with either a group or BELLYMASK (p. 243).
  
• ntr (p. 101) = 1 indicates that atoms specified in RESTRAINTMASK (p. 243) are restrained with the force constant RESTRAINT_WT.
  
• vlimit (p. 104) is an upper limit of the velocities. Do not use except at high temperatures.
• nmropt>0 (p. 95) tells the program that there will be more (NMR-type or simulated annealing)input below (&wt, &rst, etc.)
• ivcap (p. 106) 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 matcap (not present in Amber 8 manual). ivcap=2 => cap inactivated.
  
• matcap  (not present in Amber 8 manual) is the new cap atom pointer.
• igb (p. 99) Generalised Born method: igb=0 not used.

• cut2nd (p. 139) is the cut-off for non-bonded interactions calculated only once every nsnb steps
• idiel=1 (p. 137) gives a normal constant dielectric constant (epsilon).
• nrun (p. 142) 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.
• init (p. 142) is the way to obtain velocities. init=4 should be combined with ntx=5, whereas init=3 should be used with ntx=1

Sample of leap commands

Insert in leap.in (before the in/dat files are read in):

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

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

• addPath directory  -  add a directory to the path where files are looked for. Can also be done from the command line: tleap -I directory
• addIons x Na+ 203 Cl- 185 - add 203 Na+ ions and 185 Cl- ions according to electrostatic potential (replacing water molecules, if existing); may give strange results - check! Slow procedure employing a distance-dependent dielectric constant.
  
• addIonsRand x Na+ 203 Cl- 185 4 - add counter ions that are more than 4 Å apart by replacing random water molecules (need to be added before). Changepdb command fci may remove water molecules inside the protein. You probably need to say also: loadAmberParams frcmod.ionsjc_tip3p or similar to get ion parameters.
  
• alias command alias  - make an alias of a command
• bond atom atom - make a bond between two atoms (atom format example: x.22.SG)
• deleteBond atom atom - delete a bond
  
• desc unit/residue/atom - describe an object
• loadAmberParams file  - load new parameters
• loadAmberPrep file   - loads a prep file
  
• logfile filename - specify the name of a logfile
• x=loadPdb file - read in a pdb file
• quit   - quit the program
• saveAmberParm unit prmtop prmcrd - save the topology and coordinates of an unit
• saveAmberParmPert unit prmtop prmcrd - save the topology and coordinates of an unit for a perturbation calculation
• savePdb unit filename - save an object as a pdb file
• set default/object parameter object - set some parameters
  
• solvateCap solute TIP3PBOX position radius - solvate the solute/object in a sphere of water molecules centered at position (an atom or three coordinates, e.g. {0.0 0.0 0.0}) with the radius radius.
• source file - read in a command (or leaprc)  file. Can also be done directly: tleap -f file

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,ntwr=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

Old versions (prep - link - edit - parm)

1. Start with a pdb file

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

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.
 
3. If there are unnusual residues in the structure, create

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.
 
 
4. Characterise the protein using changepdb :

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).
 
5. Remove the residue conformation (found in point 4 above) with the lowest occupancy and note which residues are disordered.

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

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

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).
 
8. Change other names in the pdb file, if necessary
• Metals
• CYS->CYM, if it is a metal ligand
• CYS->CYX, if it is a Cys-Cys crosslink

 
9. Move crystal water molecules to the file watpdb

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

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

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

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
 

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

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.
 

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

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)
 
15. Ensure that the parm.in file is correct

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.
 
16. Run pol_h (it gives the position of the polar hydrogens, but not on water molecules):

$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 ).
 
17. Run gwh (it calculates the positions of the water hydrogens):

$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.
 
18. 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.

 
19. Optionally run changepdb to

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.
 
20. 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.

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

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,

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.
6. Select the radius within which all atoms should be included

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

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.
 
23. Ensure that the edit.in and pdb.in files are complete.

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)
 
24. Ensure that the parm.in file is correct

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.
 
25. Characterise the final coordinates (pdbout) with changepdb :

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).
 
26. 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.
27. Optionally change the charge of some atoms.

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).
 
28. 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.

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

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.
 
30. When all the sander jobs are completed, build a pdb file from the final coordinates:

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
 
31. Characterise the pdb file (changepdb):

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.

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

1. Run the MD simulation of the pure solvent.
2. Build a pdb file from the last mdrest file with changeparm command p
3. Insert TER with changepdb command ter
   
4. Build an off file from the pdb file.
5. Change the box size back to that in the mdrest file (last line) in the off file
6. 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