Informational and File Sharing Server

Modeling of Macromolecular Systems

Tutorial 2: Including ligands in protein simulation

Molecule parametrization

Parametrization will be demonstrated on formate anione (HCOO-) that is a very simple molecule. There is General Amber Force Field (GAFF) that has parameters for almost all the organic molecules and it's compatible with ff99SB macromolecular force field. AmberTools package contains antechamber and parmchk programs designed to be used with GAFF and very helpful for new molecule parametrization. So antechamber identifies bonds and assignes atom types and parmchk finds for missing parameters. In this tutorial we will follow the classical parametrization ptocedure so molecular geometry will be optimized by PC Gamess at the MP2/6-31G* level of theory, molecular electrostatic potential will be calculated by PC Gamess at the HF/6-31G* level of theory, charges will be fitted by RESP that ships with AmberTools package.

At the beginning you have to get file formate.pdb of formate coordinates, for example you could draw formate in ACD/ChemSketch. Then you could use Ante_R.E.D. and R.E.D scripts to facilitate geometry optimization and charge derivation.
Set R.E.D. script variables as follows: $QMSOFT = "PC-GAMESS", $OPT_Calc = "Off", $CHR_TYP="RESP-A1".
Run Ante_R.E.D. to create formate-pcg.inp input file for geometry optimization and formate-out.p2n file used in charge derivation:

perl Ante_RED.pl formate.pdb > Ante_RED.log

Set ICHARG=-1, MPLEVL=2 and METHOD=NUMERIC in formate-pcg.inp and run geometry optimization:

pcgamess -i formate-pcg.inp -o Mol_red1.log

Rename formate-out.p2n to Mol_red1.p2n.
Edit REMARK lines in Mol_red1.p2n:

REMARK TITLE Formate
REMARK CHARGE-VALUE -1
REMARK MULTIPLICITY-VALUE 1

Give equivalent names to oxygens in third column of Mol_red1.p2n to ensure equal charges on these equal atoms:

ATOM      1  C1  FMT     1     18.985 -18.115   0.023                     C
ATOM      2  H1  FMT     1     18.887 -17.159   0.581                     H
ATOM      3  O2  FMT     1     18.003 -18.782  -0.210                     O1
ATOM      4  O2  FMT     1     20.208 -18.535  -0.398                     O2

Run R.E.D. in the folder that contains Mol_red1.p2n and Mol_red1.log to calculate charges:

perl RED.pl > RED.log

Coordinates and charges will be written in file Mol_m1-o1.mol2.
Now run antechamber to assign atom types:

antechamber -i Mol_m1-o1.mol2 -fi mol2 -o formate.mol2 -fo mol2

Coordinates, charges and atom types will be written in file formate.mol2.
Then run parmchk to find all force field parameters:

parmchk -i formate.mol2 -f mol2 -o formate.frcmod -a Y

In our case parmchk couldn't find parameters describing C-H bond and O-C-H angle of formate in GAFF. But it have proposed possible ones:

BOND
c -h5  357.50   1.078       same as c2-h5

ANGLE
o -c -h5   55.400     120.000   same as hc-c -o

After analysing GAFF atom type defination you could conclude that proposed parameters are appropriate. Now you should replace equilibrium bond lengths and angle values in formate.frcmod by ones obtained from geometry optimization.

We have assigned GAFF atom types to our molecule so it's necessary to load GAFF in tleap while preparing system with formate. To avoid it you could define new atom types for formate in formate.frcmod (don`t forget also to edit atom types in formate.mol2 after that). Notice that new atom type namings should differ from ff99SB ones.

Here is the final formate.frcmod:

AMBER FF94 and FF99 Formate Params
MASS
CG 12.01
HB 1.008
OA 16.00

BOND
CG-HB 357.5 1.140
CG-OA 648.0 1.260

ANGL
HB-CG-OA 55.4 114.52
OA-CG-OA 79.1 130.97

DIHE

IMPROPER
HB-OA-CG-OA      1.1    180.0   2.

NONB
  CG 1.9080 0.086
  HB 1.3590 0.015
  OA 1.6612 0.210

And here is essential part of formate.mol2:

@ATOM
      1 C           0.0000   -0.3460   -0.0000 CG        1 FMT      0.784100
      2 H          -0.0000   -1.4890    0.0000 HB        1 FMT     -0.170900
      3 O1         -1.1480    0.1770    0.0000 OA        1 FMT     -0.806600
      4 O2          1.1480    0.1770    0.0000 OA        1 FMT     -0.806600
@BOND
     1    1    2 1
     2    1    3 1
     3    1    4 1

Finally you should create formate.lib file:

tleap -s -f leaprc.ff99SB
> FMT = loadmol2 formate.mol2
> saveoff FMT formate.lib
> quit

It is worth mentioning that atomic charges are paramount parameters in Amber force field. In this tutorial only very simple protocol of charge derivation was considered in laconic manner. To learn more about charge development read through detailed tutorials on q4md-forcefieldtools.org.

By Dmitry Nilov, last modified: 14.03.11