In general, your input file to LEaP will begin with several commands to source the relevant leaprc files. For
example the following preamble would allow you to include proteins, DNA, lipids, general components,
water, and atomic ions like Na+ or Cl-, using the current recommended force fields:
source leaprc.protein.ff14SB
source leaprc.DNA.OL15
source leaprc.lipid17
source leaprc.water.tip3p
source leaprc.gaff2
End-groups (OXT and AMT): If one wishes to have the end amino group not charged, then delete the lines in addPdbResMap. Copy file leaprc.ff03 to the local directory;
cp $AMBERHOME/dat/leap/cmd/leaprc.ff03 .Prep.in files: For special residues, a prep.in file must be constructed for the atom-types. A pdb file with the residue can be cut out from the pdbfile. The prep file can be constructed by (amber)
antechamber -i hic.pdb -fi pdb -o hic.in -fo prepi -nc 0 -c bcc -rn HIC -at amber # AMBER force field
antechamber -i hic.pdb -fi pdb -o hic.in -fo prepi -nc 0 -c bcc -rn HIC -at gaff # GAFF force fieldMissing parameters: In some cases, neither the regular amber protein force field nor the GAFF force field have parameters. In this case, the parameters must be constructed.
It is a good idea to check the resulting energies with the program changeparm
changeparm
Energies above 1000 kcal/mol should be checked.
source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
loadAmberPrep hpo.in
loadAmberParams o2.dat
loadAmberParams frcmod.ionsff99_tip3p ###parameters for counter ions ###
x=loadpdb 1icf4
bond x.156.SG x.206.SG
solvatecap x TIP3PBOX {0.0 0.0 0.0} 33
addions x Na+/Cl- 0 ###neutralize the system###
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:
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
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