Automation of the CHARMM General Force Field (CGenFF) II: Assignment of bonded parameters and partial atomic charges
Molecular mechanics force fields are widely used in computer-aided drug design for the study of drug candidates interacting with biological systems. In these simulations, the biological part is typically represented by a specialized biomolecular force field, while the drug is represented by a matching general (organic) force field. In order to apply these general force fields to an arbitrary drug-like molecule, functionality for assignment of atom types, parameters and partial atomic charges is required. In the present article, algorithms for the assignment of parameters and charges for the CHARMM General Force Field (CGenFF) are presented. These algorithms rely on the existing parameters and charges that were determined as part of the parametrization of the force field. Bonded parameters are assigned based on the similarity between the atom types that define said parameters, while charges are determined using an extended bond-charge increment scheme. Charge increments were optimized to reproduce the charges on model compounds that were part of the parametrization of the force field. A \xe2\x80\x9cpenalty score\xe2\x80\x9d is returned for every bonded parameter and charge, allowing the user to quickly and conveniently assess the quality of the force field representation of different parts of the compound of interest. Case studies are presented to clarify the functioning of the algorithms and the significance of their output data.
<!>Introduction<!>Approach and algorithms<!>Assignment of parameters by analogy<!>Total penalty score (TPS) calculation<!>Penalty matrix definition: rule file<!><!>Problem: strained rings<!><!>Additional applications of bond groups<!>Advantages and limitations<!>Extended bond-charge increment scheme<!><!>Fitting considerations and results<!>Per-charge penalties<!>Benchmarking<!><!>Case study 1: 2-(dimethylcarbamoyl)indole (2)<!>Case study 2: N-benzylacetamide (3)<!>Summary<!>
<p>Intramolecular (internal, bonded terms)
∑bondsKb(b−b0)2+∑anglesKθ(θ−θ0)2+∑dihedralsKϕ(1+cos(nϕ−δ))+∑improper dihedralsKφ(φ−φ0)2+∑Urey−BradleyKUB(r1,3−r1,3;0)2
</p><p>Intermolecular (external, nonbonded terms)
(1)∑nonbondedqiqj4πDrij+εij[(Rmin,ijrij)12−2(Rmin,ijrij)6]
</p><!><p>The problem of assigning bonded parameters and charges to an organic molecule is far from new, and a wide variety of approaches have been applied in organic force fields such as the Tripos force field,6 CVFF,7 UFF,8 CFF93,9 Allinger's MM2,10 MM311–13 and MM414 force fields and Momany and Rone's commercial CHARMm force field15 (not to be confused with the academic CHARMM force field; note the different capitalization). Of particular note are the General Amber Force Field (GAFF)16 and its associated antechamber toolkit,17 which marked the first extension of a specialized biomolecular force field (ie. AMBER18, 19) to organic molecules, and MMFF94,20–22 which provided inspiration for the present charge assignment scheme. Accordingly, the present charge assignment scheme is empirical, as are most of the aforementioned charge assignment approaches. Notable exceptions are AMBER and OPLS, where charge assignment schemes based on a QM electron density (typically AM1 or PM3) have been proposed for small organic molecules.17, 23–27 While the resulting charges capture polarization and quantum effects and therefore are typically significantly more accurate than charges generated by more empirical assignment schemes,28–30 they are also computationally significantly more expensive, which may become a bottleneck in high-throughput applications involving millions of molecules.</p><!><p>To obtain missing parameters in CGenFF, existing parameters that were determined as part of the parametrization of CGenFF are identified by analogy based on the similarity between the atom types that define the parameters. Charges are assigned using a bond-charge increment scheme that is conceptually similar to the one implemented in MMFF94.20–22 The charge increments, which function as parameters in the charge assignment algorithm, are optimized to reproduce the charges on a training set of model compounds that were part of the parametrization of the force field. A noteworthy innovation in the presented charge assignment scheme is the fact that apart from the single charge increment associated with each possible bond (as defined by 2 atom types), there are 2 charge increments associated with each possible angle and 3 charge increments with each possible dihedral angle (as defined by 3 and 4 atom types, respectively). The combined functionality in the CGenFF atom typer1 and parameter assignment scheme makes it possible to apply the force field to an arbitrary molecule in an automated fashion. Notably, as a measure of the accuracy of the approximation, a "penalty score" is returned to the user for every bonded parameter and charge, which can be used to guide selective optimization of parameters and charges.</p><!><p>The assignment of missing parameters is based on penalty scores quantifying the dissimilarity between atom types. The penalty scores are based on an n × n penalty matrix, where n is the number of atom types. It should be noted that this matrix is not necessarily symmetric; substituting a very specialized atom type with a more generic variant typically has a lower penalty score than vice versa because it carries a lower risk of catastrophically inconsistent behavior. Before discussing the construction of the penalty matrix, we give a brief overview of the computational methods for assigning parameters by analogy, which provides a foundation for assignment of the penalty scores.</p><!><p>The total penalty score between two parameters (ie. a missing parameter and an existing parameter) can be calculated as the sum of the penalties for substituting the atom types defining the one parameter with the atom types in the second parameter. Thus, for every missing parameter, the program for assignment of parameters by analogy is tasked with finding the existing parameter in the parameter file (or parameter database) with the lowest TPS with respect to the missing parameter. This is accomplished by simply iterating through the parameter database, calculating the TPS for every parameter and returning the parameter with the lowest penalty. This process is sped up by the fact that both the penalty matrix and the data structure containing the parameters are collated by atom type, reducing the number and computational expense of lookup operations in the penalty matrix. The speed gain achieved by collating the parameters is a result of the fact that a parameter in a majority of cases will only differ from the previous parameter in the collated list by its last atom type(s). Thus, in the example of an angle parameter defined by atom types A-B-C, the algorithm can be made to keep track of the penalty of the first atom type (A) and the sum of the penalties of the first two atom types (A+B). If it is determined that only the last atom type is different (C'), only this atom type needs to be looked up in the penalty matrix and its penalty can be added to the latter sum. Similarly, if the two last atom types are different (B'–C"), only those two atom types must be looked up and their sum added to the previously determined penalty for A, thus substantially improving performance. Finally, for every existing parameter, the lowest of two TPSs should be used: the TPS between the "ordered" missing parameter A-B-C and the existing parameter, and the same TPS for the (physically equivalent) "reverse ordered" missing parameter C-B-A, as the former is not guaranteed to be always lower than the latter. All resulting parameters are stored in a binary search tree, which makes it easy to determine whether a missing parameter has been encountered and assigned before.</p><!><p>The penalty scores for atom substitutions ideally should reflect the chemical characteristics on which the definitions of the atom types were based. However, as the number of atom types in CGenFF is approximately 150 and is likely to increase with the parametrization of new chemical groups, directly populating the penalty matrix manually is not feasible, or would at the very least be extremely error-prone. Instead, the full matrix is constructed from a smaller number of submatrices and offsets, which are defined in a second section of the rule file used to assign atom types, as described in the preceding article.1 This "parameter assignment" section has a formally similar layout as the atom typing rules, although it is processed quite differently. Atom types are divided in a tree-like structure, with categories and subcategories, delineated by a "cat" keyword defining the name of the category and an "end keyword". Every entry in a category can either be an atom type ("typ") or a reference to a subcategory ("sub"), the name of which is followed by a colon. The difference with the atom typing rules lies in the keywords after the colon, the first of which is "pri" (short for "priority"), followed by the penalty offset added to the atom type substitution's total penalty score when entering the subcategory through a "sub" keyword in a hierarchically higher category (this will be clarified by an example in the next paragraph). This is followed by one or more "alt" keywords detailing alternative atom types or subcategories within the same category. The number of "alt" entries is equal to the number of rules in the current category minus one; every "alt" keyword is followed by the name of another ("typ" or "sub") rule in the same category and the penalty score for substituting the current type with the named alternative type. Finally, the rule ends with an "up" keyword, followed by the penalty offset for going one level up in the hierarchy.</p><p>As an example, Table 1 contains a short extract of the rule file that handles sp3 nitrogen atom types. In the category NG3, two subcategories are defined: NG3N and NG3P for neutral and positively charged sp3 nitrogen atom types, respectively. Consider for example the protonated primary ammonium type NG3P3; substituting this atom type with NG3P2 (secondary ammonium), NG3P1 (tertiary ammonium) and NG3P0 (quaternary ammonium) types respectively incurs penalties of 1, 2 and 4. If these substitutions fail, the penalty for going up in the hierarchy is 8. Doing so brings us to the sub NG3P rule in the NG3 category, which specifies that substituting an atom in the NG3P subcategory with an atom in the NG3N subcategory carries a penalty of 2. The penalties are added, so the total penalty score so far is 8+2=10. After entering the category NG3N, the rule with the lowest "pri" penalty is picked, which by convention is the first rule in the category, in this case NG321 with a penalty of 0. Thus, the final penalty score for substituting NG3P3 with NG321 is 8+2+0=10. Similarly, the penalty score for substituting NG3P3 with NG311 is 8+2+0.5=10.5, and so forth. This scheme allows filling the full penalty matrix using a rule file of manageable size; it is left to the writer of the rule file to make a trade-off between specifying more specific (and thus potentially more accurate) values by making fewer larger categories on the one hand, and having a smaller rule file (and thus less potential for "programming" error) by creating more smaller categories on the other hand. However, as the number of penalties in a category goes with the square of the number of rules in that category, making large categories quickly becomes unwieldy. Thus, the largest category in the parameter section of the current version of the CGenFF rule file contains 25 hydrogen atom types, and the penalty scores in this category were machine-generated based on reference bond lengths in the CGenFF parameter file. The second largest and largest manually generated category contains 10 sp3 oxygen atom types.</p><p>Most of the penalties in the rule file are chosen manually; it is a relative rare occurrence that an objective measure for similarity could be found that performed satisfactory for the purpose of assigning parameters by analogy. During the process of "programming" the penalty matrix, certain constraints were respected; for instance, values at a certain level in the tree are all in the same order of magnitude. This is illustrated in Table 1: the NG3 category has an "up" penalty of 12 for both "sub" rules, while both of its subcategories have "up" penalties of 8 for all their atom types, and penalties within each category are lower than this "up" value. Also, the penalty for substituting atom A for atom B was usually kept similar to the penalty for substituting atom B for atom A. However, these rules were not always followed rigidly. For example, in cases where a very specific type is changed into a type that is clearly more general, the penalty was deliberately set lower than for the reverse substitution (ie. changing a general type into a very specific one). In the end, what mattered more than any constraint in the construction of the penalty matrix was its ability to assign sensible parameters by analogy, and a large part of the effort was focused on validating their ability to assign parameters by analogy. This validation also gave rise to the below recommendations to validate and re-optimize parameters with TPS values above certain thresholds. Most of this work was based on chemical intuition, as it is not straightforward to do this in a systematic or quantitative way. Indeed, CGenFF internal parameters are validated against MP2 target data, which is relatively expensive to compute, and it is not trivial to devise a measure to objectively compare their performance when applied on different molecules. The present paper contains two case studies; a more systematic validation study is the subject of future work.</p><p>The whole rule tree is read into a temporary data structure and checked for consistency by a tree-walking algorithm that also replaces all atom type and rule name strings with pointers. Then, a different tree-walking algorithm is started at the first "typ" rule and exhaustively visits all other rules in the tree in a similar fashion as outlined in the above example, filling in one line of the penalty matrix; the aforementioned pointers help keep the computational cost of this process well within acceptable limits. The same procedure is repeated for all "typ" rules to complete the penalty matrix. The rows and columns of the penalty matrix are then respectively sorted by atom type and by penalty.</p><!><p>The penalties for substituting the central atom in an angle or improper dihedral and the two inner atoms in a dihedral are multiplied by 10 because changing these inner atoms has a much higher impact than changing the outer atoms. The same multiplication factor is applied to the penalties for substituting atoms in bond parameters because the impact of changing atoms in a bond parameter is typically comparable to changing the inner atoms in an angle or dihedral.</p><p>The penalty matrix used for these "inner atoms" (including the atoms in a bond parameter) will henceforward be called the "bonded matrix", in contrast to the "nonbonded matrix" that is used for the two "outer atoms" in angle and dihedral parameters. Specifically, the bonded matrix focuses more on bonded properties (valence and hybridization state) while the nonbonded matrix is more aimed at describing discrepancies in nonbonded properties (electrostatic and Van der Waals). For example, in the definition of the bonded matrix in the rule file, the highest category in the hierarchy differentiates subcategories by hybridization state, which are subsequently divided by period (ie. row in the periodic table), while in the nonbonded matrix, the period has precedence over the hybridization state. The use of these two different matrices can be rationalized by considering that the potential associated with a dihedral defined by atoms 1–2–3–4 is typically the result of a combination of nonbonded 1–4 interactions and the bond order of the 2–3 bond and the hybridization states of atoms 2 and 3, which are bonded properties. A similar observation can be made about angle potentials. However, this reasoning does not apply to improper dihedrals, where better results were obtained by using the bonded matrix for all 4 atoms (although the aforementioned scaling factor only applies to the central atom).</p><!><p>The algorithm as described above was made available for public beta testing as part of versions 0.9.0 and 0.9.1 of the CGenFF program. During subsequent testing, an important shortcoming was identified; see the 2-(dimethylcarbamoyl)indole (compound 2 in Figure 1) case study below. Consider a planar 5-membered ring as depicted in Figure 2; geometry dictates that the sum of the ring's inner angles α must be 540° and that the sum of the angles α+β+β around each of the 5 trigonal planar centers must be 360°. Indeed, in practice, pyrrole and other 5-membered aromatic heterocycles have inner angles α that are relatively close to the ideal 540° / 5 = 108° and corresponding outer angles β close to (360° − 108°) / 2 = 126°. In certain cases where the parameter assignment algorithm is given a substituted 5-membered ring in which one of the outer angles corresponds to a "missing parameter", such as the C10-C2-C3 angle in compound 2, it will use an in-ring parameter as a source of analogy (or vice versa), which leads to an error of (126° − 108°) = 18° in the reference angle and an unacceptable deformation in the geometry. Indeed, in compound 2, C10 is an sp2 carbon, which exhibits significant similarity to an sp2 carbon in a 5-membered ring from a nonbonded perspective. Therefore, the parameter assignment algorithm used the angle parameter defined by atom types CG2R51 CG2R51 CG2R51 as a source of analogy for CG2O1 CG2R51 CG2R51 (and CG2R51 CG2R51 NG2R51 for the CG2O1 CG2R51 NG2R51 parameter on the C10-C2-N1 angle), giving rise to the 19.7° deviation in Table 2 and to a pronounced deviation from planarity (φ(C1-C3-N1-C10) = 26.2°). Although the penalty score for this particular case was 26.5, which indicates that validation is required, similar cases were found with penalties lower than 10, the "recommended validation threshold". As planar 5-membered rings are ubiquitous in drug-like molecules, and this phenomenon can be shown to be even more pronounced in 3- and 4-membered rings, it seriously jeopardized the applicability of the parameter assignment algorithm. We first attempted to remedy this shortcoming by increasing the penalties for substitutions between atoms types for rings with different sizes and non-cyclic atom types, but it was found that a very large increase in the penalties was required in order to eliminate the problem, and that this large increase considerably affected the quality of other parameters by analogy for ring substitutions. For example, before increasing the penalties, parameters with a 6-membered ring atom type for C2 were used as a source of analogy for the C2-C10-N11 and C2-C10-O10 angles and were given a justifiably low penalty. Penalties that are sufficiently high to cure the problem discussed above would make this impossible.</p><!><p>Similar to the above, the bond group penalties for substituting the central bond in a dihedral parameter or any bond in an angle or bond parameter are multiplied by 10 because changing these bonds has a much higher impact than changing the outer bonds in a dihedral or the bonds in an improper.</p><p>The two first bond groups in the rule file are considered equivalent for the purpose of determining differences in bond group membership between two virtual bonds. This exception makes it possible to use the bond group feature to improve treatment of conjugated double bonds, as discussed below</p><!><p>In addition to fixing the original problem with strained rings, the bond group feature allows for a better treatment of biphenyls and conjugated double bonds. Specifically, as described in the preceding manuscript (part I),1 there are 2 sets of atom types for conjugated double bonds, CG2DC1 CG2D1O CG25C1 CG251O and CG2DC2 CG2D2O CG25C2 CG252O, which are applied such that in a chain of conjugated double bonds, double bonded atoms are given atom types from the same set and single-bonded atoms get atom types from a different set. Initially,* it was necessary to apply very high penalties to substitutions involving these atom types in order to prevent single-bond parameters being assigned to double bonds and vice versa. As a side effect of these high penalties, favorable substitutions involving conjugated double bond atom types were also suppressed. For example, a missing parameter containing the virtual non-conjugated double bond CG2D1–CG2D1 could never be substituted with its "conjugated equivalents" CG2DC1–CG2DC1 and CG2DC2–CG2DC2, because penalties low enough to allow this substitution would make it equally likely to substitute the original virtual double bond with the virtual single bond CG2DC1–CG2DC2. The introduction of the first two bond groups in Table 3 (which, as mentioned above, are considered equivalent for the purpose of determining differences in bond group membership) effectively gives the parameter assignment algorithm explicit knowledge of the bond orders. Consequently, the penalties for substituting conjugated double bond atom types could be set very low, instead relying on the bond groups to penalize incorrect substitutions. Analogous observations were made for biphenyls.</p><!><p>Finally, it should be noted that, compared to the "wildcards" used in a number of force fields, the current scheme offers a more exhaustive effort to identify the "best" parameter to be used when an exact match for a parameter is not present. It also offers the ability to generate penalty scores that may alert the user to possible limitations in the assigned parameters. Nevertheless, it should be emphasized that these penalty scores are based purely on the level of chemical similarity of the atom types in a parameter. Consequently, even in cases where there is a penalty of zero (ie. when the exact parameter is already available in CGenFF), the quality of that parameter for determining the physical properties of the molecule may be poor. For example, all the parameters for a given type of linker between two rings might be available in CGenFF. However, if the nature of the rings changes significantly without a corresponding change in the atom types across the ring (eg. the C-H carbon atoms in both 2-pyridone and the pyridinium ion carry the CG2R62 type, while the former is neutral and the latter is charged), then the conformational properties of that linker are likely to change such that the associated dihedral parameters, with penalties of zero, may require optimization. This speaks to the limited transferability of empirical force field parameters and the need for the user to make decisions on the level of accuracy they need for system under study. Conversely, cases have been encountered with large penalties yet satisfactory accuracy, which can be understood by interpreting the penalty as a measure of the statistical imprecision or "spread" of the parameter assignment.</p><!><p>Charges are assigned using a bond-charge increment scheme. From an atom-centric point of view, this can be formalized as qi=qi0−∑jβij where qi is the final partial charge on atom i, qi0 is the previously assigned formal charge,1 and βij is the bond charge increment for the bond between atoms i and j, with βji = − βij. From a bond-centric point of view, this translates to the following pseudocode:

    for each atom i { qi = qi0}
    for each bond n between atoms i and j {
       qi −= βij(n) ;
       qj += Î²ij(n) }

</p><p>We extended this scheme to angle charge increments α and dihedral charge increments δ:

    for each atom i { qi = qi0 }
    for each bond n between atoms i and j {
       qi −= βij(n) ;
       qj += Î²ij(n) }
    for each angle n between atoms i, j and k {
       qi −= αij(n) ;
       qj += αij(n) − αjk(n) ;
       qk += Î±jk(n) }
    for each dihedral n between atoms i, j, k and l {
       qi −= δij(n) ;
       qj += δij(n) − δjk(n) ;
       qk += δjk(n) − δkl(n) ;
       ql += Î´kl(n) }

</p><p>This is implemented by associating one charge increment βij with every bond in the parameter file, two charge increments αij and αjk with every angle and three charge increments δij, δjk and δjk with every dihedral. For every bond, angle or dihedral in the molecules, the charge increments associated with the corresponding parameter in the parameter file are used. If there exists no such parameter, charge increments from an analogous parameter are applied, using the same algorithm for finding the best analogy as described above, except that the "nonbonded matrix" is used for all atoms. As the algorithm for assigning parameters by analogy has the ability to match an existing parameter in which the atom types are in reverse order, care was taken to reverse the order and sign of the charge increments whenever such reversal occurs. It should be noted that the present algorithm only assigns correct symmetric charges to delocalized cases like amidinium, guanidinium and carboxylate when given symmetric formal charges qi0; therefore, it was necessary to include formal charge assignment in the atom typing rules, as discussed in the preceding manuscript (Part I).1 Finally, when applying this charge assignment algorithm to a test set of molecules that were not part of the CGenFF parametrization,1 it was found that charge increments associated with high-penalty dihedrals lead to a degradation in the quality of the predicted charges. Indeed, it is a known problem in force field development that transferability of parameters rapidly degrades with increasing number of atoms in the parameter's definition, and becomes problematic for dihedral parameters;31, 32 accordingly, this non-transferabilty was shown to apply to the associated charge increments. This problem was solved by ignoring charge increments associated with dihedral parameters with a penalty score higher than an arbitrary limit, which is currently set to 50.</p><!><p>First, all bond charge increments for bonds between different atom types are optimized (the charge increment for a bond between two identical atom types is zero by definition). In total, 384 degrees of freedom were fit in this stage, and the final RMS charge deviation was 0.0394 e. The resulting bond charge increments are rounded to the third digit after the decimal point and a number of bond charge increments involving hydrogen atoms are set to CHARMM standard values in order to make the hydrogen charges consistent with the CHARMM charge assignment rules. Specifically, aliphatic hydrogen atoms are always assigned a charge of +0.09 e, except when they are located on a 5-membered ring aliphatic carbon atom directly adjacent to a positively charged nitrogen atom, where they are assigned a standard charge of +0.28 e. Similarly, hydrogen atoms bound to aromatic, methylene sp2 (=CH2) and methine sp2 (=CH–) carbon atoms are given charges of +0.115, +0.21 and +0.15 e, respectively (see Table S1, Supporting Information). It should be noted that these changes amount to further rounding of the affected bond charge increments; since a vast majority of the molecules in the target data comply with the CHARMM charge assignment rules, the optimized bond charge increments were already at or very close to their ideal values before this adjustment, as testified by the 10−5
e increase in RMSD after adjustment.</p><p>Next, charge increments for angles are optimized. As with the bonds, charge increments on symmetric angles are constrained at zero. The total number of degrees of freedom at this stage was 2112, and the final RMSD was 0.0174 e. Again, the charge increments are rounded to the third digit after the decimal point. Charge increments smaller than 0.0025 e that involve a hydrogen are further rounded to zero, as well as the charge increments on a number of select angles (see Table S1, Supporting Information), in order to enforce the aforementioned CHARMM hydrogen charge assignment rules. This raised the RMSD by 4 × 10−5
e.</p><p>Finally, the optimization of the dihedral charge increments is performed, again constraining symmetric dihedrals to zero. Additionally, instead of setting a number of selected charge increments at 0 a posteriori as was done for the angles, 7 select dihedrals were constrained a priori to 0 (see Table S1, Supporting Information). Indeed, as this is the final optimization, there are no "higher-order" charge increments that can be used to compensate for the non-perfect fit caused by zeroing charge increments a posteriori. This left 6705 degrees of freedom to be fit, resulting in an RMSD of 0.0081 e. After rounding the charge increments to the third digit after the decimal point and further rounding charge increments smaller than 0.0025 e that involve a hydrogen to zero, the RMSD increased by 1.2 × 10−4
e ; consequently, the RMSD calculated using the final charge increments was 0.0082 e.†</p><!><p>The number of degrees of freedom in the final fit (6705) is lower than the number of target data points (9204) so that the fitting problem is formally not underdefined. However, there are many redundancies (ie. chemical groups that occur in more than 1 molecule) in the target data. Additionally, bond charge increments in a ring are only defined plus or minus a constant; for example, in a 3-membered ring A-B-C, adding a constant fractional charge simultaneously to the increments βab, βbc and βca does not change the resulting charge distribution. To counteract the resulting linear dependence in the fitting problem, all bond charge increments were subject to a restraining bias towards 0; in other words, the merit function χ2 of the least-squares fit is given by
χ2=∑ atoms i(qi−qit arget)2+k∑increments nβn2
where qi is calculated from the charge increments βn as described above and k is a factor determining the strenght of the restraining bias, which we set to 10−3. This bias has a minimal effect on properly determined charge increments, but effectively restrains in-ring increments and other poorly determined degrees of freedom to their lowest possible value. In this sense, the procedure is somewhat analogous to Bayly et al.'s RESP model.23</p><p>The restraining factor also effectively prevents overfitting of the aforementioned poorly determined degrees of freedom, and this preventive effect is further amplified by the incremental fitting and rounding scheme discussed in the subsection "Charge increment optimization" above. Indeed, of the 6705 degrees of freedom, 1536 fitted to absolute values of 0.002 e or lower, which usually do not have a significant effect on the charges. Thus, it can be argued that if the poorly determined degrees of freedom that were kept near zero by the restraining bias are not counted, 5169 well-defined degrees of freedom were fit to non-zero values using 9204 data points. This still does not fully eliminate the risk of overfitting in the sense that even a well-defined degree of freedom may reflect non-representative structures and charges in the target data. This risk is hard to quantify using the current data set because it is chemically very diverse2 and contains too little redundancy for proper cross-validation. Thus, the only way to identify and cure possible cases of "well-defined overfitting" is to increase the size of the target data set, which is an ongoing effort. Indeed, as the CHARMM charge assignment methodology is laborious,2 to perform a systematic validation study is a research project in its own right, and is the subject of future work.</p><p>Special treatment was required for conjugated double bond atom types. As mentioned under the subsection "Additional applications of bond groups" above, there are 2 sets of atom types for conjugated double bonds, which alternate in a chain of conjugated double bonds. This implies that every parameter involving these conjugated double bonds has a counterpart in which the atom types are switched. For the purpose of optimizing charge increments, extra constraints were implemented so that the same increments are applied to both counterparts, in order to prevent charge assignment from arbitrarily depending on the assignment of the double bonded atom types. Similarly, charge assignments associated with parameters that remain identical (except for reversal of the order of the atom types) after swapping double bonded atom types were constrained to zero.</p><p>Figure 3 shows the agreement between the 9681 manually assigned and optimized partial charges in CGenFF, and their counterparts that were determined using the fitted charge increments. When using only the bond charge increments (green squares), the data exhibits a relatively large spread, especially considering that a 0.1 e deviation (solid diagonal lines) in the partial charge of an atom that participates in a hydrogen bond typically results in a 1–3 kcal/mol deviation in the strength of that hydrogen bond. The spread becomes gradually less as the angle and dihedral charge increments are included, in line with the halving of the RMSD for every additional term as shown above. Our final approach, including bonds, angles and dihedrals (blue diamonds), is in excellent agreement with the existing charges. A few minor outliers are visible in Figure 3 in the form of horizontal chains of data points that cross the diagonal. The chain at y=0.115 extending to the right of the diagonal is caused by the fact that the aromatic hydrogen charges on difluorobenzene, difluorotoluene and benzylphosphonate were manually optimized and hence deviate from the standard CHARMM charges to which the charge increment optimization was constrained. A similar chain representing the corresponding carbon atoms can be seen at y=−0.115 extending to the left when only taking the bond charge increments into account (green squares). Unlike the hydrogen atoms, this chain disappears when including "higher-order" increments because no constraint was applied to the charges on non-hydrogen atoms; this suggests that our charge assignment scheme would indeed be able to capture mesomeric charge effects in aromatic rings if given an appropriate training set. Another horizontal line, at y=+0.197, consists of aromatic hydrogen atoms in indole-like rings. Indeed, the CGenFF collection of model compound includes three different charge sets for indole and a number of indole-like heterocycles that were all optimized independently and have different charges as a consequense of the common excess of degrees of freedom in the charge optimization for these kind of compounds. The 0.197 charge (as well as other charges on the indole ring system) is a least-squares compromise between the aforementioned model compounds, as testified by the chain of data points being horizontally centered around the diagonal and by the acceptable reproduction of QM water interaction energies demonstrated in the 2-(dimethylcarbamoyl)indole case study below. Just like for the regular atomatic rings, there is a less clearly resolved cluster of outliers representing the corresponding carbon atoms at the negative side of the plot.</p><!><p>Finally, a method was needed to derive a per-charge penalty σi from the per-parameter penalties of the parameters that contributed charge increments to a select atom. The first attempt was to use the maximum of all these contributions. This had the undesirable effect of producing high penalties for atoms that receive a small charge contribution from a high-penalty parameter (see for instance the M+∞ penalty for H2 and C2 in Table 4), making it clear that the magnitude of the contributing charge increments should influence the per-charge penalties. Next, we considered using averages of the per-parameter penalties σn, weighted by the absolute values of the corresponing charge increments βn (equation 2 with p = 1). This approach was intially chosen in favor of sums (see below) because we wanted the magnitude and interpretation of the per-charge penalties to be on the same scale as the per-parameter penalties. However, it was observed that a large number of small penalties with large weights could strongly pull down the average. This problem, which would give rise to unintuitive results such as the two atoms in a highly polarized bond having very different penalties (compare for instance M1 for C1 and C5 in Table 4), could be somewhat alleviated by using averages that are sensitive to high numbers, such as the quadratic, cubic, quartic and 8th-power means (ie. power means with exponent p = 2, 3, 4 and 8, respectively; see equation 2), but doing so comes at the cost of a decreased sensitivity to the weights, as can be rationalized by considering that in the limit for p → +∞, the power mean becomes the maximum.
(2)σi=Mp(σn)=∑increments n|βn|σnp∑increments n|βn|p
In the end, no acceptable trade-off between sensitivity to high penalties and sensitivity to the weights was found. Moreover, the very concept of averaging is not consistent with the cumulative nature of penalties; for a constant per-parameter penalty, applying a larger number of charge increments, or a charge increment that is larger in absolute value, should yield a larger per-charge penalty. In this sense, the nearly-identical M2–M8 penalties for N and H1 in Table 4 are undesirable. Therefore, we settled on using charge increment-weighted sums of the penalties. To retain the desired sensitivity to high penalties, we use the root of the sum of the squares, thus effectively treating the per-parameter penalties as uncorrelated relative errors of the charge increments that are summed to obtain partial charges (equation 3 with q = 1).
(3)σi=∑increments n|βn|qσn2
Weighting this sum by the charge increments yielded results that were deemed too sensitive to the weights. Compare for instance C1 and C4 in Table 4; while C1 should have a higher penalty because it involves more high-value-high-penalty increments, the factor ~16 difference in the q=1 column is too high. To decrease this sensitivity, we instead tried weighting by the square root (q = ½) and the cube root (q = ⅓) of the weight. Acceptable sensitivity was obtained with both options, but the cube root had the additional advantage of bringing the final sum in an acceptable range, eliminating the need to apply an arbitrary scaling factor to compensate for the fact that the sum of the absolute values of the charge increments is much lower than 1 in a vast majority of cases. Finally, the absolute values of all charge increments were increased by δ = (0.05)6 = 1.5625 × 10−8 before being used as weight factors. As can be seen in the final expression for σi, the per-charge penalty on atom i, this is equivalent to assigning a 5% weight factor to the per-parameter penalties σn associated with increments βn=0 (equation 4). Doing so is important in the presence of dihedrals with a high penalty, for which the charge increment is set to 0 and the penalty to the maximum allowable value for dihedrals (currently 50, as mentioned earlier). Although these increments do not contribute to the partial charge, they should have some effect on the penalty to reflect the fact that a contribution to the partial charge was ignored.
(4)σi=∑increments n|βn|+δ3σn2
</p><!><p>As a simple benchmark for the performance of the methods described in the preceding and present paper, we downloaded a ready-to-dock version (reference pH = 7) of the Maybridge HitFinder™ Collection from the ZINC database33 (http://zinc.docking.org/catalogs/maybhit/). The resulting mol2 file contained 18662 structures out of which 18429 could be successfully assigned atom types; the remaining 223 compounds contained chemical groups that are not yet supported by the CGenFF force field at the time of writing, most importantly 4-membered rings. Atom typing and assignment of charges and parameter by analogy on these 18429 structures took 619 seconds on a single core of a 3.1 MHz AMD Athlon II X2 255 processor, which can be considered a below-average desktop processor at the time of writing. This is equivalent to a rate of 29.8 molecules per second, demonstrating the suitability of the current method for high-throughput applications.</p><!><p>1 bond (penalty 40 both with and without the bond group feature discussed above under the heading "Assignment of parameters by analogy"‡)</p><p>4 angles (maximum penalty 26.5 without bond groups or 71 with bond groups)</p><p>10 dihedrals (maximum penalty 88 without bond groups or 128 with bond groups)</p><p>1 improper dihedral (penalty 5)</p><!><p>Manual parameter optimization based on geometric measurements of the QM and CGenFF minimized conformations readily allowed the N1-C2-C10 and C3-C2-C10 angles to be brought within 2° of the QM values (last columns in Table 2). Further properties reported in this paragraph include this correction so that the angular deviation does not interfere with the evaluation of the dihedrals. As for the interaction energies with water, which were calculated using the MP2 geometry (Table 5), C7, N1 and O10 are significantly outside the 0.2 kcal/mol difference criterion for CGenFF, and the interaction distance for C7 is longer than expected, in line with its interaction energy being too weak. This can be explained by the fact that the current training set for the charge increments is based on the assumption that substituted aromatic systems have the same charges as unsubstituted ones, thus negating the algorithm's ability to represent mesomeric effects. However, while allowing mesomeric substituent effects to influence the charges in the aromatic rings in the training set may improve N1, it will not necessarily improve C7 because the current implementation can only capture these electronic effects up to a distance of 3 bonds. Fortunately, the hydrogen bond donated by C7 is one of the weaker interactions in the system, and most relevant intermolecular properties are dominated by the stronger hydrogen bonds. In this sense, it is comforting that the hydrogen bond accepted by O10 is reproduced significantly better. Indeed, its interaction energy is only slightly outside of the CGenFF target range, which is an encouraging result given that C10 is involved in the formation of a novel bond, as indicated by its high penalty. Finally, although N1's penalty is in line with its deviation from the QM target data, the example of C7 illustrates that even a penalty of 0 does not guarantee perfect agreement. As discussed above, the penalty is a measure of the degree of analogy between chemical groups in the training set and the user's molecule, and neither factors in the quality of the training set, nor – more importantly – the fundamental appropriateness of applying charges and parameters by analogy for the given molecule. The directions of the HF and MP2 dipole moments differ by a relatively large 16°; the CGenFF dipole moment is in better agreement with the HF dipole moment (Table 6), which is a commonly observed consequence of the fact that the charges are primarily fitted to reproduce HF interaction energies. However, the agreement with both the HF and MP2 dipole moment is acceptable given the fact that the magnitude of the CGenFF dipole moment should be overestimated by 20 to 50%.2 The potential energy scan (PES) around the C2-C10 rotatable bond is shown in Figure 5. The local minimum at 30° is 0.87 kcal/mol too low and shifted by 15° in CGenFF, while the barriers at 0° and 90° are 0.84 kcal/mol too high and 1.61 kcal/mol too low, respectively. While the agreement is not as good as in the CGenFF results following parameter optimization, these PES are expected to yield qualitatively similar results in MD simulations, which is very good considering that the penalties of the 4 dihedrals terms around this bond are 87, 88, 127 and 128, and our guidelines mandate extensive validation and/or re-optimization for any parameters with penalties higher than 50.</p><!><p>A second case study is N-benzylacetamide (compound 3 in Figure 1, Table S4 of the supporting information). In this case, one angle with penalty 28.5 was assigned by analogy, as well as three dihedrals with penalties 47, 28.5 and 28.5. The highest penalty of the charges is 39.4; a comparison between CGenFF charges and standard CHARMM charges is given in Table 7. This table shows a fair correlation between the deviation from the standard CHARMM charges and the penalty; the three highest penalties correspond to the three charges that have significant (though less than 0.05 e) deviations. Given the lack of high-penalty bond and angle parameters, it is not surprise that all minimized bond lengths and angles are in excellent agreement (see Table S5 of the supporting information). Figure 6 depicts the relaxed potential energy surfaces around the C5-C6-C7-N8 and C6-C7-N8-C9 dihedrals. The parameters by analogy for C6-C7-N8-C9 produce good agreement with the QM from 0° to ~120°, but the PES diverge by 2 kcal/mol from 120° to 180°. Based on the geometry at 180° (Figure 7), we speculate that this is because the intramolecular electrostatic interaction of the amide NH group with the benzene ring is substantially different from its interaction with an ester or protonated carboxylic acid in a C-terminal amino acid, which is the origin of the parameter by analogy for this dihedral. Manual fitting of this dihedral leads to the introduction of a dihedral term with n=1, Kφ=1.7 and δ=180 and a second dihedral term with n=3, Kφ=1.7 and δ=0. The former term's relatively high Kφ might indicate that the intramolecular nonbonded interactions are not perfectly balanced for the purpose of reproducing this vacuum MP2 PES, and we speculate that the fitted parameters might be poorly transferable to other molecules. This is in line with the observation that the transferability of dihedral parameters in general is often problematic; see section "charge assignment" above. Although the penalty of 28.5 associated with this parameter mandates validation, it is possible in some cases to obtain a similarly mediocre agreement for parameters with even lower penalties; therefore, if accuracy is of the utmost importance in the molecule(s) under study, validation is recommended regardless of the penalties. Concerning the C5-C6-C7-N8 scan, the data points with poor agreement between CGenFF and MP2 are caused by the C6-C7-N8-C9 dihedral (which is relaxed at every scan point scan) being close to 180°. Indeed, the value for the the C5-C6-C7-N8 parameter did not change significantly during manual optimization, and the improvement in agreement is purely due to the optimization of the C6-C7-N8-C9 dihedral. This demonstrates that a higher penalty is not a guarantee for poor agreement; rather, the penalty scores should be thought of as measure for the approximation error on the parameter.</p><!><p>The present paper describes the algorithms for automatic assignment of bonded parameters and charges in CGenFF, thereby allowing for automatic application of CGenFF to arbitrary molecules. Bonded parameters are assigned by substituting atom types in the definition of the desired parameter. The assignment of bonded parameters also forms the basis for the charge assignment. A penalty is associated with every substitution and the existing parameter with the lowest sum of penalties is chosen as an approximation for the desired parameter; the aforementioned sum is returned to the user as a measure for the accuracy of the approximation. Charges are assigned using a variation on the classical bond-charge increment scheme, which essentially associates a bond charge increment with each bond parameter. We extended this scheme by associating two charge increments with each angle parameter and three charge increments with each dihedral parameter. To apply this bond-charge increment scheme to a novel molecule, the aforementioned algorithm for assigning parameters by analogy is used, albeit using slightly different criteria for analogy. This allows for per-charge penalty scores to be returned, which give the user an idea of the quality of the charges. Charge increment values are optimized to reproduce the charges on model compounds that were part of the parametrization of the force field as well as possible using a restrained least-squares fitting algorithm. A final charge RMSD of 0.0082 e was obtained. The validity of the final parameters and charges and the relevance of the associated penalty scores were illustrated by 2 case studies; a more thorough validation study will be the subject of a future paper. As the methods discussed in the preceding and present paper are entirely empirical, they are excellently suited for high-throughput applications, as testified by the fact that we processed a representative database of 18429 drug-like structures at a rate of 29.8 compounds per second on a single core of a below-average desktop processor. Finally, it should be reiterated that the penalties are indicative of the analogy of the generated parameters or charges with available terms in CGenFF, rather than the intrinsic fitness of these terms for application in the molecule of interest. Consequently, large penalties may yield satisfactory accuracy, while significant deviations from acceptable accuracy may occur when small or zero penalties are reported. Indeed, given the sensitivity of chemical properties to environment, in cases where parameters are available in CGenFF for a new molecule, there may still be significant limitations in the accuracy of the resulting model. It is important that the user is aware of such limitations and perform the appropriate validation and optimization as required.</p><p>Use of the CGenFF program for automatic atom typing and assignment of parameters and charges by analogy is provided at https://www.paramchem.org/.</p><!><p>In versions 0.9.0 and 0.9.1.</p><p>All numbers presented in this paragraph are based on release 0.9.0 of the CGenFF program (force field version 2b6). As the charge increment fitting is based on explicitly parameterized compounds for which all parameters are available, the algorithmic improvements discussed in this section do not influence the results. The release of version 2b7 of the force field does, but not to a significant extent.</p><p>The assignment without and with the "bond group" feature was respectively performed with version 0.9.1 and 0.9.6 of the CGenFF program.</p><p>
          Supporting Information
</p><p>A list of bond, angle and dihedral charge increments that were constrained during the charge increment optimization.</p><p>Toppar stream files for the case studies as generated by the CGenFF program.</p><p>A table with minimized bond length and angles for N-benzylacetamide (3), analogous to Table 2.</p><p>This material is available free of charge via the Internet at http://pubs.acs.org.</p>