Ammm:Ff gau grom prdrug

Introduction

When do we need to create parameters?

• System size is too large for QM and
• Molecule has not been parameterized for the force field being used or
• Molecular environment or conformation is substantially different

Overview of parameterization:

Quantum Mechanics (QM)

• Optimized (a.k.a. minimum energy) geometry search
• Single point calculation to fit the electrostatic potential
• Torsion energy profile to calculate torsion energy term

Validation of Molecular Mechanics (MM) model

• Optimized structure comparison (to check for errors)
  
• Comparison to experimental observables (to assess quality of model) such as:
  • density
  • radial distribution function
  • melting temperature
  • conformation, solvation energy, etc.

Gaussian Overview

Gaussian 09 Online Manual

Gaussian is a well-developed quantum mechanics (QM) software package with a wide range of features. The software is not free, however, so if you require a free software distribution GAMESS is one alternative.

Input

All input parameters, atoms, and coordinates are given in the COM File, which have 5 sections:

• "Link 0" Commands
• Route section (blank line terminated)
• Title section (blank line terminated)
• Molecule specification (blank line terminated)
• Optional additional statements

Example:

%NoSave
%chk=PYC1mer3.chk
%Mem=400MB
%Nproc=2

# HF/6-31G** Opt Pop=ChelpG

Title Section

0 2
C 1.166043 -0.148206 -0.552387
C 0.225157 -1.168482 -0.483464
N -1.033322 -0.561634 -0.411853
C -0.887040  0.829496 -0.436250
C 0.474401  1.094278 -0.523043
H 2.243332 -0.273241 -0.617421
H 0.375213 -2.246760 -0.480835
H -1.724118  1.524316 -0.392990
H 0.935360  2.077302 -0.562253
H -1.893240 -1.042210 -0.355061
Cl 0.083843  0.014677  1.015435

The first two numbers in the Molecule specification section (e.g. "0 2") are the net charge and spin multiplicity of the system. The Gaussian website has good descriptions of all the keywords and various options.

If you are converting your molecular coordinates from another filetype, Open Babel can convert between Gaussian and other filetypes like PDB. Practical considerations such as memory allocations, processors, disk space usage, scratch file directory, etc. will vary between system setups. It is worthwhile to test out a few different options on your system to find the optimal settings if you will be doing lots of QM calculations.

Output

Gaussian keeps a checkpoint file in the binary CHK format, and text output is held in a LOG file. You may visualize the CHK file by converting it to an FCHK ("formatted checkpoint") file using the "formchk" built-in commnad. GaussView and gOpenMol are two visualization packages among many that support fchk files.

Example LOG file snippits:

... skipping ...
Breneman (CHELPG) radii used.
Generate Potential Derived Charges using the Breneman model, NDens= 1.
Grid spacing= 0.300 Box extension= 2.800
NStep X,Y,Z=   34     35     26   Total possible points=       30940
Number of Points to Fit=    8561

**********************************************************************

           Electrostatic Properties Using The SCF Density

**********************************************************************

      Atomic Center    1 is at   1.166043 -0.148206 -0.552387
... skipping ...
      Atomic Center   10 is at  -1.893240 -1.042210 -0.355061
      Atomic Center   11 is at   0.083843  0.014677  1.015435
   8561 points will be used for fitting atomic charges
Fitting point charges to eletrostatic potential
Charges from ESP fit, RMS=   0.00279 RRMS=   0.21577:
Charge=   0.00000 Dipole=    -0.5792    -0.4886    -1.8898 Tot=     2.0360
             1
    1  C   -0.019083
    2  C   -0.078198
    3  N   -0.194869
    4  C   -0.099409
    5  C   -0.037277
    6  H    0.096809
    7  H    0.133066
    8  H    0.133803
    9  H    0.092851
   10  H    0.239509
   11  Cl  -0.267201
-----------------------------------------------------------------
 ... skipping ...

-----------------------------------------------------------------
1\1\GINC-BME-NOVA\SP\UHF\STO-3G\C4H5Cl1N1(2)\JOHNFONNER\16-Feb-2007\0\
\# POP=CHELPG\\PYC1mer\\0,2\C,0,0.008576,1.177221,-0.006596\C,0,1.1383
94,0.368298,0.005025\N,0,0.693173,-0.957942,0.010301\C,0,-0.705552,-0.
984838,0.00126\C,0,-1.13812,0.335747,-0.008954\H,0,-0.001718,2.263623,
-0.012784\H,0,2.189699,0.651029,0.010302\H,0,-1.29065,-1.902992,0.0019
22\H,0,-2.171113,0.672212,-0.017564\H,0,1.27731,-1.753163,0.017088\Cl,
0,0.,0.,1.5\\Version=IA32L-G03RevC.02\State=2-A\HF=-660.1236865\S2=0.7
9721\S2-1=0.\S2A=0.750729\RMSD=6.593e-04\Dipole=0.2119298,-0.1685267,-
0.7346684\PG=C01 [X(C4H5Cl1N1)]\\@


Some scientists claim that hydrogen, because it is so
plentiful, is the basic building block of the universe. 
I dispute that. I say that stupidity is far more
abundant than hydrogen, and THAT is the basic building 
block of the universe.
                  --Frank Zappa
Job cpu time:  0 days  0 hours  0 minutes  6.2 seconds.
File lengths (MBytes):  RWF=     11 Int=      0 D2E=      0 Chk=      7 Scr=      1
Normal termination of Gaussian 03 at Fri Feb 16 12:30:21 2007.

Gromacs Overview

When using custom parameters in Gromacs, there are two ways modify parameters:

• Modify the PRM file for the force field of choice and create a matching PDB file with atom coordinates

(This option is probably the best choice when you plan to run simulations starting from multiple, different starting structures.)

• Work within Gromacs formats and create a TOP file with a matching GRO file.

(when only changing a few atom properties or when using the exact same atom numbering across simulations, this may be the simplest solution)

PDF File format:

"ATOM or HETATM" [atom #] [atom name] [residue name] [molecule # (1 char)] [residue #] [coordinates] [occupancy] [temperature factor]

ATOM      1  HTM PYB 0   1      22.240   0.400  17.720  1.00  0.00     
ATOM      2  CC1 PT1 0   2      22.930   2.490  17.860  1.00  0.00     
ATOM      3  CN1 PT1 0   2      22.590   1.240  18.290  1.00  0.00     
ATOM      4  N   PT1 0   2      22.210   1.320  19.600  1.00  0.00     
ATOM      5  CN2 PT1 0   2      22.700   2.500  20.100  1.00  0.00     
ATOM      6  CC2 PT1 0   2      23.130   3.260  19.050  1.00  0.00     
ATOM      7  HC1 PT1 0   2      22.990   2.880  16.870  1.00  0.00     
ATOM      8  HC2 PT1 0   2      23.450   4.280  19.120  1.00  0.00     
ATOM      9  HN  PT1 0   2      21.690   0.680  20.150  1.00  0.00     
ATOM     10  CC1 PYC 0   3      21.590   2.630  22.400  1.00  0.00     
ATOM     11  CN1 PYC 0   3      22.580   2.870  21.500  1.00  0.00     
ATOM     12  N   PYC 0   3      23.580   3.560  22.140  1.00  0.00     
...skipping...
ATOM   6075  CC2 PYU ^ 865      73.223  72.090  20.960  1.00  0.00     
ATOM   6076  HC1 PYU ^ 865      71.423  70.840  20.930  1.00  0.00     
ATOM   6077  HC2 PYU ^ 865      73.093  72.910  20.270  1.00  0.00     
ATOM   6078  HN  PYU ^ 865      75.053  70.390  22.870  1.00  0.00     
ATOM   6079  HTM PYE ^ 866      75.283  72.580  21.670  1.00  0.00     
TER
ENDMDL

Residue names must match the Gromacs parameter files for the force field you are using (more on this below!). Atom names within each residue must be unique and must match the parameter files. Gromacs connects residues based on their order in the PDB file, so keep everything sequential. If any atoms or parameters are missing grompp should give an error message (or atleast a warning).

Gromacs Parameter Files

Gromacs stores all of its parameter files in a /top/ directory. This directory contains several files for each supported force field. For example, taking a look at the OPLSAA force field, we find the following files:

ffoplsaa-n.tdb
ffoplsaa.atp
ffoplsaa.ddb
ffoplsaa.hdb
ffoplsaa.itp
ffoplsaa.n2t
ffoplsaa.rtp
ffoplsaabon.itp
ffoplsaanb.itp

Each of the files hold different parameters, some of which overlap. It can be confusing what the best way is to add a new molecule or building block. When changing these files for your molecule, I recommend changing a few things, running grompp, and checking the topol.top file to see if your changes worked. The filetypes are the following:

• ATP - Atom Type Parameters - contains atom types and their mass
• ITP - Interaction Parameters - contains bonded parameters (*bon.itp) and non-bonded parameters (*nb.itp) for all commonly used molecules. New parameters may be added here or to an RTP file. The Gromacs manual recommends putting new molecules either in these ITP files or to create a new molecule.itp file.
• RTP - Residue Database - contains charge, charge groups, bonds, angles, and dihedrals for building blocks such as amino acids. The Gromacs manual recommends putting new building blocks or amino acids in this file.
• HDB - Hydrogen Database - tells gromacs where to add hydrogens for amino acids
• TDB - Terminal Database - used by Gromacs to put capping groups at the N- or C-terminal of a protein
• DDB - Dummy types Database - I have never used it, but it contains dummy atom specifications

To talk about some caveats associated with adding parameters, let's look at an example.

Example Parameters

The following is an excerpt of parameters I have added to ffoplsaa.rtp:

---

[ PYU ]
 [ atoms ]
; opls types are given in ffoplsaa.atp file
; name    number      charge    charge_group
 CC1     opls_544    -0.215    1
 CN1     opls_543     0.107    1
   N     opls_542    -0.386    1
 CN2     opls_543     0.107    1
 CC2     opls_544    -0.215    1
 HC1     opls_547     0.141    1
 HC2     opls_547     0.141    1
  HN     opls_545     0.320    1
 
[ bonds ]
; atom1   atom2    b0         kb
   N     HN       .09937     363171
   N     CN1      .13657     282001
   N     CN2      .13657     282001
 CN1     CC1      .13644     392459
 CN2     CC2      .13644     392459
 CC1     CC2      .14223     392459
 CC1     HC1      .10715     307105
 CC2     HC2      .10715     307105
-CN2     CN1      .14567     392459  ;bond to previous monomer unit

 [ angles ]
; atom1 atom2 atom3    th0        cth
 HC1   CC1   CN1      125.93     292.88
 HC1   CC1   CC2      126.54     292.88
 HC2   CC2   CN2      125.93     292.88
 HC2   CC2   CC1      126.54     292.88
 HN    N     CN1      124.42     317.98
 HN    N     CN2      124.42     317.98
 N     CN1  -CN2      121.34     527.18
 N     CN2  +CN1      121.34     527.18
 N     CN1   CC1      107.49     527.18
 N     CN2   CC2      107.49     527.18
 CN1   N     CN2      110.12     527.18
-CN2   CN1   CC1      131.17     527.18
 CN1  -CN2  -CC2      131.17     527.18
 CN1   CC1   CC2      107.44     527.18
 CN2   CC2   CC1      107.44     527.18

 [ dihedrals ]
; atom1 atom2 atom3 atom4
 HC1   CC1   CC2   HC2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 HC1   CC1   CC2   CN2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 CN1   CC1   CC2   HC2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 CN1   CC1   CC2   CN2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 

 HC1   CC1   CN1  -CN2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 HC1   CC1   CN1   N      dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 CC2   CC1   CN1  -CN2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 CC2   CC1   CN1   N      dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 

 HN    N     CN1  -CN2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0
 HN    N     CN1   CC1    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0
 CN2   N     CN1  -CN2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0
 CN2   N     CN1   CC1    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
      
-HN   -N    -CN2   CN1    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 HN    N     CN2   CC2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
-CN1  -N    -CN2   CN1    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 CN1   N     CN2   CC2    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
 
 HC2   CC2   CN2   N      dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
-HC2  -CC2  -CN2   CN1    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0
 CC1   CC2   CN2   N      dih_PYU_ring   ; 30.33  0  -30.33  0  0  0 
-CC1  -CC2  -CN2   CN1    dih_PYU_ring   ; 30.33  0  -30.33  0  0  0
  
 N     CN1  -CN2  -N      dih_PYU_backbone  
 N     CN1  -CN2  -CC2    dih_PYU_backbone_trans  
 CC1   CN1  -CN2  -N      dih_PYU_backbone_trans  
 CC1   CN1  -CN2  -CC2    dih_PYU_backbone  

 [ impropers ]
; atom1  atom2  atom3  atom4  q0  cq
; N     HN    CN1   CN2     improper_Z_CA_X_Y
; CN1   N     CC1  -CN2     improper_Z_CA_X_Y
;-CN2  -N    -CC2   CN1     improper_Z_CA_X_Y
; CC1   CN1   CC2   HC1     improper_Z_CA_X_Y
; CC2   CN2   CC1   HC2     improper_Z_CA_X_Y

; CN1   CN2   CC1   CC2     improper_Z_CA_X_Y

Notes on parameters:

• Gromacs uses different units that some other software packages

Parameter	Units
bond length	nm
bond force constant	kJ/mol/nm^2
angle & dihedral force constants	kJ/mol/rad^2

• Dihedral Angles are expressed as a Ryckaert-Belleman function rather than a Fourier series.

To convert from a Fourier series: $V_{0}+1/2\left[V_{1}(1+cos(x))+V_{2}(1-cos(2x))+V_{3}(1+cos(3x))+V_{4}(1-cos(4x))\right]$

To a RB function: $C_{0}-C_{1}cos\left(x\right)+C_{2}{cos(x)}^{2}-C_{3}{cos(x)}^{3}+C_{4}{cos(x)}^{4}$ Use the following equations:

$C0=V_{0}+{\frac {V_{1}}{2}}+V_{2}+{\frac {V_{3}}{2}}$
$C1={\frac {1}{2}}(3V_{3}-V_{1})$
$C2=4\left(V_{4}\right)-V_{2}$
$C3=-2\left(V_{3}\right)$
$C4=-4\left(V_{4}\right)$

• When connecting building blocks, use "-" to define the connection to the previous block and use "+" only when there is no other way. Be careful not to define the same parameter twice, but if you do, expect an error/warning message.
• The RTP file has a short character limit when defining dihedral angles, so it is easiest to #define the angle in the *bon.itp file and reference it in the rtp file. For example, the "dih_PYU_ring" dihedral parameter referenced above looks like this in the *bon.itp file.

#define dih_PYU_ring                        30.33400   0.00000 -30.33400   0.00000   0.00000   0.00000

Prdrug Server Overview

The Dundee PRODRG Server

PRODRG Frequently Asked Questions

Case Study

We decided to create parameters for polypyrrole (PPy), an amorphous, semiconducting polymer that is "doped" with anions to improve conductivity. Our goal is to simulate peptide interactions on the surface of a polypyrrole mesh. To streamline the explanation of the process we used, I will only talk about undoped PPy rather than doped PPy.

Initial Decision Making

• We estimate our final system size to be on the order of 100,000 atoms, so we chose to use a fixed charge model to lower computational cost.
• Since we planned on performing liquid simulations of proteins, we chose to integrate our parameters into the OPLSAA set (Optimized Parameters for Liquid Simulations) by Jorgensen.
• To fit into this set, we decided to use the HF/6-31G* QM basis set. This basis set overestimates charge in vacuum, but this overestimation more accurately estimates charge for condensed phase materials.
• Since PPy molecules have different chain lengths and doping arrangements, we need parameters for "building blocks".

QM Parameterization

Charges, Bonds, Angles, VdW

Atomic partial charges the equilibrium values for bonds and angles may all be taken from the minimum energy structure of a molecule. For a building block, however, one should be careful to avoid bias due to "edge effects", and charges must be adjusted to zero to preserve charge neutrality. We did this by simulating small PPy oligomers of 3, 5, and 7 rings in length and looked at the structure of the center ring. Below are the minimum structures of the trimer and pentamer:

Van der Waals terms do not vary significantly for a given atom type, so VdW parameters were borrowed from similar atom types with the OPLSAA force field. The same may be said for the force constants of bonds and angles.

Torsion Profile

The torsion parameter accounts for rigidity in the bond along the backbone, and should be calculated for pretty much any new building block that doesn't closely resemble an existing one. Since the torsion parameter is usually the last parameter calculated, it also acts as an error function, so caution must be taken not to hide a poorly parameterized system through a well parameterized torsion term. For my system, I found that by using four units instead of 2 when calculating torsion, it exposed the error from the other terms.

By rotating and restraining the backbone torsion angle in 15 degree increments, we optimized the geometry of the system in Gaussian with the HF/6-31G* basis set and performed a single point energy calculation using the more accurate MP2/6-31(2d,2p) basis set. The resulting energy profile looks like the following:

The torsion energy term is calculated by creating the same torsion profile for the molecular mechanics (MM) model without any torsion parameter. This is then subtracted from the QM torsion profile, and the difference is fit to a Fourier series.

This torsion term is then added back into the model, and final comparisons can be made between QM results and MM results.

MM Validation

Methods for validation will vary significantly according to the material properties and the experimental results that are available/feasible. For all cases that I can think of, the RMSD value for the minimum energy structure is a valid check. The more validations you have the better, and the sooner you can check the better. Density is a good indicator of the quality of a force field, but it can be difficult to obtain. Also, incorrect density values do not elucidate the part of your force field that is causing the error.

Minimum Energy Structure

PyMOL can align two molecules and calculate the RMSD.

RMS: ~0.145 Å

Density Calculations

Density should be calculated using a NPT ensemble at room (or body) temperature and at 1 atmosphere (~1 bar). Notice the temperature and pressure fluctuations below.

File:Final plot.dat.bmp