VFFDT: A New Software for Preparing AMBER Force Field Parameters

Subscriber access provided by Mount Allison University | Libraries and Archives

Article

VFFDT: A New Software for Preparing AMBER Force Field Parameters for Metal-Containing Molecular Systems Suqing Zheng, Qing Tang, Jian He, Shiyu Du, Shaofang Xu, Chaojie Wang, Yong Xu, and Fu Lin J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.5b00687 • Publication Date (Web): 21 Mar 2016 Downloaded from http://pubs.acs.org on March 24, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling Manuscript for “Software Description section” in the Journal of Chemical Information and Modeling

VFFDT: A New Software for Preparing AMBER Force Field Parameters for Metal-Containing Molecular Systems Suqing Zhenga, Qing Tanga, Jian Heb, Shiyu Duc ,Shaofang Xua, Chaojie Wanga, Yong Xud, and Fu Lina* a

: School of Pharmaceutical Sciences, Wenzhou Medical University, Wenzhou, 325035, P. R. China

b

: Center for Translational Medicine, Department of Biotechnology, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, 457 Zhongshan Rd, Shahekou， Dalian, Liaoning, 116023, P. R. China c

: Engineering Laboratory of Specialty Fibers and Nuclear Energy Materials, Ningbo Institute of Materials Technology and Engineering, Chinese Academy of Sciences, 519 Zhuangshi Ave, Zhenhai, Ningbo, Zhejiang 315201, P. R. China d

: Center of Chemical Biology, Guangzhou Institutes of Biomedicine and Health, Chinese Academy of Sciences, No.190 Kaiyuan Avenue, Guangzhou Science Park, Guangzhou, Guangdong 510530, P. R. China

*: To whom all correspondence should be addressed. E-mail: [email protected];

1

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 29

Manuscript for “Software Description section” in the Journal of Chemical Information and Modeling

ABSTRACT

Force fields are fundamental to molecular dynamics simulations. However, the incompletion of force field parameters has been a long standing problem, especially for the metal-related systems. In our previous work, we have adopted the Seminario method based on the hessian matrix to systematically derive the Zinc related force field parameters for AMBER. In this work, in order to further simplify the whole protocol, we implement a user-friendly Visual Force Field Derivation Toolkit (VFFDT) to derive the force field parameters via simply clicking on the bond or angle in the 3D viewer, and further extend our previous program to support the hessian matrix input from varieties of QM packages such as Gaussian03/09, ORCA3.0, QChem, GAMESS-US, and MOPAC2009/2012. In this toolkit, a universal VFFDT XYZ file format containing the raw hessian matrix is available for all the QM packages and an instant force field parameterization protocol based on the semi-empirical quantum mechanics (SQM) method is introduced. The new function to automatically obtain the relevant parameters for zinc, copper, iron etc., which can be exported as AMBER Frcmod format is added. Furthermore, our VFFDT program can read and write the AMBER Prepc AMBER Frcmod, and AMBER Mol2 format file, and can also customize, view, copy and paste the force field parameters in the context of the 3D viewer, which provides complementary utilities to the ANTECHAMBER, MCPB, and MCPB.py in the AmberTools.

2


Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION It is well known that Molecular Dynamics (MD) simulation can provide novel insight into the microscopic process in biological and material science. The predictability of the MD simulation mainly relies on the accuracy of the force field, which is described by a group of force field parameters. Therefore, the determination of force field parameters is one of the crucial steps in a classical MD simulation. However, in a practical simulation setup, we will often encounter parameter missing problem for various small molecules or metal-containing molecules. For small organic molecules, several programs such as ANTECHAMBER1, 2 in the AMBER community, CCGEN3, 4

in the CHARMM community is frequently used to generate general force field

parameters, while, for metal-containing molecules, it is much more complicated due to the diversity of force field models such as bonded model,5-13 semi-bonded model,14, 15

non-bonded model,16-18 or polarization model.19,

20

The pros and cons of these

models have been extensively addressed in the previous literatures.21 The bonded model cannot well capture the change of coordination number around the metal during the simulation process, and also this model fails in the coordination systems with two or more water molecules as pointed out by Ryde group.21 however, the bonded model is cheap from computational perspective compared to the polarization model and compatible with the existing pairwise additive force field, and more importantly it can reproduce the geometry in the relatively rigid metal-containing systems, which are common in most cases.21-23 In this work the force field parameters will be derived in the context of bonded model in the AMBER force field, which was adopted by our previous work.22 To derive the force field parameters for the metal-containing systems, basically four categories of methods have been proposed to deal with this problem: the Potential Energy Scan (PES) method, the internal coordinate method, the Seminario method,24 and the Ayers method.25 The PES method is extensively used in the parameterization of small organic molecules or metal-containing systems. Neves et al. gave a good 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29


example to adopt PES method to parameterize the Manganese-Containing system with the semi-flexible approach.26 However, it is very tedious and computationally expensive to optimize all bond and angle related parameters for the specific metal-containing systems. The internal coordinate method is to directly derive force constants from the internal coordinate based hessian matrix after the QM calculation. It is simple and convenient. Nevertheless, the parameters derived with this method is highly dependent on the internal coordinate that user defined in the QM computations, which means that different internal coordinate schemes will obtain different force constant values.24 The Seminario method24 and Ayers method25 are proposed to adopt the Cartesian based hessian matrix from QM calculations so as to avoid the internal coordinate dependency. In the Ayers method, the hessian matrices from QM calculation and Molecular Mechanics (MM) prediction with target force field are both required to derive the bond and angle related force field parameters in order to achieve the overall consistency. Hence, this method actually needs a preliminary estimation of force field parameters in advance to get the rough hessian matrix from MM based calculation. As for the Seminario method, only one hessian matrix from QM computation is needed to derive all the bond and angle related force field parameters. Due to its simplicity, the Seminario method was frequently adopted to derive the metal related parameters based on the QM calculation, even though this method is initially proposed to deal with very small molecule such as water and nitrogen dioxide in the original paper24 and is applied to the larger systems such as peptide or nitro compounds in their recent work.27, 28 Ryde group developed the Hess2FF program29 based on the Seminario method and parameterized several metal-containing systems, and then they followed the parameterization scheme of Norrby and Liljefors (NL)30 to avoid the double-counting of nonbonded interactions and to reproduce the QM optimized geometry and hessian matrix elements from the force field parameters given by the Seminario method through the iterative optimization. Comparative 4


Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


evaluation shows that the NL method is more accurate than the Seminario method.21 Inspired by their work, we developed a protocol to systematically derive zinc related force field parameters compatible with AMBER, which was also based on the Seminario method and could subsequently be optimized to reproduce the vibrational frequencies via genetic algorithm.22 Almost at the same time, Merz group also independently parameterized the zinc-containing systems based on the Seminario method and released a powerful tool named Metal-Center Parameter Builder (MCPB),23 which was integrated in the AMBER15.31 Ryde group did systematic evaluation about the performance of zinc related force field parameters and demonstrated that Seminario method or its variants can offer a good starting point for the force field derivation or further force field optimization of metal-containing systems, albeit concluded that the parameterization of metal-containing systems is seldom an automatic process.21 Based on the Seminario method, the Hess2FF29 and MCPB toolkits23 were developed to determine the metal-related force field parameters. Both are aimed to automate the parameterization procedure for the metal-containing systems and supports the Gaussian log file or formatted checkpoint file. Relative to HessFF, MCPB is well documented, actively maintained by Merz group, and is rewritten with Python and renamed as MCPB.py32 with optimized code and fewer steps in practice relative to MCPB. But the user interface is not friendly enough, which makes it difficult for the beginners to follow the whole protocol. In this work, we will present a user-friendly Visual Force Field Derivation Toolkit (VFFDT) based on the Seminario method, which allows the users to derive the bond or angle related force field parameters for metal-containing systems or small molecules by simply clicking on a relevant bond or angle, and help users to customize, view, copy and paste the force field parameters via the interactive way. The function of VFFDT is complementary to that of MCPB, MCPB.py or ANTECHAMBER in AmberTools. In this paper, a brief description of AMBER potential function and VFFDT program will be given in the Methods section 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29


and the detailed features of our VFFDT program will be described in the Results and Discussion section.

METHODS 1. Brief review of AMBER potential function In most of typical pairwise additive force fields such as AMBER force field, the potential energy is usually decomposed into five terms: bond stretching, angle bending, torsion angle rotation, van der Waals and electrostatic interaction. For the metal-containing systems, the Van der Waals radius and well depth of common metal ion are generally available in the corresponding force field library and the optimized van der Waals parameters can be retrieved in the literature,

33-36

and the point charge

of metal in the bonded model can be easily obtained by the mature protocol such as RESP method.37 Torsion angle rotation related with metal can often be ignored.5-12 As for bond stretching and angle bending, the corresponding parameters, especially force constants, cannot be found in its intrinsic force field library in most occasions and pose a major problem for the practical molecular dynamics simulation. These parameters could be derived by the Seminario method. A brief review of this method is given as follows for the convenience of readers in the supporting information (Text S1). More details can be found in the original paper.24 2. Extraction Hessian matrix from the output of various QM packages The full hessian matrix as the only input for the Seminario method can be extracted from frequency analysis based on the QM optimized structure with the same QM method. The flag for extracting the Cartesian based hessian matrix varies in different QM packages. Our Visual Force Field Derivation Toolkit (VFFDT) supports the output file of various QM packages such as Gaussian03/09,38 QChem,39 6


Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MOPAC2009/2012,40 GAMESS-US,41 and ORCA3.0,42 of which the last three programs are available free of charge for non-profit or academic users. Although VFFDT can support those QM packages very well, it cannot cover all the QM packages due to the booming of new QM softwares nowadays. However, VFFDT offers a universal format containing both coordinate and raw hessian matrix information called VFFDT XYZ format, which is actually a modified version of XYZ format file with an extra sections for AMBER atomic type and raw hessian matrix. The detailed format will be given in the section of results and discussion. 3. Force field derivation on-the-fly with Semi-empirical based hessian matrix In order to derive the force field parameters on-the-fly, an automatic protocol called “Instant FF Derivation” has been implemented in VFFDT under the menu “Protocol”, which integrates the geometry optimization and frequency analysis with the semi-empirical quantum mechanics (SQM) method, and force field derivation based on the Seminario method. In the VFFDT program, the SQM method in Gaussian03/09 will be used to generate the hessian matrix analytically, rather than by the numerical method as in MOPAC2009/2012.

RESULTS AND DISCUSSION Our VFFDT software is developed to derive bond stretching and bond angle bending force field parameters for metal-containing molecular systems based on the Seminario method. The main interface of VFFDT is displayed in Figure 1. The basic functions of this software include: (1) 3D operation on a molecular structure, such as translation, rotation, and scaling. (2) Rendering a molecule in different styles, such as the ball-and-stick model and the space-filling model. (3) Reading the outputs of Gaussian03/09, QChem, MOPAC2009/2012, GAMESS-US, and ORAC3.0. (4) Outputting molecular structures in the SYBYL Mol2 format, the RCSB PDB format, 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 29


the standard XYZ format and the VFFDT XYZ format. (5) Setting up parameters for a Gaussian job. (6) The parameterization module for a metal-containing system. (6) View, modify, copy and paste the force field parameters in the 3D viewer. 1. Parameterization module for metal-containing system In VFFDT, there are basically three modes to derive force field parameters. (1) Users can click on the bond or angle in the viewer window (Figure 1), while the corresponding parameters related with bonds or angles will be printed in the console window (Figure S1). (2) Users can input the command in the console window with keywords “Pick Atom_ID1 Atom_ID2” or “Pick Atom_ID1 Atom_ID2 Atom_ID3”, by which the force field parameters can also be obtained (Figure 2). (3) Users can click on the menu “Protocol”, which will provide the force field derivation of several common metals or user-specified metals. In this function, all the bond or angle parameters related with metals will be derived and listed in the table (Figure 3 and Figure 2S). All the atoms involved with the bond or angle containing a metal atom can be assigned the AMBER atomic type automatically, which can be exported to AMBER Frcmod format file. It should be noted that, before deriving the force field parameter, the users can set the frequency scaling factor based on the basis set in the QM computation (Figure 4). This scaling factor is originally used to better reproduce the thermodynamic quantity and is adopted to calibrate the force constant in the VFFDT program. The default value of frequency scaling factor is 1.000 indicating that no correction will be employed. To vividly describe the function above, we take the [Zn(imidazole)4]2+ as instance and aim to derive the bond stretching force constant of Zn-N, which is often missing in a common force field library. Suppose that the output from the QM calculation of vibrational frequencies with Gaussian program has already been obtained, the procedure

is

as

follows.

(1)

Load

the

Gaussian

output

“Zn-his-his-his-his-freq-6-311++G-2d2p.out” into the viewer window of VFFDT 8


Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(Figure 1). (2) Click on the menu “Protocol Freq Scaling Factor” and input the frequency scaling factor (Figure 4). (3) Click on the icon

on the toolbar and select

two atoms in the viewer window(Figure 1S), the bond stretching force constant of Zn-N will be automatically shown in the console window, which can be open via clicking on the menu “Windows Show Log”; The bond stretching force constant can also be directly derived in the console window with the following command “>>pick 1 18”

or “>>pick 18 1” (1 and 18 is atom ID for these two atoms in

specific bond, Figure 2); (4) Click on the menu “Protocol FF For the Metal Containing System Only Zn Force Field”, a window titled “Zn Related Force Field Parameters” will pop out (Figure 3). All the bond stretching force constants related with zinc will be derived by clicking on the button “Generate All Parameters Automatically” (Figure 3 and Figure S2). Meanwhile, the AMBER atomic type can be assigned automatically if user clicks on the button “Assign AMBER GAFF Atom Type” (Figure 3). 2. Support various QM packages VFFDT supports various QM packages such as commercial softwares: Gaussian03/09 and QChem, or those free for academic use: MOPAC2009/2012, GAMESS-US, and ORAC3.0. In the calculation of vibrational frequencies, different QM packages with their specific keywords will dump different data formats for the hessian matrix. Accordingly, VFFDT have been coded separately to be compatible to the programs in the support list. For instance, frequency analysis in Gaussian09 with keyword “freq” or that in Gaussian03 with both keywords “freq” and “iop(7/33=1)” will generate the Cartesian based hessian matrix tagged with “Force constants in Cartesian coordinates:”, which is not mass-weighted. If the QM packages are not in our support list, the VFFDT offers the customized VFFDT XYZ format, which is the variant of standard XYZ format. In this format, there are three sections. The first section is for the atomic symbol and XYZ 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29


coordinate (unit: Å) for each atom, which is fully compatible with the original XYZ format; the second section defines the AMBER atomic type of each atom; and the third section contains raw hessian matrix, which could be full hessian matrix or half hessian matrix in the lower triangle where the unit of each element is Hartree/Bohr2. The example of this format has been given in Figure 5. Therefore, users can easily follow the VFFDT XYZ format to prepare an input file for VFFDT if their frequently-used QM package is not in our support list. To the best of our knowledge, the current version of MCPB or MCPB.py in AmberTools and Hess2FF program do not have this flexible option yet and supports only limited QM packages. 3. Force field derivation on-the-fly with Semi-empirical based hessian matrix An “Instant FF Derivation” protocol has also been implemented in VFFDT to quickly derive the force field with the Semi-empirical Quantum Mechanics (SQM) method. The force field parameters from this protocol may be more efficient for the metal-containing systems, since SQM such as PM643 in the Gaussian program is well parameterized, albeit more approximations have been applied in the SQM method relative to the high level QM method. Therefore, this inexpensive method could provide force field parameters for the fast molecular docking with Molecular Mechanics (MM) based scoring function in the time-urgent drug discovery project, or could offer preliminary force field parameters before further optimization to reproduce the potential energy surfaces or specific properties for systems of interest. Meanwhile, this protocol may also be useful for the parameterization of small-sized or medium-sized organic molecules, which are well described by the SQM method. As we know, ANTECHAMBER in AMBER is a versatile tool to generate General Amber Force Field (GAFF) for organic molecules without metal element. Nevertheless, GAFF will ignore some chemical environment for the specific molecule due to its general use and transferability, although the various atomic types has been defined in the GAFF aiming to cover as many chemical environments as possible. The 10


Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bond or angle related parameters, quickly derived from the protocol in VFFDT, take the chemical environment into account and treat each atom in each molecule case by case. It can help users to determine whether it is necessary to define a new atomic type to consider the special chemical environment. Thus this function implemented in VFFDT can be a complementary to ANTEAMBER for AMBER. Here we take the benzene molecule as an example to show how all the bond and angle related force field parameters are generated in VFFDT. Since the detailed tutorial has been given in the help document of VFFDT program, only a concise description of the procedures will be provided here. (1) Load Mol2 file “benzene.mol2” in the VFFDT viewer window (Figure S3). (2) Click on the menu “ProtocolInstant FF Derivation” (Figure 6), then a panel titled “Instant FF Derivation” pops out (Figure S4) and waits for the users to input the charge and multiplicity of the molecule in the viewer window. (3) Click on the command “Submit Job to Localhost” (Figure S4) and the window titled “Submit Job to Localhost” will open and show the status of the job (Figure S5). (4) Double-click on the cell with column header “JobLogDirectory” (Figure S5) once the job is finished, the Gaussian output “benenze.out” will be automatically loaded into the VFFDT viewer window. (5) Click on the menu “ProtocolFreq Scaling Factor” and input the frequency scaling factor 1.062 that is the precomputed vibrational scaling factor for PM6 method from Computational Chemistry Comparison and Benchmark DataBase (CCCBDB) on the website (http://cccbdb.nist.gov/vibscalejust.asp). (5) Click on the menu “Protocol FF for the Whole General Small Molecule” (Figure S6), the panel titled “All Related Force Field Parameters” can assign the AMBER GAFF atomic type automatically and derive all the parameters for the bond and angle terms for benzene (Figure S6). (6) Click on the command “Export as AMBER Frcmod” to export the content in the textbox to the AMBER Frcmod file, or Click the “Map Parameters in textbox to Molecule” to map the averaged parameters in the textbox or manually

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29


modified parameters in this editable textbox to the molecule in the VFFDT viewer window (Figure 7). In order to test the effects of different chemical environments on force field parameters, the bond stretching force constants of the carbon-carbon bond in the aromatic ring of benzene, phenol and aniline were derived with the same protocol as mentioned above in VFFDT and listed in the Figure 8. In the benzene, phenol and aniline molecules, the AMBER atomic type for all the carbon atoms in the aromatic ring is “ca”, hence the bond stretching force constants for all the carbon-carbon bond in the aromatic ring of three molecules are all assigned to 478 kcal/(mol·Å2) by the ANTECHAMBER program. However, the carbon atom attached to oxygen atom in phenol or nitrogen atom in aniline is different from the carbon atom in benzene. Therefore, the oxygen or nitrogen may affect the neighboring bond more or less. From Figure 8, the bond stretching force constant of the bond 1-2 (1 and 2 are atom IDs shown in Figure 8) derived from with the protocol “Instant FF Derivation” in the VFFDT is 382kcal/(mol·Å2), 369kcal/(mol·Å2) and 344kcal/(mol·Å2) for benzene, phenol, and aniline respectively, which displays that the oxygen or nitrogen atom has an impact on the bond stretching force constant of the nearby bond. It should be mentioned that those force constants of the bond 1-2 are computed with the keyword “#p opt freq=noraman pm6” and frequency scaling factor 1.062, which is from CCCBDB database as mentioned above. However, in the parameterization scheme of GAFF, geometry optimization was performed in the level of “mp2/6-31g*”, for consistency and accuracy the QM computation with keyword “#p opt freq=noraman mp2/6-31g*” and frequency scaling factor 0.943 also from CCCBDB database is used to derive the bond stretching force constant of the bond 1-2 in VFFDT, which are 397kcal/(mol·Å2), 382kcal/(mol·Å2) and 366kcal/(mol·Å2) for benzene, phenol, and aniline respectively (Figure 8). This result also indicates the different chemical environment will affect the bond stretching force constant and suggests that it may be necessary to introduce a new AMBER atomic type to describe the carbon atom 12


Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


attached to the oxygen or nitrogen atom. Nevertheless, it may influence the transferability of force field parameters if too many new AMBER atomic type is defined based on different chemical environment, which needs further validation in the future work. Thus in this work VFFDT program will provide the alternative option for the users to derive force field parameters that takes the chemical environment into account for the metal-containing systems or small molecules, that is complementary to the ANTECHAMBER programs. 4. Various VFFDT utilities complementary to the AmberTools VFFDT provides some extra utilities that can complement the MCPB, MCPB.py, ANTECHAMBER or LEAP in the AmberTools.

(1) Read and write the AMBER

Prepc, AMBER Frcmod and AMBER Mol2 Format files (Figure S7). (2) View force field parameters in the context of 3D viewer window. (3) Modify the AMBER atomic type handily (Figure S8). (4) Copy/Paste parameters from other molecule in the 3D viewer window (Figure S9). One tutorial combining the VFFDT, ANTECHAMBER and LEAP is given in the VFFDT manual to illustrate its complementary usage, which has been well illustrated in the help document of VFFDT program. 5. Software availability VFFDT is coded with C++ and OpenGL. It has been tested on Microsoft Windows platforms, including Windows XP, Windows 7, Windows 8 and Windows 10. Three tutorials and their related input files have been integrated in the installation package of VFFDT. The VFFDT program and its PDF manual can be obtained freely from the author upon request or from the following Baidu cloud disk and Dropbox shared folder. All the future version and bug report will be also updated on the pubic available links as follows. Baidu Cloud Disk: http://pan.baidu.com/share/home?uk=2642812499

Dropbox shared folder: 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 29


https://www.dropbox.com/sh/hlefz8xy8o23q44/AADZGTBN-aJ07flZ-kM0s4U6a?dl=0

CONCLUSIONS In this work, we developed Visual Force Field Derivation Toolkit (VFFDT) based on the Seminario method. VFFDT can easily derive the bond or angle related force field parameters via simple clicking on the bond or angle, support various QM packages, provide a universal VFFDT XYZ file format and offer a force field parameterization protocol on-the-fly based on the SQM method. Additionally, our VFFDT program can handle AMBER Prepc, AMBER Frcmod, and AMBER Mol2 format file and can also customize, view, copy and paste the force field parameters in the 3D viewer window. We believe that this tool can be complementary to the ANTECHAMBER, MCPB, and MCPB.py programs in AmberTools. ACKNOWLEDGMENT The authors acknowledge the support of the startup funding of Wenzhou Medical University, the National Natural Science Foundation of China (Grant No. 21502144, 51372046, 51479037 and 91226202), Zhejiang Natural Science Foundation (Grant No. LY16B070006), the Ningbo Municipal Natural Science Foundation (Grant No. 2014A610006), ITaP at Purdue University for computing resources and the key technology of nuclear energy, 2014, CAS Interdisciplinary Innovation Team. The authors thank Prof. Renxiao Wang at Shanghai Institute of Organic Chemistry, Chinese Academy of Science in China for his kind revision and suggestion. The authors also thank suggestions from the reviewers. Supporting Information

14


Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The brief review of Seminario method is given in the Text S1 and some figures used to describe the functions of VFFDT program have been in the Figure S1-S9. This material is available free of charge via the Internet at http://pubs.acs.org.

REFERENCES 1. Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 2006, 25, 247-260. 2. Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157-1174. 3. Vanommeslaeghe, K.; MacKerell, A. D., Jr. Automation of the CHARMM General Force Field (CGenFF) I: bond perception and atom typing. J. Chem. Inf. Model. 2012, 52, 3144-3154. 4. Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D., Jr. Automation of the CHARMM General Force Field (CGenFF) II: assignment of bonded parameters and partial atomic charges. J. Chem. Inf. Model. 2012, 52, 3155-3168. 5. Falconi, M.; Altobelli, G.; Iovino, M.; Politi, V.; Desideri, A. Molecular dynamics simulation of Matrix Metalloproteinase 2: fluctuations and time evolution of recognition pockets. J. Comput. -Aided Mol. Des. 2003, 17, 837-848. 6. Li, W.; Zhang, J.; Wang, J.; Wang, W. Metal-coupled folding of Cys2His2 zinc-finger. J. Am. Chem. Soc. 2008, 130, 892-900. 7. Lu, D.; Voth, G. A. Molecular dynamics simulations of human carbonic anhydrase II: Insight into experimental results and the role of solvation. Proteins: Struct., Funct., Bioinf. 1998, 33, 119-134. 8. Lu, Q.; Tan, Y.-H.; Luo, R. Molecular Dynamics Simulations of p53 DNA-Binding Domain. J. Phys. Chem. B 2007, 111, 11538-11545. 9. Merz, K. M., Jr. Insights into the function of the zinc hydroxide-Thr199-Glu106 hydrogen bonding network in carbonic anhydrases. J. Mol. Biol. 1990, 214, 799-802. 10. Suárez, D.; Merz, K. M., Jr. Molecular Dynamics Simulations of the Mononuclear Zinc-β-lactamase from Bacillus Cereus. J. Am. Chem. Soc. 2001, 123, 3759-3770. 11. Tuccinardi, T.; Martinelli, A.; Nuti, E.; Carelli, P.; Balzano, F.; Uccello-Barretta, G.; Murphy, G.; Rossello, A. Amber force field implementation, molecular modelling study, synthesis and MMP-1/MMP-2 inhibition profile of (R)- and (S)-N-hydroxy-2-(N-isopropoxybiphenyl-4-ylsulfonamido)-3-methylbutanamides. Bioorg. Med. Chem. 2006, 14, 4260-4276. 12. Vedani, A.; Huhta, D. W. A new force field for modeling metalloproteins. J. Am. Chem. Soc. 1990, 112, 4759-4767. 15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13. Santos-Martins, D.; Forli, S.; Ramos, M. J.; Olson, A. J. AutoDock4Zn: An Improved AutoDock Force Field for Small-Molecule Docking to Zinc Metalloproteins. J. Chem. Inf. Model. 2014, 54, 2371-2379. 14. Pang, Y. P. Novel Zinc Protein Molecular Dynamics Simulations: Steps Toward Antiangiogenesis for Cancer Treatment. J. Mol. Model. 1999, 5, 196-202. 15. Pang, Y. P.; Xu, K. U. N.; Yazal, J. E.; Prendergast, F. G. Successful molecular dynamics simulation of the zinc-bound farnesyltransferase using the cationic dummy atom approach. Protein Sci. 2000, 9, 1857-1865. 16. Donini, O. A.; Kollman, P. A. Calculation and prediction of binding free energies for the matrix metalloproteinases. J. Med. Chem. 2000, 43, 4180-4188. 17. Stote, R. H.; Karplus, M. Zinc binding in proteins and solution: a simple but accurate nonbonded representation. Proteins 1995, 23, 12-31. 18. Wu, R.; Lu, Z.; Cao, Z.; Zhang, Y. A Transferable Nonbonded Pairwise Force Field to Model Zinc Interactions in Metalloproteins. J. Chem. Theory Comput. 2011, 7, 433-443. 19. Sakharov, D. V.; Lim, C. Zn protein simulations including charge transfer and local polarization effects. J. Am. Chem. Soc. 2005, 127, 4921-4929. 20. Zhu, T.; Xiao, X. D.; Ji, C. G.; Zhang, J. Z. H. A New Quantum Calibrated Force Field for Zinc-Protein Complex. J. Chem. Theory Comput. 2013, 9, 1788-1798. 21. Hu, L. H.; Ryde, U. Comparison of Methods to Obtain Force-Field Parameters for Metal Sites. J. Chem. Theory Comput. 2011, 7, 2452-2463. 22. Lin, F.; Wang, R. X. Systematic Derivation of AMBER Force Field Parameters Applicable to Zinc-Containing Systems. J. Chem. Theory Comput. 2010, 6, 1852-1870. 23. Peters, M. B.; Yang, Y.; Wang, B.; Fusti-Molnar, L.; Weaver, M. N.; Merz, K. M., Jr. Structural Survey of Zinc Containing Proteins and the Development of the Zinc AMBER Force Field (ZAFF). J. Chem. Theory Comput. 2010, 6, 2935-2947. 24. Seminario, J. M. Calculation of intramolecular force fields from second-derivative tensors. Int. J. Quantum Chem. 1996, 60, 1271-1277. 25. Burger, S. K.; Lacasse, M.; Verstraelen, T.; Drewry, J.; Gunning, P.; Ayers, P. W. Automated Parametrization of AMBER Force Field Terms from Vibrational Analysis with a Focus on Functionalizing Dinuclear Zinc(II) Scaffolds. J. Chem. Theory Comput. 2012, 8, 554-562. 26. Neves, R. P. P.; Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. Parameters for Molecular Dynamics Simulations of Manganese-Containing Metalloproteins. J. Chem. Theory Comput. 2013, 9, 2718-2732. 27. Bautista, E. J.; Seminario, J. M. Harmonic force field for glycine oligopeptides. Int. J. Quantum Chem. 2008, 108, 180-188. 28. Bellido, E. P.; Seminario, J. M. Harmonic force field for nitro compounds. J. Mol. Model. 2012, 18, 2805-2811. 29. Nilsson, K.; Lecerof, D.; Sigfridsson, E.; Ryde, U. An automatic method to generate force-field parameters for hetero-compounds. Acta. Crystallogr. D Biol. Crystallogr. 2003, 59, 274-289. 16


Page 16 of 29

Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


30. Rydberg, P.; Olsen, L.; Norrby, P. O.; Ryde, U. General transition-state force field for cytochrome p450 Hydroxylation. J. Chem. Theory Comput. 2007, 3, 1765-1773. 31. Case, D.A.; Berryman, J.T.; Betz, R.M.; Cerutti, D.S.; Cheatham, III, T.E.; Darden, T.A.; Duke, R.E.; Giese, T.J.; Gohlke, H.; Goetz, A.W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T.S.; LeGrand, S.; Li, P.; Luchko, T.; Luo, R.; Madej, B.; Merz, K.M.; Monard, G.; Needham, P.; Nguyen, H.; Nguyen, H.T. Omelyan, I.; Onufriev, A.; Roe, D.R.; Roitberg, A.; Salomon-Ferrer, R.; Simmerling, C.L.; Smith, W.; Swails, J.; Walker, R.C.; Wang, J.; Wolf, R.M.; Wu, X.; York, D.M.; Kollman, P.A. AMBER, version2015; University of California: San Francisco, CA, 2015. 32. Li, P.; Merz, K. M., Jr. MCPB.py: A Python Based Metal Center Parameter Builder. J. Chem. Inf. Model. 2016. Just Accepted Manuscript. 33. Babu, C. S.; Lim, C. Empirical Force Fields for Biologically Active Divalent Metal Cations in Water. J. Phys. Chem. A 2006, 110, 691-699. 34. Li, P.; Roberts, B. P.; Chakravorty, D. K.; Merz, K. M., Jr. Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for +2 Metal Cations in Explicit Solvent. J. Chem. Theory Comput. 2013, 9, 2733-2748. 35. Li, P.; Song, L. F.; Merz, K. M., Jr. Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model. J. Chem. Theory Comput. 2015, 11, 1645-1657. 36. Li, P.; Song, L. F.; Merz, K. M., Jr. Parameterization of Highly Charged Metal Ions Using the 12-6-4 LJ-Type Nonbonded Model in Explicit Water. J. Phys. Chem. B 2015, 119, 883-895. 37. Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. Application of the multimolecule and multiconformational RESP methodology to biopolymers: Charge derivation for DNA, RNA, and proteins. J. Comput. Chem. 1995, 16, 1357-1377. 38. Frisch, M. J.; Trucks, 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. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; 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. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. 39. Kong, J.; White, C. A.; Krylov, A. I.; Sherrill, D.; Adamson, R. D.; Furlani, T. R.; Lee, M. S.; Lee, A. M.; Gwaltney, S. R.; Adams, T. R.; Ochsenfeld, C.; Gilbert, 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 29


A. T. B.; Kedziora, G. S.; Rassolov, V. A.; Maurice, D. R.; Nair, N.; Shao, Y.; Besley, N. A.; Maslen, P. E.; Dombroski, J. P.; Daschel, H.; Zhang, W.; Korambath, P. P.; Baker, J.; Byrd, E. F. C.; Van Voorhis, T.; Oumi, M.; Hirata, S.; Hsu, C.-P.; Ishikawa, N.; Florian, J.; Warshel, A.; Johnson, B. G.; Gill, P. M. W.; Head-Gordon, M.; Pople, J. A. Q-Chem 2.0: a high-performance ab initio electronic structure program package. J. Comput. Chem. 2000, 21, 1532-1548. 40. Stewart, J. J. Optimization of parameters for semiempirical methods VI: more modifications to the NDDO approximations and re-optimization of parameters. J. Mol. Model. 2013, 19, 1-32. 41. Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347-1363. 42. Neese, F. The ORCA program system. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73-78. 43. Stewart, J. J. P., Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements. J. Mol. Model. 2007, 13, 1173-1213.

18


Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48



Figure 1. The main interface for Visual Force Field Derivation Toolkit (VFFDT).

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 29


Figure 2. Derivation of bond related force field parameter via the command line “pick 1 18” in the console window.

20


Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Derivation of bond related force field parameter via clicking on the menu “ProtocolFF for the Metal Containing Region OnlyZn Force Field”.

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 29


Figure 4. The menu “ProtocolFreq Scaling Factor” and the input box for setting the vibrational frequency scaling factor.

22


Page 23 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. The example of VFFDT XYZ format. The first section is the same as original standard XYZ format; The second section is AMBER atomic type of each atom, which is assigned by the VFFDT program; The third section is about the full/half hessian matrix, which is not mass-weighted. It should be noted that the whole format adopts the free format and the units for atomic coordinate and element in the hessian matrix are Å and Hartree/Bohr2 respectively.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 29


Figure 6. The protocol “ProtocolInstant FF Derivation” used for the deriving the force field parameters on-the-fly.

24


Page 25 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. All the force field parameters will be summarized in the textbox of the panel titled “Statistics For Parameter” via clicking on the command “Statistics” on the window titled “All Related Force Field Parameters” in the Figure S6.

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 29


Figure 8. The bond stretching force constants (unit: kcal/(mol·Å2)) of carbon-carbon bond in the aromatic ring of the benzene, phenol, and aniline molecules respectively. The values in blue are given by “Instant FF Derivation” protocol with SQM method and Seminario method, those in black in the parenthesis are calculated with QM method (MP2/6-31G*) and Seminario method and those in red in the parenthesis are assigned by the ANTECHAMBER program. Note: the frequency scaling factor for the PM6 and MP2/6-31G* are 1.062 and 0.943 26


Page 27 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


respectively.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Page 28 of 29


TOC Graphic 28


Page 29 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



29


VFFDT: A New Software for Preparing AMBER Force Field Parameters

Recommend Documents