ParamIT: a Toolset for Molecular Mechanical Force Field Parameterization within the CHARMM/CGenFF Protocol

Nikolay A. Simakov and Troy Wymore

About

This is a compilation of our ParamIT toolkit presentation at the 2012 Biophysical Society Meeting. This presentation illustrates some of the major capabilities of the ParamIT. The current roadmap is to make a per request release by the end of May, 2012.

Abstract

We have developed several tools to aid the molecular mechanical force field (FF) parameterization of small molecules using the CHARMM General Force Field (CGenFF) approach(Vanommeslaeghe, Hatcher et al. 2010). The intermolecular interactions of relatively small organic molecules with biological macromolecules is critical to understanding many biological phenomena. Accurate force fields of these small organic molecules are important for long timescale molecular dynamics (MD) simulations. Given the massive variety of organic compounds, quite often the parameters for small organic molecule of interest are not present in existing molecular force fields or at least not very well developed and thus need to be derived through extensive quantum chemistry calculations. MacKerell’s group has developed a detailed protocol and FF for the parameterization of small molecules (CGenFF) that is compatible with CHARMM FF. Although, the CGenFF protocol is sufficiently detailed, the parameterization remains a labor intensive task due to large number of calculations needed to be performed in self-consistent manner. We developed tools that help the researchers in following ways: 1) automating the creation of multiple input files for quantum and molecular mechanics (MM) programs, 2) automating the output analysis and 3) substitute the use of full MM programs with a faster specialized one. The developed tools include: 1) generator of molecule-water complexes with graphical user interface (GUI), 2) semi-automatic frequency analysis using symbolic potential energy distribution matrix and comparison of optimized internal coordinates, 3) GUI for charge fitting with three modes: manual, Monte-Carlo sampling or brute force, and 4) GUI for dihedral terms fitting. The usage of these tools decreases the labor effort, lowers manual input errors and reduces the time needed for accurate MM parameterization efforts.

Molecular Mechanical Energy Function

Molecular mechanical energy function in CHARMM FF(MacKerell, Bashford et al. 1998) can be represented as a sum of non-bonded and bonded interactions. The non-bonded interaction energy is given by:

The bonded interaction energy has following form:

Implementation

To aid fast and yet accurate parameterization we are developing a complete toolkit which interfaces to Gaussian(Frisch 2009) and CHARMM program(Brooks, Brooks et al. 2009) for quantum chemical (QC) and molecular mechanical (MM) calculations respectively. ParamIt toolkit is a complete set of tools for MM parameterization and it features modular design with unified input/output enabling independent enhancements of its parts. The toolkit is mainly written in Python with computationally intensive routines written in C language. Several tools are implemented as graphical plug-ins to VMD(Humphrey, Dalke et al. 1996). NumPy, matplotlib and Python megawidgets libraries are used for graphical interface and visualization.

Complexes with Water Generator

Within the CHARMM approach the atomic partial charges are optimized against hydrogen bond energies and distances of QC optimized molecule - water complexes. Generator of molecule-water complexes provides a simple graphical interface to set-up initial conformation and generates the input for Gaussian. The graphical interface simplifies the generation of molecule-water complexes which allows the generation of a higher number of complexes. This can lead eventually to more accurate parameters or at least will inform the user about poorly parameterized groups in the molecule.

Charge fitting plug-in to VMD

Charge fitting plug-in to VMD provides a graphical interface for the optimization of atomic charge. The interface allows the manipulation of the molecular charge distribution with immediate feedback on the agreement of MM calculated dipole moments, hydrogen-bond energies and distances with previously calculated QC results. The interface features three modes to optimize charge distribution. 1) Manual – the changes are changed by the user, the plug-in recalculates the dipole moments and optimize the complexes with water, this allows the manual interactive optimization of the molecule partial charges. 2) Monte-Carlo (MC) sampling using a scoring function as an energy function of the system, the plug-in will perform a MC sampling in the space of charges. The charge set with lowest score will be loaded to the plug-in display. And 3) brute force (BF) sampling – a straight scan throughout all possible charge combinations within a specified charge range and step is performed, the results are sorted according to scoring function. Then the samples from the scan can be loaded to display for further analysis and optimization. These three modes can be used in conjunction with each other. For example one might perform BF scan first with further refinement with downhill MC (i.e. with zero temperature).

Vibrational Analysis

Vibration modes are used to obtain force constants for bonds, angles and partially dihedral terms. qcfqan.py, a MOLVIB module wrapper, simplifies the input of non-redundant internal coordinates from molecule internal coordinates (i.e. Umat) for symbolic potential energy distribution (PED) analysis. The tools also provide a capability to create Umat from the molecule’s Z-Matrix. This can serve as a good initial guess or used as is. Below is an example of input for vibrational analysis:

#Location of CHARMM binary
charmm="/usr/local/charmm/c34b2/exec/em64t/charmm"
#Gaussian output with frequency calculations
fqlog="fq.log"
#Generated CHARMM input file
charmmin="qcmolvib.inp"
#Where to store the output from CHARMM
charmmout="qcmolvib.out"
#INTERNAL COORDINATES DEFENITION (keep this label)
#func gic(T, Ai, Aj, [Ak], [Al]) 
#   gic (get internal coordinate) - return the identifier for the molecule’s
redundant IC. 
#   T – Type of IC, can be B(bond), A(angle), D(dihedral) or M(improper
dihedrals, Wilson wag)
#func umatApp([FC0,IC0,FC1,IC1,...,name]) 
#    umatApp - append a row to Umat. The entry is the combination of 
#Umat Section:
#C01
umatApp([1.0,gic(B, H06, C01), 1.0,gic(B, H07, C01), 1.0,gic(B, H08, C01),"ssC01H"])
umatApp([-1.0,gic(B, H06, C01), 2.0,gic(B, H07, C01), -1.0,gic(B, H08, C01),"saC01H"])
umatApp([-1.0,gic(B, H06, C01), 0.0,gic(B, H07, C01), 1.0,gic(B, H08, C01),"saC01H’"])
umatApp([ 1.0,gic(A, H06, C01, H07), 1.0,gic(A, H07, C01, H08), 1.0,gic(A, H08,C01, H06),
         -1.0,gic(A, H06, C01,C02),-1.0,gic(A, H07, C01, C02),
         -1.0,gic(A, H08, C01,C02), "dsC01H"])   
umatApp([ 1.0,gic(A, H06, C01, H07), 1.0,gic(A, H07, C01, H08),-2.0,gic(A, H08,C01, H06), "daC01H"])
umatApp([ 1.0,gic(A, H06, C01, H07),-1.0,gic(A, H07, C01, H08), 0.0,gic(A, H08,C01, H06), "daC01H'"])
umatApp([ 2.0,gic(A, H06, C01, C02),-1.0,gic(A, H07, C01, C02),-1.0,gic(A, H08,C01, C02), "rC01H"])
umatApp([ 0.0,gic(A, H06, C01, C02), 1.0,gic(A, H07, C01, C02),-1.0,gic(A, H08,C01, C02), "rC01H'"])
umatApp([1.0,gic(D, H06, C01, C02, N03), 1.0,gic(D, H07, C01, C02, N03),
         1.0,gic(D, H08, C01, C02,N03), "iC01H"])
#C04
umatApp([1.0,gic(B, H10, C04), 1.0,gic(B, H11, C04), 1.0,gic(B, H12, C04),"ssC04H"])
umatApp([-2.0,gic(B, H10, C04), 1.0,gic(B, H11, C04),  1.0,gic(B, H12, C04),"saC04H"])
umatApp([ 0.0,gic(B, H10, C04),-1.0,gic(B, H11, C04), 1.0,gic(B, H12, C04),"saC04H’"])
umatApp([ 1.0,gic(A, H10, C04, H11), 1.0,gic(A, H11, C04, H12), 1.0,gic(A, H12,C04, H10),
         -1.0,gic(A, H10, C04,C02),-1.0,gic(A, H11, C04, C02),-1.0,gic(A, H12, C04, C02), "dsC04H"])     
  
umatApp([ 1.0,gic(A, H10, C04, H11), 1.0,gic(A, H11, C04, H12),-2.0,gic(A, H12,C04, H10), "daC04H"])
umatApp([ 1.0,gic(A, H10, C04, H11),-1.0,gic(A, H11, C04, H12), 0.0,gic(A, H12,C04, H10), "daC04H'"])
umatApp([ 2.0,gic(A, H10, C04, C02),-1.0,gic(A, H11, C04, C02),-1.0,gic(A, H12,C04, C02), "rC04H"])
umatApp([ 0.0,gic(A, H10, C04, C02), 1.0,gic(A, H11, C04, C02),-1.0,gic(A, H12,C04, C02), "rC04H'"])
umatApp([1.0,gic(D, H10, C04, N03, C02), 1.0,gic(D, H11, C04, N03, C02), 
          1.0,gic(D, H12, C04, N03,C02), "iC04H"])
#H09
umatApp([1.0,gic(B, H09, N03), "sH09"])
umatApp([1.0,gic(A, H09, N03, C02),-1.0,gic(A, H09, N03, C04), "rH09"])
umatApp([1.0,gic(M, H09, C02, C04, N03), "wH09"]) 
#O05 
umatApp([1.0,gic(B, O05, C02), "sO05"]) 
umatApp([1.0,gic(A, O05, C02, C01),-1.0,gic(A, O05, C02, N03), "rO05"]) 
umatApp([1.0,gic(D, O05, C01, N03, C02), "wO05"]) 
#Rest 
umatApp([1.0,gic(B, C02, C01), "B_C02_C01"]) 
umatApp([1.0,gic(B, N03, C02), "B_N03_C02"]) 
umatApp([1.0,gic(B, C04, N03), "B_C04_N03"]) 
umatApp([1.0,gic(A, N03, C02, C01), "A_N03_C02_C01"]) 
umatApp([1.0,gic(A, C04, N03, C02), "A_C04_N03_C02"]) 
umatApp([1.0,gic(D, C04, N03, C02, C01), "D_C04_N03_C02_C01"])

Parameterization of Bonds Angles and Dihedrals

Bonds, angles and partially dihedral terms are parameterized against QC calculated equilibrium geometry and mode of vibration. mmfqan.py uses CHARMM program for MM energy minimization and symbolic PED analysis of MM calculated modes of vibration. It compares the QC and MM optimized geometries and modes of vibration. Small output from mmfqan.py is shown below:

Comparing all ICs:
IC                Ref     Val     Dif   Acc
Bonds
H11-C04           1.0934  1.1121  0.019 True
H06-C01           1.0931  1.1102  0.017 True
H08-C01           1.0936  1.1102  0.017 True
H10-C04           1.0894  1.1142  0.025 True
N03-C02           1.3674  1.3397 -0.028 True
N03-H09           1.0101  0.9932 -0.017 True
N03-C04           1.4499  1.4450 -0.005 True
C02-O05           1.2327  1.2230 -0.010 True
C02-C01           1.5157  1.4810 -0.035 False
C01-H07           1.0905  1.1119  0.021 True
C04-H12           1.0944  1.1142  0.020 True
Angles
H06-C01-H08       108.74  109.58   0.84 True
H06-C01-C02       113.09  110.22  -2.87 True
H06-C01-H07       109.43  108.80  -0.63 True
H08-C01-C02       109.13  110.22   1.09 True
H07-C01-H08       107.96  108.80   0.84 True
H07-C01-C02       108.35  109.18   0.82 True
O05-C02-N03       123.03  121.28  -1.75 True
C01-C02-N03       114.96  117.05   2.08 True
C01-C02-O05       121.99  121.68  -0.31 True
C02-N03-H09       118.84  119.72   0.88 True
C02-N03-C04       122.11  121.89  -0.22 True
H09-N03-C04       118.98  118.39  -0.59 True
H10-C04-H11       109.41  108.34  -1.07 True
N03-C04-H11       110.81  110.69  -0.12 True
H11-C04-H12       109.13  108.34  -0.78 True
N03-C04-H10       106.91  110.67   3.76 False
H10-C04-H12       109.11  108.03  -1.08 True
N03-C04-H12       111.42  110.67  -0.75 True
Dihedrals
H09-N03-C02-O05   178.39 -180.00   1.61 True
H09-N03-C02-C01    -0.05   -0.00   0.05 True
C04-N03-C02-O05     1.52    0.00  -1.52 True
C04-N03-C02-C01  -176.91  180.00  -3.09 True
C02-N03-C04-H11  -126.03 -180.00 -53.97 False
C02-N03-C04-H10    -6.87  -59.87 -53.00 False
C02-N03-C04-H12   112.26   59.87 -52.40 False
H09-N03-C04-H11    57.11   -0.00 -57.11 False
H09-N03-C04-H10   176.27  120.13 -56.14 False
H09-N03-C04-H12   -64.60 -120.13 -55.54 False
N03-C02-C01-H06   -22.12  -60.54 -38.42 False
N03-C02-C01-H08    99.05   60.54 -38.51 False
N03-C02-C01-H07  -143.62  180.00 -36.38 False
O05-C02-C01-H06   159.43  119.46 -39.97 False
O05-C02-C01-H08   -79.40 -119.46 -40.06 False
O05-C02-C01-H07    37.93   -0.00 -37.93 False
Impropers
C02-C01-N03-O05     0.81   -0.00  -0.81 True 
#  FreqQC sPEDQC           PerQC FreqMM sPEDMM           PerMM
 1   26.6 iC04H               98   63.7 iC01H               99
 2   51.8 iC01H               95   98.7 iC04H               95
 3  144.6 D_C04_N03_C02_C01  100  183.3 D_C04_N03_C02_C01   98
 4  279.9 A_C04_N03_C02       53  263.2 A_C04_N03_C02       50
          A_N03_C02_C01       36        A_N03_C02_C01       29
 5  386.6 wH09                91  433.5 A_N03_C02_C01       49
 6  414.8 A_N03_C02_C01       39  543.6 rO05                51
          rO05                31        B_C02_C01           22
 7  583.7 wO05                46  573.3 wH09                88
          rO05                17        wO05                19
 8  610.3 wO05                36  715.4 wO05                63
          rO05                23        wH09                15
 9  839.6 B_N03_C02           29  771.5 B_N03_C02           30
          rC04H               21        sO05                20
          B_C02_C01           17        B_C02_C01           19
10  965.1 rC01H               49  948.3 B_C04_N03           34
          B_C02_C01           20        rC01H'              26
11 1022.5 rC01H'              62  995.6 rC01H'              33
          wO05                19        B_C04_N03           28
12 1081.6 B_C04_N03           58 1056.2 rC04H               44
                                        daC04H'             19
                                        rC04H'              15
13 1117.1 rC04H'              82 1080.5 rC01H               43
                                        wO05                18
14 1128.7 rC04H               38 1086.5 rC04H'              51
          A_C04_N03_C02       19        rC04H               17
                                        A_C04_N03_C02       15
15 1238.0 rH09                31 1252.5 rH09                38
          B_N03_C02           24        B_C02_C01           27
16 1368.9 dsC01H              89 1384.3 dsC01H              94
17 1401.5 dsC04H              88 1412.3 daC01H              73
                                        daC04H              17
18 1444.0 daC01H              88 1415.8 daC04H              78
                                        daC01H              16
19 1460.0 daC01H'             81 1418.4 daC01H'             91
20 1466.8 daC04H'             61 1427.3 daC04H'             76
                                        rC04H               18
21 1474.2 daC04H              87 1474.0 dsC04H              43
                                        rH09                28
22 1516.9 rH09                36 1580.9 dsC04H              46
          B_N03_C02           21        rH09                19
          rC04H               15        B_C04_N03           17
23 1696.7 sO05                78 1676.1 sO05                67
24 2926.4 ssC04H              94 2852.7 ssC04H             100
25 2935.4 ssC01H              98 2913.9 saC04H’             74
                                        saC04H              25
26 2997.6 saC04H’             99 2914.7 saC04H              75
                                        saC04H’             25
27 3017.2 saC01H’            100 2917.0 ssC01H             100
28 3042.4 saC01H              98 2975.1 saC01H’            100
29 3053.0 saC04H              95 2975.3 saC01H             100
30 3463.8 sH09               100 3325.1 sH09                99 

Dihedral fitting plug-in to VMD

Conclusion

• ParamIT improves the parameterization process by automating many routine steps. The interactivity of graphical interfaces enhances and accelerates the optimization of partial charges and dihedral term. Overall, ParamIT significantly reduces the manual labor and sources of human input error.
• The toolkit can be used for complete parameterization or in conjunction with ParamChem, a web-service for parameterization by analogy, for partial reoptimization.

References

• Brooks, B. R., C. L. Brooks, et al. (2009). "CHARMM: The biomolecular simulation program." Journal of Computational Chemistry 30(10): 1545-1614.
• Frisch, M. J. T., G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. (2009). Gaussian 09, Revision A.02. Gaussian, Inc., Wallingford CT.
• Guvench, O. and A. D. MacKerell, Jr. (2008). "Automated conformational energy fitting for force-field development." Journal of molecular modeling 14(8): 667-679.
• Humphrey, W., A. Dalke, et al. (1996). "VMD: Visual molecular dynamics." J. Mol. Graphics 14(1): 33-38.
• MacKerell, A. D., D. Bashford, et al. (1998). "All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins†." The Journal of Physical Chemistry B 102(18): 3586-3616.
• Vanommeslaeghe, K., E. Hatcher, et al. (2010). "CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields." Journal of computational chemistry 31(4): 671-690.