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

  1. 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!
  2. 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.):

  1. antechamber -i lig.pdb -fi pdb -o lig.in -fo prepi -nc 1 -c bcc -rn LIG -at gaff
  2. \rm ANTECH* ATOMTYPE.INF PREP.INF
  3. Note: that you often need to symmetrise some charges, e.g. of H atoms bound to the same C.
  4. 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)
    ridft -resp -proper > logm
    
    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.
  • 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:

    1. Select program
    2. Press enter (new)
    3. Name of turbomole output (logm)
    4. Do not use Boltzman weights (press enter)
    5. Press enter (special requirements if # grid pints > 99999)
    6. Press enter (output: resp.pot)
    7. Press enter (resp.in / resp.in1)
    8. 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.

  1. Optimize the metal site with TPSS/def2-SV(P) and calculate the hessian (see Turbomole )
  2. Make atmtyp file (see template).
  3. 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
  1. Go to the folder where the hessian is calculated
  2. Run describe describe
  3. Select t (c/t=Turbomole)
  4. Press enter (i.e. use the default option "control")
  5. Press enter (for the output)
  6. Select temperature (press enter for default)
  7. Select scaling (with DFT usually close to 1)
  8. Force constant analysis (press enter)
  9. Select scaling (with DFT usually close to 1)
  10. If there are any unexpected bond, remove them (check what bond it is and that it is wise to remove)
  11. 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).
  12. Generate amber file (select "a")
  13. Enter the name of the atmtyp file (atmtyp). See template above.
  14. Give the residue a unique name
  15. Write logfile (press enter)
  16. Do not ignore sp3 dihidrals (press "n" or enter)
  17. Use ideal angles (press "i" or enter)
  18. Allow period=6 torsions (press "y" or enter)
  19. Select central atom (write number of the atom in the list printed to screen; typically this is the metal atom)
  20. Do not add atoms (press "n" or enter)
  21. Select new atom and repeat or STOP (press enter).
  22. Use default cut-off (press enter)
  23. 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

results matching ""

    No results matching ""