Molecular mechanics
Table of Contents
AMBER
Construct a solute-solvent system
The procedure below first goes through making an MD for an organic solute, using cysteine (cysSH.pdb) as example in zwitterionic form (this procedure is the same as for any simple solute). More advanced solute-solvent systems (e.g. a solvated transition metal) follow a similar strategy with some important differences (noted below when required). The procedure is written with inspiration from an AMBER tutorial and the procedures described by Ulf Ryde (look for the sander.in files).
Two important concepts are explained below: periodic boundary conditions and non-bonded (electrostatic and van-der-Waals) interactions:
Periodic boundary conditions: Though some solvent molecules will be at the boundary between solute and solvent and others will be within the bulk of the solvent, a large number will be at the edge of the solvent and the surrounding vacuum. This is obviously not a realistic picture of a bulk fluid. In order to prevent the outer solvent molecules from boiling off into space, and to allow a relatively small number of solvent molecules to reproduce the properties of the bulk, periodic boundary conditions are employed. In this method the particles being simulated are enclosed in a box which is then replicated in all three dimensions to give a periodic array.
Electrostatics and van-der-Waals interactions: By employing a technique known as the Ewald sum, or its more modern equivalent the Particle Mesh Ewald (PME) method, it is possible to obtain the infinite electrostatics, even though we use a cutoff when simulating periodic boundaries in sander. However, the van-der-Waals interactions are still needed, which means we cannot make the cut-off too small. For production calculations, the ideal range is between 8 and 10 angstroms. A value of 8 is generally considered reasonable. Running with a larger cut-off will not cause harm, but it will make the calculation more expensive. One should never reduce the cut-off below 8 angstroms for periodic boundary (PME) calculations.
The steps in running the MD is provided below
1. Setup of simple system
- Obtain a pdb and arrange it. As an example, we use pdb called "cysSH.pdb" below, which contains a cysteine amino acid in zwitter ionic form. We name the cysteine residue and the residue CYZ (residue name). The pdb can be found at pdb.
- Use antechamber to construct an parameter file. These are called ".in" in the old format and ".frcmod" in the new format.
- Note: antechamber can give issues for more advanced systems - a more safe procedure is given under the advanced setup section. In case of problems, do this before the tleap step below.
antechamber -i cysSH.pdb -fi pdb -o cyz.in -fo prepi -nc 0 -c bcc -rn CYZ -at amber # AMBER force field (alternative -c resp)
antechamber -i cysSH.pdb -fi pdb -o cyz.in -fo prepi -nc 0 -c bcc -rn CYZ -at gaff # GAFF force field (alternative -c resp)
NOTE: "nc" indicates the total charge (should be changed for charged residues).
Note: antechamber can give issues for more advanced systems - a more safe procedure is given under the advanced setup section. In case of problems, do this before the tleap.
- Run
parmchk2 -f prepi -i cyz.in -o cyz.frcmod
to check for missing parametes. - Clean up the directory
rm ANTECH* ATOMTYPE.INF PREP.INF
.
pdb of cystein zwitterion
REMARK 4 COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
REMARK 888
REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
TITLE Displayed atoms
MODEL 1
ATOM 1 N CYZ A 1 -0.623 -1.017 0.831 1.00 0.00 N1+
ATOM 2 H CYZ A 1 -0.671 -0.583 1.742 1.00 0.00 H
ATOM 3 H1 CYZ A 1 -1.519 -1.441 0.638 1.00 0.00 H
ATOM 4 H2 CYZ A 1 0.060 -1.756 0.921 1.00 0.00 H
ATOM 5 CA CYZ A 1 -0.374 0.005 -0.185 1.00 0.00 C
ATOM 6 HA CYZ A 1 -1.235 0.074 -0.850 1.00 0.00 H
ATOM 7 CB CYZ A 1 -0.396 1.407 0.451 1.00 0.00 C
ATOM 8 HB2 CYZ A 1 -1.371 1.585 0.903 1.00 0.00 H
ATOM 9 HB3 CYZ A 1 0.377 1.471 1.217 1.00 0.00 H
ATOM 10 SG CYZ A 1 -0.094 2.738 -0.748 1.00 0.00 S
ATOM 11 HG CYZ A 1 -0.129 3.918 -0.123 1.00 0.00 H
ATOM 12 C CYZ A 1 0.941 -0.294 -0.923 1.00 0.00 C
ATOM 13 O CYZ A 1 1.625 -1.275 -0.637 1.00 0.00 O
ATOM 14 OXT CYZ A 1 1.278 0.597 -1.899 1.00 0.00 O1-
ENDMDL
END
2. Generation of AMBER coordinate and topology files with tleap
- Run tleap
tleap -s -f leap.in
with leap.in file is given below
leap files (example)
source leaprc.gaff # replace with "source leaprc.ff14SB" if amber was chosen in the [setup](#1-setup-of-simple-system) above
loadAmberPrep cyz.in
x=loadpdb cysSH.pdb
#
solvateOct x TIP3PBOX 15
saveamberparm x prmtop prmcrd
#
quit
Example with path to amber (aurora)
source /sw/pkg/amber/Amber16_Ambertools17/amber16/dat/leap/cmd/leaprc.water.tip3p
source /sw/pkg/amber/Amber16_Ambertools17/amber16/dat/leap/cmd/leaprc.gaff2
loadAmberPrep h2s.prepi
loadAmberParams /sw/pkg/amber/Amber16_Ambertools17/amber16/dat/leap/parm/frcmod.tip3p
loadOff /sw/pkg/amber/Amber16_Ambertools17/amber16/leap/lib/tip3pbox.off
HOH = TIP3P
#
x = loadpdb h2s.pdb
solvateOct x TIP3PBOX 15
saveamberparm x prmtop prmcrd
savepdb x prm.pdb
#
quit
pdb of cystein zwitterion
REMARK 4 COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
REMARK 888
REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
TITLE Displayed atoms
MODEL 1
ATOM 1 N CYZ A 1 -0.623 -1.017 0.831 1.00 0.00 N1+
ATOM 2 H CYZ A 1 -0.671 -0.583 1.742 1.00 0.00 H
ATOM 3 H1 CYZ A 1 -1.519 -1.441 0.638 1.00 0.00 H
ATOM 4 H2 CYZ A 1 0.060 -1.756 0.921 1.00 0.00 H
ATOM 5 CA CYZ A 1 -0.374 0.005 -0.185 1.00 0.00 C
ATOM 6 HA CYZ A 1 -1.235 0.074 -0.850 1.00 0.00 H
ATOM 7 CB CYZ A 1 -0.396 1.407 0.451 1.00 0.00 C
ATOM 8 HB2 CYZ A 1 -1.371 1.585 0.903 1.00 0.00 H
ATOM 9 HB3 CYZ A 1 0.377 1.471 1.217 1.00 0.00 H
ATOM 10 SG CYZ A 1 -0.094 2.738 -0.748 1.00 0.00 S
ATOM 11 HG CYZ A 1 -0.129 3.918 -0.123 1.00 0.00 H
ATOM 12 C CYZ A 1 0.941 -0.294 -0.923 1.00 0.00 C
ATOM 13 O CYZ A 1 1.625 -1.275 -0.637 1.00 0.00 O
ATOM 14 OXT CYZ A 1 1.278 0.597 -1.899 1.00 0.00 O1-
ENDMDL
END
3. Setup of advanced system
For simple systems (if the first step) with antechamber worked without problems) this can be skipped we can go directly to running the MD. Otherwise, we need to construct Non-bonded (and perhaps also bonded) parameters ourselves. We do this by constructing prepi and frcmod files (or prep.in and ".dat" files in the old format).
- Non-bonded parameters (electrostatics and van-der Waal): a guide for construction of resp charges (electrostatics) is found under RESP. This guide also has some hints for vDW parameters for non-stndard atoms. Note that antechamber can perhaps be used to generate the prepi (or prep.in) file, but sometimes it is not possible - if not, start with the bonded parameters since the "in" file can also be generated in this step. Alternatively, it must be constructed manually!
- Bonded parameters: can be constructed the approach by Seminario or Per-Orla Norby or PON
4. Running an equilibration and production MD
On aurora (cpu version) submit the run script below
#!/bin/bash
#SBATCH -N 1
#SBATCH -n 20
#SBATCH --exclusive
#SBATCH -t 130:00:00
ml purge
ml iomkl/2017b Amber/16.12-AT-17.08
cd $SLURM_SUBMIT_DIR
mpirun -bind-to core -np 20 pmemd.MPI -O -i sander.in1 -o sander.out1 -p prmtop -c prmcrd -r mdrest1 -ref prmcrd
mpirun -bind-to core -np 20 pmemd.MPI -O -i sander.in2 -o sander.out2 -p prmtop -c mdrest1 -r mdrest2 -ref prmcrd
mpirun -bind-to core -np 20 pmemd.MPI -O -i sander.in3 -o sander.out3 -p prmtop -c mdrest2 -r mdrest3 -ref prmcrd
mpirun -bind-to core -np 20 pmemd.MPI -O -i sander.in4 -o sander.out4 -p prmtop -c mdrest3 -r mdrest4 -ref prmcrd
mpirun -bind-to core -np 20 pmemd.MPI -O -i sander.in5 -o sander.out5 -p prmtop -c mdrest4 -r mdrest5 -x mdcrd5.nc
The individual input files are with a few comments below:
sander.in1
Minimization (SHAKE off)
&cntrl
irest=0,ntx=1,
imin=1,maxcyc=500,drms=0.0001,ntmin=2,
ntc=2,ntf=1,
cut=8.0,tishake=1,
ntpr=100,ntwx=0,ntwv=0,ntwe=0,
&end
sander.in2
20 fs MD at constant pressure with SHAKE for H-atoms
&cntrl
irest=0,ntx=1,
nstlim=10000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=1,
cut=8.0,tishake=1,
ntb=2, ! ntb=0 periodic boundaries. ntb=1 constant volume. ntb=2 constant pressure.
ntp=1,pres0=1.0,taup=1.0,
ntpr=5000,ntwx=0,ntwv=0,ntwe=0,
&end
sander.in3
Simulation 1 ns with SHAKE for H-atoms
&cntrl
irest=1,ntx=5, ! ntx=5 coord. vel. read from the prest file. ntx=1 only coord.
nstlim=500000,dt=0.002,
temp0=300.0,ntt=3,gamma_ln=2.0,
ntc=2,ntf=1,
cut=8.0,tishake=1,
ntb=2,
ntp=1,pres0=1.0,taup=1.0,
ntpr=5000,ntwx=5000,ntwv=0,ntwe=0,
&end
sander.in4
Equilibration, 1 ns NPT (constant pressure)
&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,tishake=1,
ntpr=500,ntwx=0,ntwv=0,ntwe=0,iwrap=0,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ntr=0
&end
sander.in5
Production 20 ns. Constant pressure
&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,tishake=1,
ntpr=500,ntwx=100,ntwv=0,ntwe=0,iwrap=0,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ntr=0
&end
3. Writing out pdb files (with cpptraj)
Run cpptraj -i cpptraj.in
An example file of cpptraj is given below
parm prmtop # topology file
trajin mdcrd5.nc 1 last 500 # input -see syntax below
autoimage # centering for Periodic Boundary Conditions
mask (:1<:15.0) maskpdb mask.pdb # select all residues within 15 Å of residue 1
#trajout out_frame_autoimage.pdb pdb multi # output -see syntax below
run
note: maestro cannot read the mask.pdb file
Convert with VMD
vmd -e input
The input is given before
mol load pdb "mask.pdb"
set a [atomselect top "resid 1 to 3155"]
atomselect0 writepdb "mask_conv.pdb"
exit
A script to do this over many files should be constructed...
syntax for trajin
trajin <filename> {[<start> [<stop> | last] [<offset>]]} | lastframe
<filename>: Trajectory file to read in.
[<start>]: Frame to begin reading at (default 1)
[<stop> | last]: Frame to stop reading at; if not specified or 'last' specified, end of trajectory
[<offset>]: Offset for reading in trajectory frames (default 1)
[lastframe]: Select only the final frame of the trajectory
syntax for trajout
trajout <filename> [<format>] [append] [nobox] [novelocity]
[notemperature] [notime] [noforce] [noreplicadim]
[parm <parmfile> | parmindex <#>] [onlyframes <range>] [title <title>]
[onlymembers <memberlist>]
[start <start>] [stop <stop>] [offset <offset>]
[ <Format Options> ]
<filename>: Trajectory file to write to.
[<format>]: Keyword specifying output format
[append]: If <filename> exists, frames will be appended to <filename>
[nobox]: Do not write box coordinates to trajectory
[novelocity]: Do not write velocities to trajectory
[notemperature]: Do not write temperature to trajectory
[notime]: Do not write time to trajectory
[noreplicadim]: Do not write replica dimensions to trajectory
[parm <parmfile>]: Topology filename/tag to associate with trajectory (default first topology)
[parmindex <#>]: Index of Topology to associate with trajectory (default 0, first topology)
[onlyframes <range>]: Write only the specified input frames to <filename>
[title <title>]: Output trajectory title
[onlymembers <memberlist>]: Ensemble processing only; only write from specified members (starting from 0)
[start <start>]: Begin output at frame <start> (1 by default)
[stop <stop>]: End output at frame <stop> (last frame by default)
[offset <offset>]: Skip <offset> frames between each output (1 by default)
Restraints
From the AMBER 16 manual, there are several ways to restrain atoms in an MD. Both use a special syntax for the mask explained below.
restraintmask: String that specifies the restrained atoms when ntr=1.
belymask: String that specifies the moving atoms when ibelly=1.
A couple tips:
1) Use residue ranges when possible. For example,
:1-10 # is equivalent to :1,2,3,4,5,6,7,8,9,10
or combinations of ranges :1-5,6,8-10
2) Use relational operators (!, &, |) . Specifically the 'not' operator (!) is quite helpful if you want to select all residues except certain ones
!:1-5,6,8-10 # will match all residues except numbers 1, 2, 3, 4, 5, 6, 8, 9, and 10.
3) Use distance mask selections.
(.20 < :10.0) # select all residues within 10 angstroms of atom number 20.
Another example:
(:20 < .10.0) # This will select all aitoms within 10 angstroms of residue number 20.
Any and all of these can be combined. For instance,
(.20 < :10) & :WAT
# will select all atoms that are in water residues within 10 angstroms of atom number 20.
4) The command ambmask can be used to specify the number of atoms affected by a given mask
ambmask -p prmtop -c inpcrd -prnlev [0-3] -out [short| pdb| amber] -find [maskstr]
5) Include the belly or restrainmask in the sander inpit file (sander.inX).
ibelly=1,bellymask=":WAT | @H="
or
ntr = 1, ! Turn on positional restraints
restraint_wt = 10, ! 10 kcal/mol/A**2 restraint force constant
restraintmask = '@CA,C,O,N' ! Restraints on the backbone atoms only
RESP charges
Try first the automatic precedure (for ligand LIG with pdb file lig.pdb). Arrange lig.pdb in a way that is easy to read (e.g. H atoms after the C atom they are bonded to etc.):
antechamber -i lig.pdb -fi pdb -o lig.in -fo prepi -nc 1 -c bcc -rn LIG -at gaff
\rm ANTECH* ATOMTYPE.INF PREP.INF
- Note: that you often need to symmetrise some charges, e.g. of H atoms bound to the same C.
- Note: We have sometime seen issues due to the pdb format (e.g. for pdb automatically generated with openBabel -here all connectivity (including the line "MASTER") should be deleted.
Charges with Gaussian
Procedure taken from Ulf Ryde's page
Optimise the molecule with QM methods, typically B3LYP/6-31G* or AM1 for a large molecule.
Calculate the electrostatic potential. For the Amber-99SB force field, this should be done at the HF/6-31G level (B3LYP/6-31G if you have metal sites). Note: Other force fields require other levels of theory.
Example of a Gaussian input file:
#P B3LYP/6-31g* Opt
Title, e.g. B3LYP, 6-31G*, Date
-1 1
c 33.34700000000000 7.55100000000000 14.58000000000000
...
h 36.15300000000000 11.00700000000000 16.90000000000000
--Link1--
%Chk=checkpoint.chk
#P B3LYP/6-31G* SCF=Tight Geom=AllCheck Guess=Read
Pop=MK IOp(6/33=2,6/41=10,6/42=17)
- Run antechamber on the output file of this calculation:
antechamber -fi gout -fo prepi -c resp -i file.out -o res.in -rn RES -at amber -pf y
where, you should insert the name of the output file from the previous point (gout), and the residue name "res", (three letters) twice, first in lower-case and second in capital letters.Use -at amber, if want to use only standard Amber atom types. The default is to get GAFF atom types, which is normally better. The -pf y you will delete intermediate files. Alternatively, you can delete them by: \rm ANTECH* ATOMTYPE.INF PREP.INF
Note: Often Amber gives you strange atom names. This should be corrected (a procedure is found at Ulf Ryde's page)
Note: you often need to symmetrise some charges, e.g. of H atoms bound to the same C.
If antechamber fails, go to the procedure below:
Charges with Turbomole
This procedure is copied from protein setup pages
- Run a geometry optimization on the non-standard residues [typically B3LYP/def2-SV(P) or TPSS/def2-SV(P)]. Make sure that the structure is not too different from the crystal structure or other trustworthy experimental structures (if there are large deviations, consider using a single point calculation or only optimized hydrogens to generate charges)
- With Turbomole (must be version > 7.0) insert "$esp_fit resp in the "control" file 'and run turbomole (Note: It can be convenient to ensure that the order of atoms is the same as in the pdb file)
The use of "-proper" indicates that no wave function optimization is required. Look for "GRID" in logm to check whether the Electrostatic Potential (which will be used to fit the resp charge) is generated.ridft -resp -proper > logm
Use now the script "changepot" (generates resp.pot, resp.in1, resp.in and qin). Note: changepot might give an error due to a missing line in logm. Insert this line with:
awk '/charge for group # 1/{x++} x==2{sub(/charge for group # 1/,"writing ESP data for Amber to file esp.dat \n&")}1' logm > logm_2 && mv logm_2 logm
changepot:
- Select program
- Press enter (new)
- Name of turbomole output (logm)
- Do not use Boltzman weights (press enter)
- Press enter (special requirements if # grid pints > 99999)
- Press enter (output: resp.pot)
- Press enter (resp.in / resp.in1)
- Press enter (qin)
The files resp.in and resp.in1 are input files for the fitting of resp charges. qin are initial charges. Fitting of resp charges are usually done in two steps, first a fit involving all atoms (check this) with loose convergence. This is done with resp.in. In a next step, all heavy atoms are fixed while hydrogens in methylene (CH2) and methyl (CH3) groups are constrained to be identical and the fit is done with more tight convergence (see below).
Check if there are any polar hydrogens bound to the same atom (e.g. water). These should be constrained to have the same charge (this is done in resp.in, see below for the format):
Format of resp.in files:
&cntrl ihfree=0, qwt=0.0005, iqopt=2 # qwt = optimization parameter &end 1.0 Mol1 2 63 # charge and tot. # of atoms 7 0 # Atom nr. (here 7 = N) and command (-1 = freeze, 0 = optimized freely, n = fix to the same as atom #n e.g. 43 means fix to same value as atom number 43) 6 0 6 0
Run resp (part of the AMBER suite)
resp -O -i resp.in -o resp.out -q qin -e resp.pot
Open resp.in1: Fix all charges except carbon atoms with more than one hydrogen (usually CH3 or CH2 groups). Fix the the hydrogens to have the same charge). qwt should be 0.001. Save qout in qin1:
cp qout qin1
and run
resp -O -i resp.in1 -o resp.out1 -q qin1 -e resp.pot
The final charges are now in resp.out1
van-der Waal parameters
For non-standard atoms, the CNS (crystallography) force field offers a good choice.
write here how they are included (under NONB or NONBOND)
Seminario parameterization
Use the method by Seminario (JM Seminario, Int. J. Quant. Chem. Quant. Chem. Symp. 30, 1996, 59-65). The in-house program describe (by Ulf Ryde) can extract the required force-constants. Normally, only force constant /angles and dihedrals are required for the bonds directly bond to the metal. For other atoms in the site we use standard force fields.
- Optimize the metal site with TPSS/def2-SV(P) and calculate the hessian (see Turbomole )
- Make atmtyp file (see template).
- Use describe to set up the force field (from folder where the hessian is calculated). See separate guide below
atmtyp template
LPMO active site [Cu(H2O)(HIC)His]2+ (edh 26/11 2018)
CUL ! name for the ligand
N n1 ! Bond to Cu. Need special type. Format: PDB NAME ATOM TYPE
H2 H2 ! AMBER atmtyp
H3 H3 ! AMBER atmtyp
CA CX !...(same for the rest)...
HA H1
CB CT
HB2 HC
HB3 HC
CG CC
ND1 n2 ! Bond to Cu
CE1 CR !
HE1 H5 !...(same for the rest)...
NE2 N*
CZ CT
HZ1 H1
HZ2 H1
HZ3 H1
CD2 CW
HD2 H4
C C
O O
HA HC ! introduced as cap
CB CT ! AMBER atmtyp
HB2 HC !...(same for the rest)...
HB3 HC
CG CC
ND1 NA
HD1 H
CE1 CR
HE1 H5
NE2 n3
CD2 CV
HD2 H4
CU cu
O ol ! bond to Cu (H2O)
H1 hl ! H2O
H2 hl ! H2O
H H5 ! introduced as cap
- Go to the folder where the hessian is calculated
- Run describe
describe
- Select t (c/t=Turbomole)
- Press enter (i.e. use the default option "control")
- Press enter (for the output)
- Select temperature (press enter for default)
- Select scaling (with DFT usually close to 1)
- Force constant analysis (press enter)
- Select scaling (with DFT usually close to 1)
- If there are any unexpected bond, remove them (check what bond it is and that it is wise to remove)
- Add extra pairs (if they are all there, no need to add any). Otherwise write "XX YY" to add the bond suggested by "Missing bond XX YY" (XX and YY are the atom numbers).
- Generate amber file (select "a")
- Enter the name of the atmtyp file (atmtyp). See template above.
- Give the residue a unique name
- Write logfile (press enter)
- Do not ignore sp3 dihidrals (press "n" or enter)
- Use ideal angles (press "i" or enter)
- Allow period=6 torsions (press "y" or enter)
- Select central atom (write number of the atom in the list printed to screen; typically this is the metal atom)
- Do not add atoms (press "n" or enter)
- Select new atom and repeat or STOP (press enter).
- Use default cut-off (press enter)
- Use default cut-off for hydrogen bonds (press enter)
Parameters for metal site
Here I will describe how to use the MCPB.py framework in AMBER