# Increasing the Accuracy and Robustness of the CHARMM General Force Field with an Expanded Training Set

Anastasia Croitoru, Anmol Kumar, Jean-Christophe Lambry, Jihyeon Lee, Suliman Sharif, Wenbo Yu, Alexander Mackerell, Alexey Aleksandrov

<figure><img src="image_1.png" /></figure>

## To cite this version:

Anastasia Croitoru, Anmol Kumar, Jean-Christophe Lambry, Jihyeon Lee, Suliman Sharif, et al. Increasing the Accuracy and Robustness of the CHARMM General Force Field with an Expanded Training Set. Journal of Chemical Theory and Computation, 2025, 21 (6), pp.3044-3065. 10.1021/acs.jctc.5c00046 . hal-05231720

HAL Id: hal-05231720

https://hal.science/hal-05231720v1

Submitted on 31 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Increasing the Accuracy and Robustness of the CHARMM General
Force Field with an Expanded Training Set

Anastasia Croitoru,¹,² Anmol Kumar,² Jean-Christophe Lambry,¹ Jihyeon Lee,² Suliman Sharif,² Wenbo Yu,² Alexander D. MacKerell Jr.,²* and Alexey Aleksandrov¹*

¹Laboratoire d'Optique et Biosciences (CNRS UMR7645, INSERM U1182), Ecole Polytechnique, Institut polytechnique de Paris, F-91128 Palaiseau, France

²Department of Pharmaceutical Sciences, School of Pharmacy, University of Maryland, 20 Penn Street, Baltimore, Maryland 21201, USA

*Corresponding authors: [alex@outerbanks.umaryland.edu](mailto:alex@outerbanks.umaryland.edu) and [alexey.aleksandrov@polytechnique.edu](mailto:alexey.aleksandrov@polytechnique.edu)

Keywords: CHARMM; CGenFF; force field development; medicinal chemistry; drug-like molecules; ZINC library
ABSTRACT

Small molecule empirical force fields (FFs), including the CHARMM General Force Field (CGenFF), are designed to have wide coverage of organic molecules and to rapidly assign parameters to molecules not explicitly included in the FF. Assignment of parameters to new molecules in CGenFF is based on a trained bond-angle-dihedral charge increment linear interpolation scheme for the partial atomic charges along with bonded parameters assigned based on analogy using a rules-based penalty score scheme associated with atom types and chemical connectivity. Accordingly, the accuracy of CGenFF is related to the extent of the training set of available parameters. In the present study that training set is extended by 1,390 molecules selected to represent connectivities new to CGenFF training compounds. Quantum mechanical (QM) data for optimized geometries, bond, valence angle, and dihedral angle potential energy scans, interactions with water, molecular dipole moments, and electrostatic potentials were used as target data. The resultant bonded parameters and partial atomic charges were used to train a new version of the CGenFF program, v5.0, which was used to generate parameters for a validation set of molecules, including drug-like molecules approved by the FDA, which were then benchmarked against both experimental and QM data. CGenFF v5.0 shows overall improvements with respect to QM intramolecular geometries, vibrations, dihedral potential energy scans, dipole moments and interactions with water. Tests of pure solvent properties of 216 molecules show small improvements versus the previous release of CGenFF v2.5.1 reflecting the high quality of the Lennard-Jones parameters that were explicitly optimized during the initial optimization of both the CGenFF and the CHARMM36 force field. CGenFF v5.0 represents an improvement that is anticipated to more accurately model intramolecular geometries and strain energies as well as non-covalent interactions of drug-like and other organic molecules.
INTRODUCTION

Empirical force fields (FFs) represent the cornerstone of condensed phase molecular simulations.¹ These include studies of small chemical species, such as ionic liquids, biological macromolecules including highly heterogeneous assemblies, as well as complex condensates, among many others. While FFs for biological molecules are well established, with example FFs including the additive models CHARMM,²⁻⁵ AMBER,⁶⁻⁸ OPLS,⁹⁻¹¹ and GROMOS¹²⁻¹⁴ along with polarizable FFs such as the AMOEBA¹⁵⁻¹⁷ and Drude¹⁸⁻²² FFs, modeling the diverse range of molecules found in medicinal chemistry, ionic liquids, and polymers among others remains a significant challenge due to the vast chemical diversity of these classes of molecules. To meet this challenge small molecule force fields such as GAFF,⁷,²³ CGenFF,²⁴⁻²⁶ MMFF,²⁷,²⁸ and OPLS⁹,²⁹ were developed largely covering compounds associated with medicinal chemistry as well as other organic molecules. These FFs, whose origins date back 20 years or more, focused on the appropriate coverage and systematic parametrization in the context of computationally accessible models. Such FFs as well as more recent variants such as OpenFF,³⁰,³¹ were largely optimized targeting quantum mechanical (QM) data for the bonded terms and the partial atomic charges. Lennard-Jones (LJ) parameter optimization in GAFF, CGenFF, and OPLS primarily targeted pure solvent and crystal experimental data supplemented by QM data in certain cases. These efforts yielded models that have been applied to a wide range of systems allowing for calculations with sufficient accuracy on large systems (>100,000 atoms) and long time scales (> 1 µs) taking advantage of the simplicity of the underlying potential energy function. While more recently machine learning (ML) models have been developed and have shown success especially in the context of intramolecular energies and geometries,³²⁻³⁴ these models are significantly slower than empirical FFs. Furthermore, these models have limitations in the treatment of long-range interactions and charged species, and they lack transparency regarding the contributions of individual atoms or functional groups to calculated energies and forces.³⁵,³⁶ Thus, empirical FFs still remain the primary tools of computational studies of condensed phase systems associated with their computational speed, accurate treatment of condensed phases, through accurate treatment of long-range nonbond interactions, and are readily accessible to a wide chemical space including charged species. Additionally, empirical FFs benefit from decades of validation, with extensive comparisons to diverse experimental data that both define their limitations and underscore their reliability.

Beyond the computational requirements associated with the molecular simulations themselves are the costs associated with applying a given model to a new molecule and improving the accuracy of the FF itself. For example, with GAFF and related methods, the parametrization of the bonded portion of the model involves primarily a rules-based approach that is essentially instantaneous while the electrostatic parameters require a QM calculation that takes on the order of minutes. Recent ML models based on deep neural networks (DNNs) allow for instantaneous assignment of parameters, where the parameters are the weights in the trained DNN.³²⁻³⁴ Similarly, the CGenFF parameter assignment protocol, incorporating an ML model based on a linear interpolation scheme for the electrostatic parameters,²⁴ also generates both the bonded and electrostatic parameters instantaneously making these methods of utility for high-throughput computational studies common in computer-
aided drug design (CADD), as it eliminates the need to store an underlying model for each molecule. With respect to improving the accuracy and coverage of FFs, this requires an extension of the training set of molecules, addition of new atom types for novel topologies, and optimization of the model to accurately represent these additional molecules. If new atom types are required, it may be necessary to optimize the associated van der Waals parameters, typically using the LJ formalism, which represents the most challenging aspect of FF optimization. In contrast, ML models only require expansion of the training set.³⁷ However, while intuitively straightforward, for DNN-based ML models this requires generating significant numbers of conformations for each new molecule included in the training set due to the discrete nature of the parameters in the DNN. Physics-based FFs, on the other hand, require a relatively small number of conformations for each new molecule in the training set, because their physically grounded potential energy functions inherently allow accurate interpolation between regions of configurational space explicitly sampled in the model compounds. This inherent “feature” of a physics-based FF also enhances its transferability. The physical model captures the correct underlying molecular behavior (e.g., the harmonic nature of bond stretching or valence angle bending) through the simple form of the energy function, which is based on harmonic and sinusoidal terms. Furthermore, the inclusion of distance-dependent nonbonded interactions allows for effective representation of long-range and ionic effects.

Force field methods have proven highly effective in medicinal chemistry, enabling accurate modeling of protein-ligand interactions, prediction of binding poses, and calculation of binding free energies. Among these, the CHARMM/CGenFF force field has been widely incorporated into drug discovery protocols and tools, facilitating the discovery and optimization of ligands.³⁸⁻⁴⁶ To enhance its accessibility, CGenFF is integrated into user-friendly web-based platforms such as the CGenFF program and its web interface, and CHARMM-GUI tools,⁴⁷ enabling researchers to efficiently parameterize and simulate drug-like molecules. These protocols have been successfully applied in drug design projects, leading to the discovery of novel molecules, including ligands for bromodomains,⁴⁸⁻⁵¹ epitranscriptomic targets such as the m7G-RNA writer methyltransferase⁵² and the m6A-RNA reader YTHDC1,⁵³ as well as inhibitors of the serine protease TMPRSS2, a promising antiviral target.⁵⁴ These achievements demonstrate the growing role of CGenFF in drug discovery while also emphasizing the need to further expand and refine the CGenFF force field to encompass a broader spectrum of drug-like molecules, including those with complex chemical scaffolds containing non-canonical ring systems, such as spirocycles and bridged rings.

CGenFF is a physics-based force field model in which bonded parameters are assigned based on analogy, while the charges are assigned with an ML model using a linear interpolation scheme that is also assigned based on analogy. When a specific bonded parameter is absent in the CGenFF force field, the most analogous parameter is identified using a search guided by atom-type based penalty scores. The penalty scores, described in the original CGenFF paper,²⁶ involve a set of rules in which the penalty increases as the chemical nature of alternate atom types deviates from that of the requested atom type. This scheme yields the most analogous parameter available in the CGenFF parameters as well as assigns a penalty to the new parameter that the user may use as guidance concerning the need for optimization of parameters.⁵⁵ This by-analogy approach is also applied to determine the
bond, angle and dihedral parameters from which the increments for the ML charge assignment model are also obtained. As with the bonded terms, non-zero penalties are assigned to the charges for which exact connectivity matches in the force field are not available (i.e., penalties are zero for all bonded and charge parameters that have exact matches in the training set). Thus, the quality of the parameters assigned to a new molecule by the CGenFF program depends on the coverage of the CGenFF parameter set itself. However, since the penalties are based on analogy and not directly related to chemical accuracy, in many cases transferred bonded parameters or charges with high penalties are often sufficiently accurate for molecular simulations of new molecules, though testing of high-scoring parameters is warranted.⁵⁵ To maximize the accuracy of CGenFF based on the assumption that parameters explicitly optimized for a given model compound perform satisfactorily in previously unseen molecules with similar connectivities, it is necessary to expand the CGenFF training set. Such expansion will also lead to an overall decrease in the penalty scores assigned to the parameters and charges in new molecules.

In the present study the training set for CGenFF is extended with respect to both the bonded parameters and the partial atomic charges. This involves increasing the number of molecules in the training set from ~930 to over 2,300 molecules. The new molecules were selected from the ZINC20 database based on their atomic connectivities not present in the current CGenFF training set, their relatively small size (molecular weight 150 to 250 Daltons) and, with six exceptions, having a maximum of one rotatable bond between non-hydrogen atoms. Consistent with the development of CGenFF, parameters with limited analogy in the current training set are explicitly optimized within the selected molecules, and the final optimized parameters are used to train CGenFF, resulting in version 5.0. As presented below, this involves refitting the terms in the charge ML linear interpolation scheme based on bond, valence angle, and dihedral angle charge increments used to predict the charges. The extended set of bonded parameters directly supplements the existing collection, with parameters that are either transferred to a new molecule when an exact match is found or assigned by analogy using available parameters and employing the CGenFF penalty scheme.

In addition to expanding the training set with a collection of molecules selected from the ZINC20 database, this study also extends CGenFF to incorporate fused 4-membered rings. This extension involved adding new atom-typing rules to incorporate these molecules, followed by optimization of the corresponding bonded parameters for newly included molecules. All the additional molecules and parameters were then used to retrain the CGenFF ML charging engine, yielding version 5.0, which was then applied to the new training set molecules added to the force field as well as additional validation molecules including drug-like molecules approved by the US Food and Drug Administration (FDA). Results were compared to both experimental and quantum mechanical data. The resulting version 5.0 of CGenFF represents a significant improvement in the accuracy of the force field as well as an increase in its coverage, which are expected to improve the utility of CGenFF in CADD and other areas of computational chemistry and biophysics.

# METHODS AND MATERIALS
CHARMM potential energy function

The potential energy function of the additive all-atom CHARMM force field, as recently detailed in our work on nonstandard amino acids,⁵⁶ is that used for the CGenFF force field. The CHARMM additive potential energy is:

$$U = U_{inter} + U_{intra} \quad [Eq 1]$$

The intermolecular or non-bonded energy, $U_{inter}$, is due to electrostatic and van der Waals (vdW) interactions:

$$U_{inter} = \sum_{elec}^{nonbonded} \frac{q_i q_j}{4\pi\epsilon_0 r_{ij}} + \sum_{vdW}^{nonbonded} \epsilon_{ij} \left( \left( \frac{R_{min,ij}}{r_{ij}} \right)^{12} - 2 \left( \frac{R_{min,ij}}{r_{ij}} \right)^{6} \right) \quad [Eq 2]$$

The electrostatic term is described by Coulomb's law with $q_i$ and $q_j$ being the respective partial atomic charges on atoms i and j, and $r_{ij}$ is the distance between atoms i and j. The vdW term is treated by the LJ 6-12 potential in which $\epsilon_{ij}$ is the well depth, $R_{min,ij}$ is the radius at which the LJ potential has a minimum. In the additive CHARMM force field, the LJ parameters for pairs of atoms i and j are calculated using the Lorentz-Berthelot combination rule:⁵⁷

$$\epsilon_{ij} = \sqrt{\epsilon_i \epsilon_j} \quad \text{and} \quad R_{min,ij} = \frac{R_i + R_j}{2} \quad [Eq 3]$$

The intramolecular or bonded part of the potential energy function in Equation 1, $U_{intra}$, contains contributions for the bond, valence angle, dihedral angle, improper angle, and Urey-Bradley terms. It is given by:

$$U_{intra} = \sum_{bonds} K_b (b - b_0)^2 + \sum_{angles} K_{\theta} (\theta - \theta_0)^2 + \sum_{Urey-Bradley} K_{UB} (r_{1-3} - r_{1-3;0})^2 + \sum_{dihedrals} \sum_{n=1}^{N} K_n (1 + \cos(n\varphi - \delta_n)) + \sum_{impropers} K_{\varphi} (\varphi - \varphi_0)^2$$

, [Eq 4]

where $b_0$, $\theta_0$, $r_{1-3;0}$, and $\varphi_0$ are the bond, valence angle, Urey-Bradley, and improper angle equilibrium values, respectively; the K's are the corresponding force constants; and n and $\delta_n$ are the dihedral multiplicity and phase. The dihedral term is expressed through a Fourier series with N multiplicities. In CGenFF, the use of multiplicities is restricted to fewer than seven multiplicities (N < 7), typically excluding multiplicity 5. The phase of the dihedral
term typically assumes one of two possible values: 0° or 180° as this allows the same parameters to be applied to diastereoisomers around a chiral center. While the contribution of the UB term can enhance the accuracy of fitting vibrational properties, its parametrization poses a greater challenge due to the necessity of accounting for the cross-dependency between bond stretching and valence angle bending. Consequently, the UB term was omitted in this work. Additionally, the improper dihedral term was applied selectively, only in cases requiring special treatment for out-of-plane distortions.

## Model compounds

Compounds for expansion of the training set were selected from the ZINC20 database⁵⁸ based on their drug-like properties and their ability to cover a broad region of chemical space represented in CGenFF but not already included in its existing training set. The atom names from the ZINC20 library were retained for all selected molecules. These compounds were chosen to include common ring systems found in drug-like molecules, with the following criteria: they are not radicals, can be purchased, have molecular weights between 150 and 250 Daltons, typically contain no more than one rotatable bond between non-hydrogen atoms, and are composed only of the atoms C, N, O, H, S, and F. The selection resulted in 285,041 molecules. From this set, a representative subset was chosen as follows. First, initial parameters for each molecule were generated with the CGenFF program. Secondly, compounds were selected to cover the full chemical space of the 285,041 molecules, based on bonded parameters. In other words, all parameters required for the 285,041 molecules were included in the selected subset.

For bonded term parametrization, molecules containing two individual ring systems were split into their individual ring components (e.g., mono- or multicyclic fused rings), which we refer to as “fragments,” as described in the next paragraph. The set was further supplemented by neutral compounds, generated from initially ionized ZINC20 molecules by adding or removing hydrogen atoms to render them neutral. A total of 890 neutral compounds were added resulting in a final set of 7,387 molecules that covered all the bonded parameters present in the original 285,041 compounds.

To simplify the parameter optimization process “parent” compounds with a rotatable bond between two ring systems were dissected into two individual fragments, each containing individual or fused rings. During the construction of these fragments, a methyl group was added at the severed bond. In compounds where the rotatable bond included a nitrogen atom in a ring, only the fragment with the nitrogen was retained as the second fragment was typically already available. This fragment-based approach offers several significant advantages. It reduced the number of atoms involved in QM calculations, making them more computationally efficient. Additionally, it minimized the number of free parameters in each optimization and reduced the influence of multiple functional groups on the target QM data. For charge fitting, this methodology enables the examination of atom types in the context of various functional groups in individual molecules thereby avoiding secondary interactions in water-model compounds interactions. This process yielded 480 fragments following exclusion of fragments that generated errors. Errors included cases where CGenFF assigned an incorrect total charge, an issue to be addressed in future revisions of CGenFF, or where QM geometry optimization failed to converge. Once the individual-
fragment parameters were optimized, the new parameters for connectivity between the fragments were optimized in the complete parent compound. For the parent compounds the target QM data were computed on the complete molecules. Some selected compounds contained multiple ionizable groups within the same molecule, such as carboxylic acid and amine moieties. In these cases, QM gas-phase minimization could lead to proton transfer, resulting in tautomeric states that are less likely to occur under aqueous conditions. To address this, constraints were applied to preserve proton positions consistent with the protonation states specified in the ZINC20 database. Protonation states in ZINC20 were assigned using the ChemAxon tool,⁵⁹ which determines biologically relevant tautomers at physiological pH. Final compounds were selected based on bonded parameter-based sorting, using the penalty for each parameter. The penalty was zero for all parameters explicitly present in the initial CGenFF training set. For the parameters with penalties, only those with penalty >10 were considered for optimization. This process resulted in 180,992 bonded terms with high penalties that required optimization. From this list, the unique parameters were identified, and compounds containing those parameters were selected.

In CGenFF, bonded parameters are optimized in a hierarchical approach, where once a parameter is optimized, it is fixed and reused in subsequent compounds without further refinement. Thus, it was necessary to establish an optimization sequence for the compounds. Highest priority was given to neutral molecules, ranked by the number of unique parameters requiring optimization. For compounds with the same number of parameters, the compounds with the smaller number of atoms had higher priority. As the selection process progressed, compounds without unique parameters requiring optimization were removed. This process was repeated in an iterative fashion yielding a final ranked set of the 1,390 molecules for bonded parameter optimization. Before optimizing bonded parameters, all 1,390 compounds were subjected to optimization of the partial atomic charges. In CGenFF, residue names for these molecules follow a six-letter code, with the first two letters indicating the year of parametrization and the rest identifying the molecule. The SMILES strings of these compounds are publicly available on the MacKerell lab GitHub repository (https://github.com/mackerell-lab/cgenff5).

## Parametrization Protocol

The parametrization protocol is shown in Figure 1. For a given model compound, the CGenFF program was employed to assign CGenFF atom types²⁶ and to obtain initial guesses of the partial atomic charges and bonded parameters.²⁴, ²⁵ The partial charges that were assigned a zero penalty by the CGenFF program, i.e. already optimized in CGenFF, were not considered for optimization in the present study, except for selected zero-penalty atoms covalently bonded to high-penalty atoms. As stated above, only bonded parameters with penalties >10 were optimized. The LJ parameters were taken directly from CGenFF and were not further optimized in this work.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The force field parametrization workflow employed in this study. Initial partial charges and bonded parameters were generated with CGenFF 2.5.1 program and iteratively modified until MM agrees with QM minimum energy geometry, water interaction energies and distances, dipole moment, ESP, conformational energies, and PES scans. For anionic species, the 6-31+G(d) basis set was used for all QM computations.</figcaption></figure>

The initial geometries for the model compounds were taken from the ZINC20 library. The fragments described above were constructed using RDKit⁶⁰ or manually from parent molecules using the Avogadro program;⁶¹ missing hydrogens were added using the Open Babel software.⁶² For flexible molecules, conformational analysis was performed as described below. The lowest energy conformations identified from this analysis were then used to initiate the parametrization workflow. Optimized QM structures were used to optimize the atomic charges. The MM model with optimized charges was then used to optimize the bonded terms. Adiabatic QM potential energy surface (PES) scans were performed along the degrees of freedom for which bonded parameters were optimized. These bonded parameters were optimized to minimize differences between QM and MM geometries and the PES. All optimizations of FF parameters were performed with the program *developer* that was previously used to parametrize a large number of modified nucleotides,⁶³ non-standard amino acids,⁶⁴ and a number of important ligands in our previous studies.⁶³

## Conformational analysis

Conformational analysis employed the Conformer-Rotamer Ensemble Sampling Tool (CREST) version 2.12.⁶⁵ This tool utilizes semi-empirical tight-binding methods in conjunction with a meta-dynamics driven search algorithm. Subsequently, the resulting CREST conformations were initially optimized using the Psi4 program,⁶⁶ with further optimization performed with Gaussian 16.⁶⁷ The compound structures were optimized with the MP2/6-31G(d) model chemistry,⁶⁸,⁶⁹ or MP2/6-31+G(d) model chemistry for anionic molecules. All QM optimizations were performed to default tight tolerances. The resultant conformations were filtered based on their
relative energies and root mean square (RMS) deviations between coordinates. Conformations were considered identical if the RMS deviation (RMSD) for non-hydrogen atoms was less than 0.1 Å, and the absolute energy difference was also less than 0.1 kcal·mol⁻¹. The RMS deviation between QM-optimized conformations was computed using the obrms tool of the Open Babel.⁶²

## Optimization of partial atomic charges

Consistent with the development of the CHARMM force field, atomic charges were optimized targeting interactions between the model compound and individual water molecules and the dipole moment of the model compound. Additionally, QM electrostatic potentials (ESPs) were included as additional target data in the charge fitting as described in our previous study.⁷⁰ The weights of the dipole moments and ESPs were smaller than that used for water interactions (see below). Gaussian⁰⁹⁷¹ was used for the QM calculations of water interactions, ESPs, and dipole moments.

Atoms of the model compound, excluding carbons, were probed by individual water molecules placed in idealized linear orientations.⁷² For the MM calculations the CHARMM TIP3P model was used.⁷³ Different orientations of the water molecule were considered around the interaction axis: the complex was calculated every 45° or 90° of the water probe rotation for polar atoms, and one or two orientations for non-polar atoms. Example configurations of model compound-water interactions are shown in Figure 2. Each water-model compound dimer was optimized only by varying the interaction distance to find the energy minimum. The QM-optimized gas-phase geometry was used for the model compounds and TIP3P model geometry for the water molecule. The angle defining the orientation of the water molecule around the interacting axis was held fixed during the optimization. QM calculations were done at the HF/6-31G(d) level.²⁶,⁷² Following the CHARMM standard protocol the QM interaction energies were scaled (made more favorable) by an empirical factor of 1.16 for neutral polar compounds and the HF/6-31G(d) minimum interaction distance was corrected by subtracting 0.2 Å for all polar interactions involving neutral compounds.⁷² In the case of sulfur atoms, the model compound-water interactions were calculated at the MP2/6-31G(d) level including the basis set superposition error (BSSE) correction of Boys and Bernardi⁷⁴ and without applying the scaling and offset rules. Unfavorable QM water interactions were not considered in the fitting process. The molecular dipole moment, defined by the charge distribution in the MM model, was used to provide additional target data for the optimization of the atomic charges of neutral compounds only.⁶³ The dipole moment was calculated in vacuum at the MP2/6-31G(d) model chemistry using the QM-optimized conformation.²⁶ Following the standard CHARMM protocol, to implicitly account for the induced electronic polarization in an aqueous environment, the reference dipole moment being targeted in the MM optimization was increased by 30% with respect to the QM value.²⁶ Both the magnitude and direction of the dipole moment were targeted.⁶³

QM water interaction data may not be sufficient to define partial atomic charges on all atoms for large compounds, since only water interactions with a subset of atoms at the molecular surface are probed. To supplement the interaction data, ESP calculations were performed at the MP2/6-31G(d) model chemistry and at
MP2/6-31+G(d) for anions.⁷⁰ At each iteration during the charge optimization, the RMSD between QM and MM ESPs was evaluated and added with the corresponding weight to the target function. However, the weight for the ESPs (1.0 kcal⁻¹·mol·e⁻¹) was kept small relative to the weight for the water interactions (10.0 kcal⁻¹·mol). The dipole moment contributions were weighted (3.0 Debye⁻¹) more than the ESPs but less than the water interactions, as the reproduction of water-compound interaction energies and geometries is important to balance the solvent-solvent, solvent-solute, and solute-solute interactions.

The target function for fitting of the partial atomic charges included the following terms: the RMSD between MM and QM minimum probe-water interaction energies, the RMSD between QM and MM probe-water minimum interaction distances, the norm of the difference between the MM and QM dipole moments, the RMSD between QM and MM ESPs, and a restraint term to limit deviations from the initial charges assigned by CGenFF. The magnitude of the charge restraint term was scaled based on the CGenFF penalty for each atom. The calculation of this term is defined by:

$$F_{restr} = \sqrt{\frac{\sum_i W_i \cdot (q_i - q_i^{initial})^2}{\sum_i W_i}}, \quad [Eq 5]$$

Here, $q_i^{initial}$ and $q_i$ are the initial CGenFF and current charges for atom $i$. The weights $W_i$ are determined by:

$$W_i = \frac{1}{(\alpha + P_i)}, \quad [Eq 6]$$

where $P_i$ is the CGenFF penalty for atom $i$, and $\alpha = 0.1$ is a small number introduced to avoid divergence when $P_i=0$. This weighting ensures that the charge restraint term is more pronounced for atoms with low penalties, preventing substantial deviations from the initial CGenFF charges. Conversely, it has minimal impact on the optimization process for charges with higher penalties, where the initial charge assignments are considered potentially less accurate and more likely to require refinement.

Charges of symmetrical atoms had identical values during the charge optimization. Charges of aliphatic hydrogen atoms were not optimized, in accord with the standard CHARMM method with aliphatic hydrogen atoms having a charge of +0.09e. Charges on hydrogen atoms bonded to aromatic carbons were permitted to deviate from the initial CGenFF charge. However, owing to a small penalty associated with partial charges on these atoms, the optimized values exhibit only minor deviations from their initial values.

## Optimization of bonded parameters

Three variations of PES scans were used as target data to optimize three classes of bonded terms. The first category encompasses harmonic terms: bonds, valence angle, and improper dihedral angles. Non-rotatable dihedral angles, characterized by a deep minimum and elevated energy for small deformations (e.g. dihedrals in aromatic rings or amide groups), are also part of this class. The second category involves rotatable dihedrals, exhibiting a shallow energy surface with relatively small barriers between minima. Consequently, these dihedrals may undergo significant fluctuations during simulations, and thus accurate treatment of the full range of accessible conformations around these dihedrals is important for conformational sampling in molecular simulations. The third category comprises dihedrals within non-conjugated rings. Although full rotations around these dihedrals are
prohibited due to structural constraints, they may exhibit shallow PES with multiple minima. Rotations around such dihedrals induce changes in other ring dihedral angles, resulting in diverse ring conformations or puckering states. This category also includes dihedrals not confined to rings, where complete revolutions are not possible without disrupting the chemical structure. An example of such dihedrals involves nitrogens in pyramidal conformations, where the nitrogen is one of two central atoms in the dihedral, and pyramidal inversions are feasible.

For the dihedral PES, each molecule is characterized by *M* rotatable dihedrals and *N* dihedral angles describing ring conformations. These dihedrals define the molecule's overall conformation, when energy minimization is performed along all other degrees of freedom. In this study, we considered only one-dimensional PES scans, wherein scanning a dihedral angle involves keeping all other soft and in-ring dihedral angles constant at their values in the energy-minimum conformation. For the first category of terms, corresponding to rigid degrees of freedom, all soft and in-ring dihedral angles are constrained to their respective energy-minimum values.

Optimization of bonded parameters was performed with the developer program.⁷⁵ Each conformation for the MM calculation was extracted from the QM scan and minimized with a harmonic restraint force constant of 5·10⁴ kcal·mol⁻¹·Å⁻² or 5·10⁴ kcal·mol⁻¹·radian⁻² on the target bond and valence angle, respectively. QM PES energies were scaled down by 0.889, except for soft-dihedral PES, in accord with scaling applied to vibrational frequencies performed in the original CGenFF optimization.²⁶, ⁷⁶ At each iteration of bonded parameter optimization, adiabatic PES scans were performed with CHARMM⁷⁷ using the new set of MM parameters. The target function included RMSD between QM and empirical PES energies, QM and MM geometries, and restraints on the parameters themselves. Restraints on the dihedral force constants were implemented to limit their magnitudes and eliminate unnecessary multiplicities. The adjustment of MM parameters continued until the target function could no longer be minimized further.

## Optimization of harmonic parameters.

Parameters for the intramolecular terms described by harmonic potentials; bonds, valence angles, and improper dihedrals, as well as non-rotatable dihedral angles were optimized using the following protocol. The initial guess for force constants was provided by the CGenFF program. The initial equilibrium values for bonds, valence angles were taken directly from MP2/6-31G(d) geometries (MP2/6-31+G(d) for anions). During the PES scans performed by varying one stiff degree of freedom the potential energy may become very high, even for relatively small deformations. Such high-energy regions of PES are not sampled during typical MD simulations. To ensure the parametrization focuses only on relevant PES regions, we employed a method from our previous work.⁷⁸ This method limits the extent of deformations and corresponding energies by estimating force constants based on initial distortions and then correcting the maximum allowed distortions using the equation:

$$ \Delta x' = \sqrt{2 \Delta E_{max}/k}, \quad [Eq 7] $$

where *k* = 2 (*E*(*Δ*x*) − *E*₀)/Δ*x*². Δ*x* and Δ*x'* are the initial and adjusted maximum distortions, respectively; *E*₀ and *E*(*Δ*x*) are the minimum energy and energy of the deformed structure. Δ*E*max defines the highest energy on
the scanned PES. To optimize each bonded term five points were used for each PES scan equally spaced in the range of $x \in [x_0 - \Delta x', x_0 + \Delta x']$, including the minimum energy structure at $x = x_0$. 2.0 kcal·mol⁻¹ was used for $\Delta E_{max}$ in Equation 7. All PES scans were performed at the MP2/6-31G(d) model chemistry and MP2/6-31+G(d) for anions with Gaussian 16.

For the harmonic terms, MM calculations involve a two-step protocol described in our previous work, which separates the dependence of the cost function on equilibrium values from that on force constants.⁷⁵ In the first step, the conformation was optimized with constraints applied only to rotatable dihedral angles. Bond lengths, valence angles, and dihedral angles from this MM-optimized conformations were used for MM PES scans. The value $x_k^{MM}$ for scanned degree of freedom, k, was allowed to deviate from $x_k^{QM}$ and is given by:

$$x_k^{MM} = x_{0,k}^{MM} + (x_k^{QM} - x_{0,k}^{QM}), \quad [Eq 8]$$

where $x_{0,k}^{MM}$ and $x_{0,k}^{QM}$ are values in the MM and QM optimized structures without any restraints; $x_k^{QM}$ is the value in the corresponding QM structure. In addition, the cost function incorporates the deformation term introduced in our prior research,⁷⁵ which ensures accurate optimization of force constants for valence angle terms that share the same atomic center.

## Optimization of flexible dihedral parameters

To parametrize rotatable dihedral terms, adiabatic PES scans were performed for each torsion, $\chi_k$, in which the torsion angle was scanned in the range from -180° to 180° in 10° increments. During these scans, the compound was energetically optimized along all degrees of freedom, except for the soft dihedral angles. The scanned soft dihedral $\chi_k$ was constrained to the target value, while all other soft dihedrals [$\chi_{i \neq k}$] including those that define ring conformations were constrained to values corresponding to the minimum energy geometry of the model compound. This approach ensured smooth energy profiles by preventing abrupt changes caused by conformational variations along dihedral angles other than the one being scanned. Each conformation for the MM calculations was extracted from the QM scan and minimized with a harmonic restraint with the force constant of 5·10⁴ kcal·mol⁻¹·radian⁻² on the target torsion. All other rotatable dihedrals were restrained with the same force constant to the values corresponding to the QM minimum energy geometry. Using these dihedral restraints, we ensured that the QM and MM structures for each dihedral PES scan are close to each other, i.e. that we compare the same region on QM and MM PES. The dihedral parameters were optimized to achieve a minimum deviation between the QM and MM surfaces only in the low-energy regions, defined as those with energies below $E_{cut} = 10$ kcal·mol⁻¹ above the minimum energy. Conformations were excluded only if both QM and MM relative energies exceeded $E_{cut}$.

The third category of terms includes in-ring dihedrals defining the conformations of non-conjugated rings. To optimize these in-ring dihedral terms, one-dimensional adiabatic PES scans were conducted, akin to the approach used for rotatable dihedrals described above. However, in contrast, only the lower-energy region of the PES was explored. The exploration began at the dihedral angle, $\chi_{0,k}$ corresponding to the QM-optimized structure, with the energy $E_0$. Conformations corresponding to neighboring points ($\chi_{0,k} - 10°$ or $\chi_{0,k} + 10°$) on the grid
were then optimized to obtain their energies $E_i$ and $E_{i-1}$. If the energy difference ($E_i - E_{i-1}$) for a neighboring grid point $i$ was lower than a specified threshold $E_{cut}$, the next grid point (i.e. $\Delta x_{i+1} = x_i + \Delta x$ or $\Delta x_{i-1} = x_i - \Delta x$ where $\Delta x = 10^\circ$) was included for further exploration. The optimization was terminated when the conformational energy exceeded the threshold of $E_{cut} = 7 \text{ kcal} \cdot \text{mol}^{-1}$. This protocol allowed coverage of the low energy region of the PES, including the initial point ($E_0$) and all grid points with relative energies below $E_{cut}$.

It's important to note that this approach introduces the approximation that in-ring transitions can occur along a single dihedral, which is typically not the case. To mitigate this approximation, we have also incorporated puckering states and their corresponding energies, identified through a conformational search, into the parametrization process as described below.

## Inclusion of conformational minima in the optimization

During each iteration of the dihedral parameter optimization, when present, multiple conformations (called target structures) derived from conformational analysis were optimized using updated parameters. The relative energies of these target structures were incorporated into the cost function. Optimization began with the QM-optimized structures. The steepest descent (SD) algorithm was applied for 200 steps, followed by 200 steps with the Powell algorithm,⁷⁹ and finally 500 steps with the ABNR algorithm.⁸⁰ To prevent significant deviations from QM-optimized structures, harmonic restraints were imposed on all soft dihedrals with a force constant of $5 \cdot 10^4 \text{ kcal} \cdot \text{mol}^{-1} \cdot \text{radian}^{-2}$. These restraints were defined relative to the dihedral values in the target QM geometry. This approach can be viewed as a multi-dimensional dihedral PES scan, where only the PES minima are considered.

Upon completion of the MM optimization, the energy was computed, considering all terms except for the harmonic restraints. To compare the energies of target structures, we employed an exponential weighting RMS deviation formula:

$$F_{CM} = \sqrt{\frac{\sum_i w_i \cdot (E_i^{MM} - E_i^{QM})^2}{\sum_i w_i}}, \quad [\text{Eq 9}]$$

where weights $W_i = \exp(-\alpha \cdot \min(E_i^{MM}, E_i^{QM}))$ and $\alpha = 0.2$. This formula increases the contribution of low-energy conformational minima, emphasizing those most frequently sampled in simulations, thereby focusing on the low-energy PES region.

In the optimization protocol, parameters associated with dihedral angles including rotatable dihedral angles, along which large conformational differences can occur were first optimized. Following this, the bonded terms associated with other degrees of freedom were adjusted. The steps were repeated iteratively until no significant improvement was obtained.

## 4-membered ring-containing compounds

Following the parametrization involving the 1,390 model compounds described above, it was determined that the CGenFF coverage should be extended to include fused ring systems involving 4-membered rings.
Parameter optimization focused on the bonded terms. Model compound QM geometry optimizations were performed using Psi4⁸¹ at the MP2/aug-cc-pVDZ model chemistry to provide target data for optimization, including minimized geometries and PES scans as implemented in the FFParam package.⁵⁵, ⁸² In addition, vibrational calculations were performed with the MP2/6-31G(d) model chemistry following geometry optimization at the same level of theory. All MM model compound geometry optimizations, PES and vibrational calculations were performed using CHARMM.⁷⁷, ⁸³

## Retraining of the CGenFF program

Optimization of the parameters for the 1,390 compounds and the 4-membered ring test compounds yielded an extended set of parameters for training the CGenFF program. The results presented below include the outcomes of the direct optimization of those parameters associated with the training set compounds with respect to the QM target data. However, in practice CGenFF is based on the parameters generated by the CGenFF program. Accordingly, using the extended set of parameters, including bonded terms and partial atomic charges, the bond-angle-dihedral charge increments in the ML model were retrained based on the extended set of bonded parameters and partial atomic charges. This process resulted in CGenFF version 5.0. Following this, MM data for the 1,390 training set compounds were recalculated and compared with the corresponding target QM data. Furthermore, CGenFF was used to generate parameters for >900 drug-like molecules approved by the FDA employed in our previous study on the development of a Deep Neural Network to predict electrostatic parameters for the Drude polarizable FF.⁸⁴ Additionally, a collection of molecules for which experimental pure solvent properties and hydration free energies are available, and which have been previously investigated in several studies,⁸⁵, ⁸⁶ was included for validation. Expanding the set of molecules beyond those used in the training is crucial, as it also allows for analysis of the impact of the extended training set on the penalties obtained from the CGenFF program.

Training of the CGenFF program involves the optimization of charge increments related to bonds, angles and dihedrals using a restrained least-squares fitting algorithm as described previously.⁴³,⁴⁴ The in-house developed code used for training of previous CGenFF versions was extended to accommodate the dramatically increased amount of target data considering memory management. 42 independent least-squares fittings were conducted with a range of restraining biases ($X \cdot 10^{-Y}$ with X being 1, 2 and 5, and Y being an integer number between 4 to 17). From the resulting ML models, the one that yields partial atomic charges with the highest correlation with the charges generated by the previous CGenFF version was kept as the final optimal solution. After a preliminary round of training, some problematic atomic charge assignments for a few molecules were identified and found to be caused by inconsistency in resonance structure determination in some rare cases between the training stage, which is based on the input topology file, and the prediction stage, which relies on the input molecular structure and the CGenFF program. A few specific rules were introduced to the training code to fix this mismatch by resetting the formal charges. The final retrained CGenFF program, version 5.0, was tested using the ZINC 10K dataset used in the original CGenFF paper²⁴ and penalty scores were examined and compared with previous versions.
## Analysis of FDA approved drug-like molecules

To assess the performance of retrained CGenFF version 5.0 against the older version 2.5.1, QM target data for the 921 FDA-approved drug molecules were generated. These drug molecules are larger (200–700 Da) than those in the training dataset and contain multiple ring systems connected by various linkers with rotatable bonds. Starting from QM optimized structures at MP2/6-31+G(d), a maximum of 10 conformers were generated for each compound using the CREST software,⁶⁵ ensuring that each conformer is a local minimum on the PES. QM energy and molecular dipole moment of all conformers were obtained by performing single point calculations at the MP2/6-31+G(d) model chemistry. For the compounds containing iodine atoms, the Def2TZV basis-set was used. Topologies for 921 compounds were generated using both CGenFF version 5.0 and version 2.5.1. The MM dipole moment was calculated for the lowest-energy conformer of each compound without optimization and compared to the corresponding QM dipole moment. To evaluate geometric accuracy, the lowest energy conformer for each compound was optimized using parameters obtained from CGenFF 5.0 and 2.5.1 followed by RMSD calculations between non-hydrogen atoms of the QM structure and the MM-minimized structure. To compare the PES, rotatable bonds were identified using the RotatableBondSmarts library in RDKit.⁶⁰ Harmonic restraints were applied to the atoms involved in the rotatable bonds and their immediate neighbors during MM minimization of each conformer. This approach allows computation of the RMSD between the relative QM and MM energies for all conformers of each compound:

$$RMSD = \sqrt{\frac{\sum_{i=1}^{N} (E_{QM}^{i} - E_{MM}^{i})^{2}}{N}}, \quad [Eq 10]$$

where $E_{QM}^{i}$ and $E_{MM}^{i}$ are the energies of conformers in kcal·mol⁻¹, relative to the energies of the lowest QM and MM energy conformers for each compound, respectively.

## Pure solvent simulations

Condensed phase pure solvent properties were calculated using a cubic cell containing 216 copies of the molecule being studied. The cell was initially constructed by placing a copy of the respective molecule on each grid point of a cubic 6×6×6 lattice, whose grid spacing was determined by the molecule size. In the presence of periodic boundary conditions,⁵⁷ each system was minimized for 500 steepest descent steps and for 500 steps using the ABNR⁸⁷ using the CHARMM program.⁷⁷, ⁸³

MD simulations using OpenMM 8.0.0⁸⁸ were performed on the pure solvents to equilibrate to the relevant temperature for 1 ns in an NPT ensemble using a 1 fs timestep with the Langevin integrator⁸⁹ with a friction coefficient of 20 ps⁻¹. The Particle Mesh Ewald (PME) method⁹⁰ was used for the treatment of the Coulomb interactions with a real space cutoff of 12 Å and a tolerance error of 0.0005. For the LJ interactions, a force switching function⁵⁷ was applied over the range of 10–12 Å, and an isotropic long-range correction was used to account for LJ interactions beyond the cutoff distance. The lengths of covalent bonds to hydrogen atoms were constrained to their equilibrium values. An isotropic Monte Carlo barostat⁹¹, ⁹² to maintain constant pressure at 1 atm with a coupling frequency of 20 steps was also used. This equilibration was followed by a production
simulation of 10 ns with the same parameters except the timestep was increased to 2 fs. The molecular volume, $\langle V_{liq} \rangle$ was calculated from the average volume of the 10 ns simulation volume of the simulation box during the 10 ns production simulation divided by the number of monomers in the system, $N_{mol}$. Heats of vaporization, $\Delta H_{vap}$, were calculated using:

$$ \Delta H_{vap} = -\frac{\langle U_{liq} \rangle}{N_{mol}} + \langle U_{gas} \rangle + RT, \quad [Eq 11] $$

where T is the temperature, R is the gas constant, $\langle U_{liq} \rangle$ is the average potential energy of the liquid obtained from the 10 ns production MD simulations, and $U_{gas}$ is the gas phase energy. The gas-phase energy was computed from separate 200 ps vacuum MD simulations of individual monomers for all 216 molecules in the cubic box, with the final energy calculated as the average of the 216 individual potential energies averaged over the 200 ps simulation. These gas-phase simulations were performed using Langevin dynamics with a 2 fs integration time step, a friction coefficient of 20 ps⁻¹, and no truncation of nonbonded interactions.

## Free energies of hydration

Hydration free energies (HFE) were computed for 316 neutral molecules with known experimental data using a replica exchange free energy perturbation (FEP) method implemented in OpenMM.⁹³,⁹⁴ For each molecule, the HFE was computed using the CGenFF 2.5.1 and 5.0 parameters for the molecules and the CHARMM TIP3P model for water.⁷³ Calculation of HFE using FEP involved introduction of a series of intermediate states between two thermodynamic end states. To describe these intermediate states, a single non-physical parameter $\lambda$ going from 0 to 1 is used to describe the Hamiltonian $H(\lambda)$, with $\lambda = 0$ and $\lambda = 1$ representing the two end state Hamiltonians. The HFE calculations requires free energy calculations to be performed in the gas phase and the bulk aqueous phase. The initial FEP setup script in FFParam⁵⁵ utilizes Packmol to generate a simulation box with 580 TIP3P water model solvating the concerned molecule. The simulation box was equilibrated for 1 ns followed by 5 ns of production run. The replica-exchange alchemical free energy simulations employed for the HFE calculations followed a two-stage alchemical protocol. In the first stage, charges were linearly scaled to annihilation using coupling parameters ranging from 0.0 to 1.0 in steps of 0.04, while LJ interactions were annihilated in the second stage using a softcore potential.⁹³,⁹⁴ Simulations were run for 1 ns/replica using OpenMM 8⁸⁸ and openmmtools 0.23.1.⁹⁵ Bulk phase simulations were carried out at 298 K at 1 atm using Langevin dynamics with a friction constant of 5 ps⁻¹ applied to the non-hydrogen atoms of all molecules at 298 K, 1 atm pressure. The PME method with a real space cutoff of 12 Å was employed⁹⁶ to treat long-range electrostatics. LJ interactions were truncated at 12 Å with a VFSwitch smoothing function with long-range contributions included through an isotropic correction. A 2 fs integration time step was used, and all covalent bonds involving hydrogen atoms were constrained using SHAKE.⁹⁷ The bulk phase and gas phase free energies were calculated using the multistate Bennett acceptance ratio (MBAR)⁹⁸ with absolute HFEs obtained as the difference between the bulk and gas phase free energies.⁹⁴
# RESULTS

The present study involved a substantial expansion of the training set of molecules and parameters in CGenFF and for use in the CGenFF program. In addition, CGenFF was expanded to include fused rings that include a 4-membered ring along with optimization of the bonded parameters, which were also included with the new version of CGenFF. The expansion of the training set primarily included a collection of diverse low molecular weight compounds selected from the ZINC20 database that required parameters not already present in CGenFF. These molecules, including fragments obtained from selected compounds with multiple ring systems, were then subjected to a systematic optimization of the bonded parameters and the partial atomic charges. Initial guesses of the parameters were from version 2.5.1 of the CGenFF program. All bonded parameters with a penalty > 10 not previously in CGenFF were optimized, as were all partial atomic charges. For the 4-membered ring systems, all bonded parameters previously not present in CGenFF were optimized. Once the parameter optimizations were complete, the CGenFF program linear interpolation scheme for calculating partial atomic charges was retrained using the extended set of optimized parameters. This effort yielded version 5.0 of the CGenFF program. Accordingly, results below include 1) those based on parameters from CGenFF version 2.5.1, 2) those using parameters fully optimized in the present study, and 3) those generated by the retrained version 5.0 of the CGenFF program. In addition, an external set of FDA-approved compounds and compounds for which experimental pure solvent properties and hydration free energies are available were used for validation of parameters generated using the CGenFF program.

## Model compounds used for parameter optimization

This study involved the parametrization of 1,390 small molecules as detailed in Table 1 along with additional compounds that included 4-membered rings described below. The chemical structures and supplementary results for each molecule are provided in the MacKerell lab GitHub repository. The ensemble of 1,390 compounds includes 773 molecules sourced directly from the ZINC20 database, 139 neutralized ZINC20 molecules, and 478 fragments extracted from parent compounds in the ZINC20-direct set into the dataset.

<table><caption>**Table 1.** Molecules parametrized in this work and their selected properties.</caption><thead><tr><th>Compound set</th><th>Number</th></tr></thead><tbody><tr><td>ZINC20-direct</td><td>773</td></tr><tr><td>Neutralized</td><td>139</td></tr><tr><td>Fragments</td><td>478</td></tr><tr><td>Total</td><td>1,390</td></tr></tbody></table><table><thead><tr><th colspan="2">Properties</th></tr></thead><tbody><tr><td>Rotatable bonds</td><td>&lt;2</td></tr><tr><td>Molecular Weight</td><td>&lt; 250 Daltons</td></tr><tr><td>H-bond acceptors</td><td>&le; 6</td></tr><tr><td>H-bond donors</td><td>&le; 5</td></tr><tr><td>LogP</td><td>-3.94 to 7.50</td></tr><tr><td>Total charge</td><td>-3 to 2</td></tr></tbody></table>
Partial atomic charge optimization

*Case study: Optimization of atomic charges for 2300R6 and 23003N.* In this section, charge optimization is exemplified using two neutral molecules with CGenFF residue names: 2300R6 (ZINC000238583940) and 23003N (ZINC000001568873), whose initial CGenFF models showed significant deviations from the target QM data. 2300R6 corresponds to 3-fluoro-2-mercapto-5-(trifluoromethyl)pyridine, and is a calcium channel blocker used in pain treatment.⁹⁹ 23003N corresponds to (5R)-1,3,8-trimethyl-1,3,6,8-tetrazaspiro[4.4]nonane-2,4,7,9-tetrone, and it contains a spirocyclic scaffold of increasing interest for medicinal chemistry¹⁰⁰ that had not been explicitly optimized for CGenFF prior to this study.

For each compound, the QM-minimized geometry was used to compute the dipole moment and compound-monohydrate water minimum interaction energies and distances with QM, using one to four water orientations per probed atom, as shown in Figure 2. This minimized geometry was also used to compute the MM water interactions, dipole moment and the RMSD between QM and MM ESPs. Charge optimization started from initial CGenFF 2.5.1 partial charge predictions, which were iteratively refined to maximize the reproduction of the QM target data. Initial and optimized charges, as well as MM and QM data for 2300R6 and 23003N are given in Tables S1 and S2. For 2300R6, six atoms had charge penalties > 10 with N10 having a penalty of 24.4. The water interaction energy for N10 deviated by 8.14 kcal·mol⁻¹ from the QM value. Optimization of the charges for the entire molecule resulted in the charge of N10 changing from 0.369e to -0.132e, and the water interaction energy difference improved from 8.14 kcal·mol⁻¹ to 0.93 kcal·mol⁻¹. However, the primary contribution to this improvement was not the change in the charge of N10, but rather the adjustment in the charge of S12, which decreased from -1.702e to -0.535e. Notably, the overall RMSD of all water interaction energies decreased from 7.80 to 0.60 kcal·mol⁻¹. In addition, substantial improvements were observed in both the magnitude and orientation of the dipole moment, as illustrated in Figure 2A. For 23003N (Table S2) ten atoms had penalties > 10 with values > 50 occurring for N5, C6 and N15. Although the penalty for H20 linked to N5 was low, 3.7, the difference in the interaction energy of the N5-H20 moiety relative to the QM data was significant 9.25 kcal·mol⁻¹, which improved to 0.04 kcal·mol⁻¹ after charge fitting. This involved a small change in the charge in H20 from 0.356e to 0.358e, while N5 underwent a larger charge adjustment from -0.281e to -0.564e. The overall agreement between all the QM and MM water minimum interaction energies improved, with the RMSD decreasing from 3.74 to 0.32 kcal·mol⁻¹. Again, the improvement in the dipole moment magnitude and orientation was substantial as shown in Figure 2B. The QM dipole moment and MM dipole moments with the optimal charges for 2300R6 were 5.8 and 5.8 D, respectively, while those values for 23003N were 1.2 and 1.2 D. Furthermore, the RMSD between QM and MM ESPs decreased significantly from 29.0 to 2.8 kcal·mol⁻¹·e⁻¹ for 2300R6, and from 19.0 to 1.5 kcal·mol⁻¹·e⁻¹ for 23003N.

An important point is that the CGenFF penalty is not necessarily an indication of the extent of charge modification required for each atom. For instance, for atoms S12 and N10 of 2300R6 with the penalties 5.4 and
24.4, respectively, the absolute differences between the initial and optimized charges were both large; 1.17e and 0.50e, respectively (Table S1). With compound 23003N the absolute change in the charge of C6 with a penalty of 96.6 was 2.04e while that of N15 with a penalty of 83.8 was only 0.17e (Table S2). These examples highlight that while large changes in charges are observed for some high-penalty atoms, the magnitude of the CGenFF penalty, which is based on analogy with known atom types and connectivities, does not necessarily reflect the accuracy of the parameter itself. The overall improvement in agreement with QM target data depends not on the magnitude of individual charge modifications per atom but rather on the ensemble of charge modifications across the entire molecule.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Interaction orientations of water with parametrized molecules (A) 3-Fluoro-2-mercapto-5-(trifluoromethyl)pyridine (CGenFF: 2300R6) and (B) (5R)-1,3,8-trimethyl-1,3,6,8-tetrazaspiro[4.4]nonane-2,4,7,9-tetrone (CGenFF: 23003N). Insets show chemical structures of the compounds. Model compound gas phase MP2/6-31G(d) minimized geometries were used for charge optimization. Note that only a single water molecule is interacting with the compound during each calculation; all water molecules are shown simultaneously for convenience. Water interactions with aliphatic hydrogens are not shown for clarity. Dipole moments computed</figcaption></figure>
with QM MP2/6-31G(d), MM initial and optimized charges are shown by arrows colored in cyan, gray and orange, respectively. In panel (B) the MM dipole with optimal charges aligns closely with QM dipole moment.

*Global results for the partial atomic charge optimization.* The partial atomic charges were optimized for all the compounds excluding those compounds with maximum charge penalties < 10, yielding a total of 1,350 compounds. Optimization targeted the QM water-compound interactions, the dipole moment magnitude and its orientation, and ESPs with varying weights as described above. Results are present for the direct optimization of the charges. Subsequently, once the optimized charges were available for the 1,350 compounds, along with all the optimized bonded parameters, the partial atomic charge prediction algorithm in CGenFF was retrained using the charges and bonded parameters from all 1,390 compounds along with the compounds previously included in the training set. The results obtained using charges generated with retrained CGenFF are also provided. As described below, during the initial retraining outliers were identified that were associated with atom type assignments that were subsequently corrected. The retrained results presented below are based on these corrected assignments.

<table><caption>**Table 2.** Statistics for intermolecular parameter development and agreement with respect to selected target data for all model compounds.</caption><thead><tr><th rowspan="2">Property</th><th rowspan="2">N points</th><th colspan="2">RMSD</th></tr><tr><th>initial/optimized/retrained</th><th>initial/optimized/retrained</th></tr></thead><tbody><tr><td><sup>a</sup>norm of μ</td><td>932</td><td>2.44/0.68/1.14</td><td>1.65/0.43/0.81</td></tr><tr><td><sup>b</sup>direction of μ</td><td>932</td><td>47.9/4.8/20.1</td><td>32.2/1.7/10.8</td></tr><tr><td><sup>c</sup>water-compound E<sub>int</sub></td><td>34,313</td><td>1.67/0.56/0.68</td><td>1.34/0.50/0.60</td></tr><tr><td><sup>d</sup>water-compound d<sub>min</sub></td><td>34,313</td><td>0.33/0.18/0.19</td><td>0.24/0.15/0.16</td></tr><tr><td><sup>e</sup>φ<sub>elec</sub></td><td>1,350</td><td>7.20/2.31/3.11</td><td>5.97/2.14/2.76</td></tr></tbody></table>

<sup>a</sup>The magnitude of the dipole moment (μ) is given in Debye; <sup>b</sup>angle (°) between the *ab initio* and empirical dipole moment vectors; <sup>c,d</sup>probe water-model compound minimum interaction energies, *E*<sub>int</sub>, and distances, *d*<sub>min</sub>, are in kcal·mol<sup>-1</sup> and Å, respectively; <sup>e</sup>electrostatic potential is in kcal·mol<sup>-1</sup>·e<sup>-1</sup>.

The statistics for water-compound interactions over all the compounds are provided in Table 2. A total of 34,313 monohydrate probe water-model compound interaction complexes were used as target data. This includes 9,659 probe water-model compound interactions for charged and 24,654 for neutral compounds. The RMSD for interaction energies with the initial, optimized and retrained charges was 1.67, 0.56, and 0.68 kcal·mol<sup>-1</sup>, respectively, while the mean absolute error (MAE) was 1.34, 0.50, and 0.60 kcal·mol<sup>-1</sup>, respectively. In this work, several probe water orientations were considered for a compound atom capable of forming hydrogen bonds (H-bonds).<sup>82</sup> The advantage of probing using different orientations of the probe water is that it provides additional target data on the interaction between the water dipole moment and the dipole moment of the molecule, as well as interactions with neighboring groups within the compound. Certain water orientations exhibit significantly more favorable interaction energies due to interactions with other groups within the molecule, making them more challenging to replicate using the simple additive form of the force field (e.g. N10 atom in Figure 2A and Table
S1). Additionally, the dataset includes numerous charged molecules, with 9,659 probed water interactions, all characterized by highly favorable water-interaction energies. For these charged species, the QM-MM differences in interaction energies computed in vacuum are larger in absolute magnitude compared to neutral compounds, yet the optimized force field still provides accurate results. This contributes to the RMSD between QM and MM interaction energies of 0.56 kcal·mol⁻¹ for optimized and 0.68 kcal·mol⁻¹ for retrained charges obtained in this work being slightly higher relative to the 0.34 kcal·mol⁻¹ reported in the original 2010 CGenFF paper.²⁶ Figure 3A compares QM and MM interaction energies for all probe waters. The slope of the interaction energies computed with the three types of charges is close to one, indicating that the force field model can effectively reproduce solvent interactions across a broad spectrum of molecules. However, the R² values are 0.87, 0.98, and 0.97 for initial, optimized, and retrained charges, underscoring the significance of optimization and retraining on a larger dataset. The frequency of smaller energy absolute differences and RMSD also increases when using optimized and retrained charges, as depicted in Figure 3B and 3C.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** QM and CHARMM water minimum interaction energies for the compound-water monohydrates computed with initial, optimized, and retrained CGenFF charge sets. QM minimum interaction energies for the neutral compounds were scaled by 1.16 as described in the text. A) *Ab initio* vs CHARMM energies for all 34,313 probed water-interactions for both charged and neutral molecules. The linear regression line between QM and</figcaption></figure>
CHARMM data is shown by the solid line. The diagonal dashed lines represent deviations of ±2 kcal·mol⁻¹ from the regression line. B) Probability distribution of absolute differences between *ab initio* and CHARMM energies for all water interactions. C) Probability distribution of RMSD between *ab initio* and CHARMM energies computed for each molecule.

The RMSD for minimum-energy interaction distances was 0.33 Å with the initial charges, which decreased to 0.18 Å with the optimized atomic charges and 0.19 Å with the retrained charges. The consistency in interaction distances aligns well with earlier findings for 2010 CGenFF, showing an RMSD of 0.20 Å for distances.²⁶ Similar to the water interaction energies, the frequency of smaller absolute differences and RMSD for interaction distances increased when employing optimized and retrained charges, as demonstrated in Figure 4A and 4B. The larger differences with respect to the QM data in the present study versus the 2010 CGenFF study can be attributed to the inclusion of a large number of interactions with water molecules with non-polar atoms, such as with methyl groups. This discrepancy is likely due to the empirical minimum energy offset distance of 0.2 Å for the HF/6-31G(d) model chemistry being associated with interactions involving polar neutral moieties.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of minimum interaction energy distances for compound-water monohydrates obtained with initial, optimized and retrained CGenFF atomic charges. (A) Probability distribution of QM and CHARMM absolute distance differences for all water interactions used for charge optimization and (B) probability distribution of distance RMSD per molecule. QM minimum interaction distances for the neutral compounds were offset by -0.2 Å for non-hydrogen atoms of neutral molecules as described in the text.</figcaption></figure>

The dipole moment optimization was included only for neutral compounds and the QM dipole moments were scaled by ~30% during optimization to overestimate the gas phase dipole moments. This was performed to account for the implicit inclusion of the contribution of electronic polarization of molecules in an aqueous environment consistent with the standard CHARMM protocol.²⁶ The dipole moments with the optimal and retrained charges show significant improvement compared to those calculated using the initial CGenFF charges. The RMSD between scaled QM and MM dipole moments for neutral compounds is 2.44, 0.68, and 1.14 Debye for the initial, optimal, and retrained charges, respectively. The agreement of the orientation of the dipole moments is also enhanced; the average angle between the QM and MM dipole moments across the compounds is 47.9°, 4.8°, and 20.1° for the initial, optimal, and retrained sets of charges, respectively. The RMSD for both the magnitude and direction of the dipole moment in this study is comparable or better than the agreement of 1 Debye and 8.5°
achieved with the 2010 CGenFF force field.²⁶ The initial charges yield a QM vs MM dipole moment magnitudes slope of 0.68, whereas optimized and retrained charges produce a slope close to one, as illustrated in Figure 5A. This suggests that the optimized and retrained force field model can accurately replicate the QM dipole moment for a variety of molecules. Additionally, the R² improves from 0.47 with initial charges to 0.98 and 0.90 with optimized and retrained charges, respectively. Figures 5B and C also show an increase in the frequency of dipole moment angles and their magnitudes exhibiting smaller errors when using optimized and retrained charges.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison between the scaled QM and CHARMM dipole moments for 932 neutral molecules and ESP for all 1350 molecules. A) Scaled QM (MP2/6-31G(d)) vs CGenFF dipole moment magnitude. The linear regression is shown by dotted lines. B) Probability distribution of absolute differences between scaled QM and CGenFF dipole moment magnitude, and C) angle between the two computed with initial, optimized and retrained charges for the neutral fragments. D) Probability distribution of RMSD between QM (MP2/6-31G(d)) and CHARMM electrostatic potentials per molecule obtained with initial, optimized and retrained charges.</figcaption></figure>

ESPs were incorporated as additional target data during the charge optimization as previously described.⁶⁴

⁷⁰ This is particularly important for charged molecules where fewer atoms are available for interactions with water, and the dipole moment is not included in the target data. The RMSD between MM and QM ESPs averaged over all 1,350 compounds was 7.20, 2.31, and 3.11 kcal·mol⁻¹·e⁻¹ with the initial, optimized, and retrained CGenFF charge sets, respectively. As shown in Figure 5D, the ESP RMSD shows significant improvement compared to
the initial values. Specifically, 80% of the molecules exhibit an ESP RMSD from QM of less than 2 kcal·mol⁻¹·e⁻¹, while the ESP RMSD with the initial charges ranges from 1 to 11 kcal·mol⁻¹·e⁻¹.

The optimized charges along with the optimized bonded parameters were then used for pure solvent simulations as part of the validation which showed selected compounds to have substantial disagreement with the experimental heat of vaporization (Figure S1). Subsequent analysis of the QM and MM water-compound minimum interaction energies identified some outliers associated with compounds for which problematic pure solvent properties were observed (Figure S2). This was associated with certain carbon atoms having partial charges exceeding an absolute value of one (Figure S1). These large-magnitude charges were associated with “buried” carbons not directly on the surface of molecules for which the CGenFF penalties were large such that the restraints during the automated optimization were small, which allowed the large-magnitude charges to emerge from the optimization. Though these charges may not pose a problem for the compounds explicitly parametrized, they hinder transferability to other molecules with the retrained CGenFF program. Consequently, a second round of optimization was undertaken for 102 problematic molecules involving human intervention to ensure accurate charge distributions. Comparing the results from both rounds of optimization showed that outliers from the first fully automated optimization were largely eliminated (Figure S2). The elimination of the anomalous partial atomic charges associated with various carbon atom types for the validation molecules is evident when the distribution of the charges on individual atoms (Figure S3).

The relationship of the penalty term to the quality of the partial atomic charge can be considered for the full set of compounds (Figure S4). As discussed above, penalties reflect the analogy of the new molecule's topology to those parametrized in the CGenFF force field, where larger values signify low similarity. As shown in Figure S4, the correlations are small. When the differences in all the charges in the molecule based on the RMSD are compared with the maximum penalty per molecule (Figure S4A), the correlation is negligible with R² = 0.07. However, some correlation exists between the penalty values and the magnitude of the changes in the charges for the individual atoms (Figure S4B) with R² = 0.22. Accordingly, there is some relationship between the magnitude of the charge penalty and the accuracy of individual charges. Thus, the charge penalties may indicate the need to undertake charge optimization although in many cases those charges may be of sufficient accuracy. Thus, in high-throughput scenarios where charges need to be rapidly generated for large numbers of molecules, the charges directly from CGenFF may generally be considered acceptable with the new retrained version v5.0 offering additional improvements over the earlier version of CGenFF v2.5.1. However, when specific molecules are being investigated, checking the charges and performing optimization when large penalties are encountered is suggested.

## Optimization of the bonded terms

The parametrization of the bonded terms was performed targeting the QM potential energy surfaces. Optimized parameters included the equilibrium values and the force constants. In the case of the dihedrals the phase and multiplicity of non-rotatable dihedral angles were taken from the CGenFF parameter assigned by analogy. For dihedral angles in conjugated systems the multiplicity was set to two with a phase angle of 180°. For
improper terms, the equilibrium values were set to zero. Bonded parameters exceeding a CGenFF penalty of 10 were subjected to optimization. In selected instances, parameters with lower penalties underwent optimization to enhance agreement with the target QM data. Table 3 gives an overview of the number of bonded parameters considered in this work, including those explicitly optimized due to high penalties and those parameters directly adapted without further modification.

<table><caption>**Table 3.** Bonded terms introduced into CGenFF in this work</caption><thead><tr><td>Term</td><td>Optimized</td><td>Adapted</td><td>Total</td></tr></thead><tbody><tr><td>Bond</td><td>422</td><td>156</td><td>578</td></tr><tr><td>Angle</td><td>2,090</td><td>1,852</td><td>3,942</td></tr><tr><td>Dihedral</td><td>8,329</td><td>4,491</td><td>12,820</td></tr><tr><td>Improper</td><td>139</td><td>138</td><td>277</td></tr></tbody></table>

*Case study: Optimization of bonded term parameters for molecule 2300P8.* Diazepines are well-known to be of importance for medicinal chemistry¹⁰¹ and in this section the bonded terms optimization is exemplified for 2300P8 (ZINC000143190822), a diazepine that corresponds to (5aS,9aR)-3-nitro-5a,6,7,8,9,9a-hexahydro-1H-1,5-benzodiazepine. 2300P8 has an interesting scaffold containing rotatable dihedrals both inside and outside the ring and having low analogy with CGenFF, with 60 of its bonded terms having penalties up to a maximum of 181. It also has multiple conformational minima for different puckering states of the ring, for which a QM-MM agreement was obtained after parameter optimization as observed in Figure 6. MM and QM PES scans were performed for 25 bonded terms for which the penalty was higher than 10. These one-dimensional PES scans were conducted near the global minimum conformation, with restraints applied to the five rotatable dihedrals to prevent conformational changes along these soft degrees of freedom. Example PES scans for each type of bonded term are included in Figure 6, where the agreement between both equilibrium values and energies are obtained between QM and MM after parameter optimization. Equilibrium values and PES scan energy RMSD for all 25 optimized bonded terms are given in Table S3 showing the overall improvement in the treatment of the compound. While additional restraints on selected dihedrals were used during the parameter optimization, the final MM conformations of local minima obtained in the absence of restraints were in good agreement with the QM geometries as shown for 2300P8 in Figure 7.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Potential energy scans (PES) for selected degrees of freedom in molecule 2300P8 ((5aS,9aR)-3-nitro-5a,6,7,8,9,9a-hexahydro-1H-1,5-benzodiazepine). Example PES scans are presented for a bond (A), for a valence angle (B), for an improper torsion (C), for a soft dihedral angle associated with the rotation of the nitro group (D), for an in-ring soft dihedral (E) and for a stiff dihedral (F). The degree of freedom along which the scan was performed is highlighted in yellow. Rotatable dihedrals were fixed unless the PES scan was along that specific torsion.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Minimum energy conformations of 2300P8 ((5aS,9aR)-3-nitro-5a,6,7,8,9,9a-hexahydro-1H-1,5-benzodiazepine) obtained with MP2/6-31G(d) and MM. For MM calculations, restraints on soft dihedrals were applied during optimization of the bonded terms but were excluded in the final energy minimization, resulting in the conformations shown. The 3D representations in the green rectangle correspond to the lowest-energy conformation. Conformations obtained with both the initial and optimal force field parameters are included.</figcaption></figure>

**Global results for bonded parameter optimization.**

The bonded parameter optimization procedures exemplified in the preceding section were applied to all 1,390 compounds using the automated approach. Table 4 compares overall agreement between the QM and MM data for equilibrium structures. In instances where multiple conformations were revealed during conformational analysis, only the lowest energy conformation was included in this analysis. As may be seen significant improvements were obtained for all the analyzed metrics. The minimum-energy structures of the compounds showed systematic improvements as indicated by the improved RMSD values, though the initial CGenFF parameters generally gave satisfactory structures. For the individual internal terms, the agreement with the target data is generally quite good. In CGenFF when optimizing specific molecules, the target level of agreement is 0.02 Å and 2° for the bonds and valence angles, respectively. As is evident from Table 4, this was generally achieved. The level of improvement seen in Table 4 was largely expected as all the parameters included were not previously optimized.

<table><caption>**Table 4.** Comparison between MM and QM optimized geometries for equilibrium structures.</caption><thead><tr><th rowspan="2">Property</th><th rowspan="2">N</th><th>MAE</th><th>RMSD</th></tr><tr><th>initial/optimal/retrained</th><th>initial/optimal/retrained</th></tr></thead><tbody><tr><td><sup>a</sup>RMSD (Å)</td><td>1,390</td><td>-</td><td>0.12/0.09/0.11</td></tr><tr><td><sup>b</sup>bond (Å)</td><td>433</td><td>0.031/0.004/0.005</td><td>0.045/0.009/0.010</td></tr><tr><td><sup>b</sup>angle (°)</td><td>2,111</td><td>3.6/0.9/0.9</td><td>5.0/2.1/2.1</td></tr></tbody></table>
<table><tr><td><sup>b</sup>dihedral (°)</td><td>7,423</td><td>9.4/4.8/4.9</td><td>16.6/8.7/8.8</td></tr><tr><td><sup>b</sup>improper torsion (°)</td><td>139</td><td>2.3/1.7/1.7</td><td>4.5/3.8/3.7</td></tr></table>

<sup>a</sup>RMS deviation between QM and MM optimized equilibrium structures for all atoms; <sup>b</sup>statistics are provided only for the degrees of freedom associated with the terms optimized in this study.

To quantify the overall improvements in the optimized parameters to model local distortions of the molecules that can occur during molecular simulations the energies of the targeted PES were analyzed. Table 5 summarizes the agreement between QM and MM energies for PES scans of targeted terms, comparing results obtained with the initial and optimized parameters. As with the geometries, improvements were obtained in all cases. The improvements in the bond and angle terms were significant but smaller than that with the dihedrals. This is consistent with the use of vibrational spectra during the previous optimization of CGenFF yielding good quality force constants requiring only minimal improvement in the present optimizations. The quality of the initial parameters also indicates the reliability of the CGenFF parameter assignment approach based on analogy, which is especially effective for the bond and angle terms. With the dihedrals the improvements are substantial with the improved correlations being notable. This outcome was expected, as the range of CGenFF coverage for dihedrals is limited by the sheer number of possible dihedral parameters associated with the chemical space covered by the force field. Similar improvements occurred with the improper as with the dihedral terms. Again, this is associated with the limited number of explicitly optimized improper parameters within the coverage of CGenFF prior to this study.

<table><caption>**Table 5.** Comparison between empirical and QM energies of PES scans in kcal·mol⁻¹</caption><thead><tr><th rowspan="2">Term</th><th><sup>a</sup>PES data</th><th><sup>b</sup>RMSD</th><th><sup>c</sup>MAE</th><th><sup>d</sup>R<sup>2</sup></th></tr><tr><th>points</th><th>initial/optimal/retrained</th><th>initial/optimal/retrained</th><th>initial/optimal/retrained</th></tr></thead><tbody><tr><td>bonds</td><td>2,110</td><td>0.60/0.22/0.21</td><td>0.27/0.13/0.13</td><td>0.79/0.97/0.97</td></tr><tr><td>angles</td><td>10,445</td><td>0.85/0.25/0.25</td><td>0.27/0.13/0.13</td><td>0.71/0.96/0.96</td></tr><tr><td>stiff dihedrals</td><td>29,080</td><td>1.58/0.30/0.30</td><td>0.49/0.16/0.16</td><td>0.53/0.95/0.95</td></tr><tr><td>rotatable dihedrals</td><td>21,349</td><td>3.42/1.48/1.50</td><td>2.15/0.86/0.88</td><td>0.69/0.93/0.93</td></tr><tr><td>in-ring dihedrals</td><td>6,734</td><td>4.04/1.75/1.74</td><td>2.25/0.78/0.79</td><td>0.79/0.94/0.94</td></tr><tr><td>improper angles</td><td>695</td><td>1.27/0.28/0.28</td><td>0.50/0.15/0.15</td><td>0.61/0.95/0.95</td></tr></tbody></table>

<sup>a</sup>Number of data points in all the PES scans performed; <sup>b,c,d</sup>RMS deviation between QM and MM energies, Mean Absolute Error, and linear regression, R², respectively.

Conformational analysis was performed for each molecule within the dataset to determine the ability of the optimized CGenFF parameters to reproduce conformational minima. Of the 1,390 molecules parametrized, 1,058 exhibit multiple conformational minima, while 332 are characterized by a single minimum. Over all molecules with multiple minima there are a total of 7,456 distinct conformational minima, averaging 7 stable
conformations per molecule. These diverse conformations were explicitly integrated into the force field parametrization of dihedral and improper terms to enhance the accuracy of CGenFF in sampling of stable conformations. The results of this parametrization are summarized in Table 6. Note that the CGenFF energies associated with these conformations, beyond being influenced by the optimized dihedral parameters, are also impacted by the current CGenFF parameters that were not optimized in this study as well as long-range nonbond energetic contributions. Overall, the expected improvements in the relative energies of the conformational states were obtained. This included substantial improvements in the correlation of the relative energies indicating the improvements in the rank ordering of the minima.

<table><caption>**Table 6.** Comparison between MM and QM relative energies of multiple conformational minima</caption><thead><tr><th></th><th>RMSD</th><th>R<sup>2</sup></th></tr><tr><th></th><th>initial/optimal/retrained</th><th>initial/optimal/retrained</th></tr></thead><tbody><tr><td>All<sup>a</sup></td><td>3.88/2.72/2.72</td><td>0.70/0.77/0.78</td></tr><tr><td>Per molecule<sup>b</sup></td><td>2.44/1.55/1.55</td><td>-</td></tr></tbody></table>

## Addition of 4-membered rings

Four-membered rings that are bridged, contain a double bond, or include non-carbon atoms with functional groups beyond carbonyl were not previously considered in CGenFF. To address this, the atom typing scheme in the CGenFF rules file was extended by developing new rules to assign the proper atom types. Subsequently, the bonded parameters for the 4-membered ring training set molecules were optimized, focusing on those parameters not already present in CGenFF. The parametrization included the training set molecules shown in Figure 8, with the initial parameters optimized based on benzocyclobutene (CGenFF: 4mem_1) and then extended to the other compounds. An additional 31 molecules were used as the test set to check both the addition to the CGenFF rules file and the optimized parameters (Figure S5).

The CGenFF program had previously been extended to include cyclobutene and related species, leading to the creation of the CG3C41 atom type, which was applied to many of the model compounds in this study. To include new 4-membered ring-containing compounds, the rules file was extended. The rule “typ CG2RC0 : ring 4 ne ( or ( bo 1 ) ( bo 2 ) ( or ( arom 6 ) ( bo 2 ring 6 ) ) )” was introduced, which first perceives the parent sp2 carbon in the 4-membered ring and determines whether it includes single or double bonds, or bridges to an aromatic or sp2 6-membered ring. This rule assigns the CG2RC0 type to the carbon atoms common to the 4-membered ring. Testing on both the training and validation set of compounds confirmed that the new rule correctly identified the molecular connectivities and assigned atom types as expected.

<footnote><sup>a</sup>All conformational minima; <sup>b</sup>Average RMSD between QM and MM minima, calculated across all molecules.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Training set compounds for the 4-membered rings.</figcaption></figure>

Selected bonded parameters previously not present in CGenFF were optimized for the model compounds in the training set (Figure 8). Equilibrium values of the bond, valence angle, and dihedral terms were optimized to match QM geometries obtained at the MP2/6-31G(d) model chemistry. Force constants were optimized targeting QM MP2/6-31G(d) PES for compound 4mem_1 and vibrational frequencies for compounds 4mem_2 and 4mem_4. Comparison of the intramolecular geometries of the 4-membered rings for the four training compounds in Tables S4 to S7 shows general consistency. Some large differences in the valence angles (> 2°) are associated with the limited number of parameters available for optimization along with the strained nature of the 4-membered rings associated with their fusion to various 6-membered rings and for a 4-membered/3-membered fused ring system (Figure 8). This is further evident from the dihedrals in compounds 4mem_2 and 4mem_4, where the non-planar QM-optimized structures (Tables S5 and S7) make it challenging to reproduce the geometries with the MM model. Optimization of the force constants yielded overall good agreement between the QM and MM PES for 4mem_1 (Figure S6). The overall shapes of the QM and MM surfaces are generally comparable, with similar locations of the minima. The largest differences occur for the C2-C1-C3 valence angle, as the parameters for this term were not optimized, being already available in CGenFF, which also contributed to larger discrepancies in other PES analyses. Importantly, the PES associated with the dihedrals are in good agreement with the QM data indicating that out-of-plane distortions of the ring are being satisfactorily modelled. Further optimization of the force constants involved the vibrational frequencies of 4mem_2 and 4mem_4, as shown in Figures S7 and S8, respectively. The MM model adequately reproduces all the low mode frequencies that exclude stretching modes involving hydrogens. With compound 4mem_4, some degradation of the frequencies occurred, associated with
optimization of parameters in other model compounds. For the valence angle parameter along the linear CCC angle in the fused 4-membered/4-membered ring in compound 4mem_4, the force constant was set to zero.

The new CGenFF rules and optimized parameters were then tested on 31 bridged 4-membered ring validation molecules (Figure S5), which included various substituents, nitrogen atoms within the ring, and diverse fused systems. In all cases the extended CGenFF rules were able to properly assign atom types to each system based on manual inspection. The MM minimum-energy geometries of the four training and 31 validation compounds were subsequently compared to their MP2/6-31G(d) QM-optimized structures, with the def2-svp basis set employed for compound 4mem_11. Following alignment of the QM and MM optimized structures the RMSD for all atoms was less than 0.1 Å for all but two molecules, which had the RMSD values > 0.6 Å (Table S7). These results indicate that the new CGenFF rules and bonded parameters appropriately treat a wide range of molecules containing 4-membered rings.

## Validation of Extended CGenFF program

### QM data for FDA compounds

To independently evaluate CGenFF 5.0 against version 2.5.1, for 921 FDA-approved drugs QM target data was generated on up to ten minimum-energy (stable) conformations identified through conformational analysis for each molecule. The data included energies, dipole moments and minimum water interaction energies and distances as described in the methods section. Table 7 summarizes the agreement between QM and both versions of CGenFF. As indicated by the low charge and bonded terms penalties for CGenFF 2.5.1, with average values of 88 and 47, respectively, the FDA molecules already have a high analogy with CGenFF 2.5.1. The presence of new parameters marginally improves the penalties, which decrease to 39 and 18, respectively. Consequently, both versions are expected to perform at a similar level. Nonetheless, the retrained model shows improvement in the water minimum interaction energies and the dipole moments, indicating potential improvements in the interactions of drug-like molecules with condensed phase environments. Probability distribution of the penalties and differences between the QM and MM results for various metrics along with dipole moment correlation analysis in Figure S9 also show improvement with some exceptions for the smaller penalties. This analysis indicates that the earlier version of CGenFF was already capable of well reproducing the QM data for the FDA compounds. This is due to those compounds containing common topologies such that they were largely well modeled in the earlier version of CGenFF. Testing a set of molecules with a lower analogy with CGenFF 2.5.1 could show additional improvements in version 5.0, which we leave for future studies.

<table><caption>**Table 7.** Statistics for FDA molecule parameters generated with CGenFF 2.5.1 (initial) and 5.0 (retrained), showing the agreement with respect to selected QM target data.</caption><thead><tr><th></th><th>Initial</th><th>Retrained</th></tr></thead><tbody><tr><td>Average Charge Penalty</td><td>88</td><td>39</td></tr><tr><td>Average Bonded Penalty</td><td>47</td><td>18</td></tr><tr><td>Water Interaction Energy RMSD (kcal · mol<sup>-1</sup>)</td><td>4.75</td><td>4.13</td></tr></tbody></table>
<table><tr><td>Water Interaction Distance RMSD (Å)</td><td>1.02</td><td>1.02</td></tr><tr><td>Dipole Moment Magnitude Absolute Difference (Debye)</td><td>1.69</td><td>1.29</td></tr><tr><td>Dipole Moment Angle (°)</td><td>32</td><td>28</td></tr><tr><td>Global Minimum Coordinate RMSD (Å)</td><td>0.19</td><td>0.21</td></tr><tr><td>Multiple Conformations Energy RMSD (kcal · mol<sup>-1</sup>)</td><td>3.08</td><td>3.04</td></tr></table>

## Pure solvent properties

The 1,390 optimized molecules utilized for retraining the CGenFF program did not include new atom types, therefore no additional LJ parameters were introduced into the program. Nevertheless, given that the retrained CGenFF program will affect some of the parameters assigned to compounds and that the total collection of parameters for a molecule significantly influences its properties, the pure solvent properties such as heat of vaporization and molecular volume were assessed for the retrained CGenFF. A set of 377 molecules previously employed to validate the LJ parameters for Drude FF⁸⁵ were used to evaluate the performance of the extended version of CGenFF. The range of charge and bonded parameter penalties for these molecules with the new version of CGenFF ranged from 0 to 264, indicating that for a number of molecules parameters transferred by analogy were being used in these molecules. Table 8 presents a statistical analysis of molecular volumes and heat of vaporization derived from pure solvent MD simulations. It is evident that both the original and retrained CGenFF versions adequately reproduce the pure solvent properties, with the retrained CGenFF showing marginally better performance. The MAE of approximately 2.8% for molecular volume and ~9% for heat of vaporization are slightly poorer than those reported in the 2010 CGenFF paper (2.1% and 7% for 111 and 95 compounds, respectively)²⁶ and in the LJ parameter optimization work for the Drude polarizable force field (2% and 5.8%).⁸⁵ Despite the lack of optimization, the parameters generated by CGenFF exhibit overall good quality relative to experimental data, as depicted in Figure 9, with an R² value exceeding 0.8 for the heat of vaporization. This analysis indicates the high quality of the earlier version of CGenFF for the treatment of pure solvent properties with the new version of CGenFF showing small improvements due to the present parametrization effort. This is again associated with the compounds used for the pure solvent properties including the most common types of functional groups and which are already being accurately treated by CGenFF 2.5.1.

<table><caption>**Table 8.** Comparison of CGenFF 2.5.1 and 5.0 for predicting experimental properties of pure solvents across 377 molecules. Values in parentheses represent the differences weighted by the experimental value.</caption><thead><tr><th rowspan="2">Property</th><th colspan="2">Molecular Volume<br>(Å<sup>3</sup>)</th><th colspan="2">Density<br>(g · cm<sup>-3</sup>)</th><th colspan="2">Heat of Vaporization<br>(kcal · mol<sup>-1</sup>)</th></tr><tr><th>2.5.1</th><th>5.0</th><th>2.5.1</th><th>5.0</th><th>2.5.1</th><th>5.0</th></tr></thead><tbody><tr><td>MAE</td><td>5.3 (2.80%)</td><td>5.2 (2.77%)</td><td>0.03 (2.71%)</td><td>0.03 (2.68%)</td><td>1.0 (9.1%)</td><td>0.9 (8.7%)</td></tr><tr><td>RMSD</td><td>7.5 (4.07%)</td><td>7.3 (4.05%)</td><td>0.04 (3.75%)</td><td>0.04 (3.75%)</td><td>1.4 (13.9%)</td><td>1.4 (13.2%)</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Molecular volumes (left panel) and heat of vaporization (right panel) from pure solvent simulations of 377 molecules. The experimental molecular volume of each molecule was derived by dividing its molecular mass by its experimental density.</figcaption></figure>

Free energies of hydration

A subset of the compounds subjected to pure solvent simulations was utilized for HFE calculations
provided availability of experimental HFE and if the compound was neutral and did not contain lone pairs (e.g. Cl
or Br atoms) as these are not implemented in the FEP method in OpenMM and OpenMMTools.⁸⁸,⁹⁵ This yielded
a total of 314 molecules. Comparison of calculated HFE with experimental data shows an R² value of 0.82 for the
initial CGenFF v2.5.1 and 0.84 for retrained CGenFF v5.0 (Figure 10) indicating the high quality of the parameters
generated by both versions of CGenFF for adequately reproducing the experimental HFE data. MAE and RMSD
of the retrained CGenFF v5.0 parameters (1.24 and 1.54 kcal·mol⁻¹, respectively) are also better than the initial
CGenFF parameters (1.30 and 1.63 kcal·mol⁻¹, respectively) with a total of 176 of the 314 molecules having
improved agreement with v5.0. Three individual molecules, 2,4-dimethyl-3-pentanone, cyclopentanone and
tetrahydroinden, show improvement of over 2 kcal·mol⁻¹ in HFE value with CGenFF v5.0 versus v2.5.1 (Table
S9). Their topologies, as reported in Table S10, show the partial charges of atoms while bonded parameters have
also changed for cyclopentanone and tetrahydroinden. Notably, charge as well as parameter penalty values have
also improved in CGenFF v5.0 for these molecules further reflecting improved coverage of chemical space in the
new version. Overall, these results show that the combination of improved electrostatic and bonded parameters
lead to a small, but notable improvement with the optimized v5.0 of CGenFF. The high reproducibility of both the
experimental pure solvent properties and HFEs indicates a good balance of electrostatic and LJ parameters in both
versions of CGenFF.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparison of experimental HFE values with calculated HFE values for 314 molecules using initial and retrained CGenFF parameters.</figcaption></figure>

## Penalty Analysis

Central to CGenFF is the generation of penalties based on analogy assigned to new parameters, including both bonded and partial atomic charges. The penalties are based on individual rules for the bonded terms and partial atomic charges that relate the extent of analogy of a new parameter to available parameters to the magnitude of the penalty.²⁶ Accordingly, the availability of a greater number of parameters available in CGenFF v5.0 will lead to an increase in the number of parameters assigned a zero penalty and a decrease in the magnitude of the penalties assigned to parameters for new terms. This includes the charge penalties for which the bond, angle, and dihedral charge increments are based on bonded parameters already in CGenFF. For the new parameters that are assigned penalties it is anticipated that the magnitude of the penalties will decrease in CGenFF v5.0 given the likelihood that highly analogous atom groups are available in the extended FF. To verify this a subset of ~10K compounds curated from the Maybridge HitFinder Collection of the ZINC database, which was used previously to benchmark the CGenFF program,²⁴ was used to further validate the new CGenFF version in the current work using both CGenFF v2.5.1 and v5.0. As shown in Figure 11 the expected decrease in the penalty scores is observed. With the charge penalties the number of molecules with high penalties is much lower with the overall decrease evident. Similar trends occur with the bonded parameters though to a lesser extent. The number of ZINC molecules for which all bonded parameters have zero penalties increased from 164 for v2.5.1 to 454 for v5.0. Overall, the lower penalties indicate the improved coverage of parameter space associated with the large increase in the number of training set compounds in v5.0 although this does not necessarily indicate significantly improved accuracy, as evidenced with the small changes in the metrics analyzed for the FDA compounds (Table 7). However, as
discussed above this is associated with those compounds being largely covered by the parameters available in v2.5.1 associated with their more common chemical structures. Greater improvements in accuracy may be anticipated with more diverse chemical connectivities as seen with the results for the 1,390 training set compounds reported above in Tables 5 and 6.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Log plots of the probability of atomic penalty scores of the compounds from the ZINC 10K set. Distributions are shown for the (A) charge, (B) bond, (C) valence angle, (D) dihedral angle and (E) improper torsion penalties. A penalty of 10, which was suggested as the cutoff for good quality parameters is indicated by a red dashed line. For bonded parameters, with more parameters being parametrized, fewer parameters will be generated through analogy. In order to make a fair comparison, only bonded parameters generated by version 5.0 for the 10K ZINC validation set are considered when collecting penalty scores from version 2.5.1.</figcaption></figure>

## Conclusion

The presented work significantly extends the training set of molecules and associated parameters used in CGenFF. Model compounds were chosen to represent atomic connectivities new to CGenFF, thus maximizing the number of new parameters included in the training set and increasing the force field's explicit coverage of parameters in the chemical space that CGenFF encompasses. The large-scale parameter optimization followed a hierarchical approach where parameters in molecules with the largest number of new parameters were optimized first. Once a parameter was optimized, it was transferred to other molecules, consistent with the CGenFF
methodology. Target data included QM results for geometries, potential energy scans, interactions of model compounds with water, dipole moments, and electrostatic potentials. The optimization led to significant improvements for the 1,390 molecules studied. In addition, CGenFF coverage was expanded to include a wider range of 4-membered rings including optimization of the new bonded parameters in those molecules. The extended set of parameters was then used to retrain the bond-angle-dihedral increment charging algorithm in the CGenFF program, resulting in CGenFF version 5.0. The new CGenFF version showed improved performance for the 1,390 compounds with smaller improvements with drug-like molecules from the FDA versus v2.5.1. In addition, pure solvent properties and HFE of small molecules showed relatively minor improvements, highlighting the inherent quality of CGenFF for modeling condensed phase properties. The smaller differences in the improvements with the FDA compounds as well as in the pure solvent and HFE results are likely due to those chemical connectivities being largely covered in CGenFF v2.5.1 along with the maintenance of the previously optimized LJ parameters. We anticipate that the improvements seen with the 1,390 training set compounds with this extension of CGenFF will be more pronounced with compounds with more diverse chemical connectivities than those common in approved drugs and small model compounds previously subjected to experimental analyses.

## Supporting Information

Additional figures and tables for: water-compound interactions; pure solvent MD simulations; partial charge distribution; penalty and charge correlation; example of equilibrium values for a molecule with optimized bonded terms; 4-membered ring validation molecules and corresponding comparison of MM vs QM internal coordinates; PES scans energy plots; vibrational analysis; additional data for FDA compounds validation and HFE calculations.

The GitHub repository (https://github.com/mackerell-lab/cgenff5) contains topology and parameter files for the 1,390 molecules for which optimization was performed and developer program; QM and MM coordinates, 2D representations and names for all the molecular datasets; magnitude, and angles of the molecular dipole moment for each molecule, individual empirical and QM water interaction energies and distances for compounds with optimized parameters and the FDA compounds; PES scan plots for all optimized bonded terms; all data points used to generate the figures and tables of the article. Additional data can be provided upon request.

For academic users with an .edu email address, individual molecules may be uploaded online into CGenFF at https://cgenff siksbio.com/ and the CGenFF program can be obtained at no charge for academic users from Silcsbio LLC (silcsbio.com).

## Acknowledgements

This work was supported by grants ANR-18-CE44-0002 to AA, and NIH R35GM131710 to ADM. This project was provided with computer and storage resources by GENCI at CINES thanks to the grants 2023-A0150814640 and 2023-AD010814466 on the supercomputer Adastra.
Conflict of Interest

ADM is co-founder and CSO of SilcsBio LLC.

References

(1) Mackerell, A. D., Jr. Empirical force fields for biological macromolecules: Overview and issues. *J Comput Chem* **2004**, *25* (13), 1584-1604. DOI: 10.1002/jcc.20082.

(2) Yu, Y.; Kramer, A.; Venable, R. M.; Simmonett, A. C.; MacKerell, A. D., Jr.; Klauda, J. B.; Pastor, R. W.; Brooks, B. R. Semi-automated Optimization of the CHARMM36 Lipid Force Field to Include Explicit Treatment of Long-Range Dispersion. *J Chem Theory Comput* **2021**, *17* (3), 1562-1580. DOI: 10.1021/acs.jctc.0c01326.

(3) Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B. L.; Grubmuller, H.; MacKerell, A. D., Jr. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. *Nat Methods* **2017**, *14* (1), 71-73. DOI: 10.1038/nmeth.4067.

(4) Hart, K.; Foloppe, N.; Baker, C. M.; Denning, E. J.; Nilsson, L.; Mackerell, A. D., Jr. Optimization of the CHARMM additive force field for DNA: Improved treatment of the BI/BII conformational equilibrium. *J Chem Theory Comput* **2012**, *8* (1), 348-362. DOI: 10.1021/ct200723y.

(5) Guvench, O.; Mallajosyula, S. S.; Raman, E. P.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T. J.; Jamison, F. W., 2nd; Mackerell, A. D., Jr. CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling. *J Chem Theory Comput* **2011**, *7* (10), 3162-3180. DOI: 10.1021/ct200328p.

(6) Debiec, K. T.; Cerutti, D. S.; Baker, L. R.; Gronenborn, A. M.; Case, D. A.; Chong, L. T. Further along the Road Less Traveled: AMBER ff15ipq, an Original Protein Force Field Built on a Self-Consistent Physical Model. *J Chem Theory Comput* **2016**, *12* (8), 3926-3947. DOI: 10.1021/acs.jctc.6b00567.

(7) 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* (9), 1157-1174. DOI: 10.1002/jcc.20035.

(8) Yildirim, I.; Stern, H. A.; Kennedy, S. D.; Tubbs, J. D.; Turner, D. H. Reparameterization of RNA χ Torsion Parameters for the AMBER Force Field and Comparison to NMR Spectra for Cytidine and Uridine. *J Chem Theory Comput* **2010**, *6*, 1520-1531. DOI: 10.1021/ct900604a.

(9) Dodda, L. S.; Cabeza de Vaca, I.; Tirado-Rives, J.; Jorgensen, W. L. LigParGen web server: an automatic OPLS-AA parameter generator for organic ligands. *Nucleic Acids Res* **2017**, *45* (W1), W331-W336. DOI: 10.1093/nar/gkx312.

(10) Robertson, M. J.; Tirado-Rives, J.; Jorgensen, W. L. Improved Peptide and Protein Torsional Energetics with the OPLS-AA Force Field. *J Chem Theory Comput* **2015**, *11* (7), 3499-3509. DOI: 10.1021/acs.jctc.5b00356.

(11) Kony, D.; Damm, W.; Stoll, S.; van Gunsteren, W. F. An Improved OPLS-AA Force Field for Carbohydrates. *J Comput Chem* **2002**, *23*, 1416-1429. DOI: 10.1002/jcc.10139.

(12) Horta, B. A.; Merz, P. T.; Fuchs, P. F.; Dolenc, J.; Riniker, S.; Hunenberger, P. H. A GROMOS-Compatible Force Field for Small Organic Molecules in the Condensed Phase: The 2016H66 Parameter Set. *J Chem Theory Comput* **2016**, *12* (8), 3825-3850. DOI: 10.1021/acs.jctc.6b00187.

(13) Reif, M. M.; Winger, M.; Oostenbrink, C. Testing of the GROMOS Force-Field Parameter Set 54A8: Structural Properties of Electrolyte Solutions, Lipid Bilayers, and Proteins. *J Chem Theory Comput* **2013**, *9* (2), 1247-1264. DOI: 10.1021/ct300874c.

(14) Pol-Fachin, L.; Rusu, V. H.; Verli, H.; Lins, R. D. GROMOS 53A6GLYC, an Improved GROMOS Force Field for Hexopyranose-Based Carbohydrates. *J Chem Theory Comput* **2012**, *8* (11), 4681-4690. DOI: 10.1021/ct300479h.

(15) Liu, C.; Piquemal, J. P.; Ren, P. AMOEBA+ Classical Potential for Modeling Molecular Interactions. *J Chem Theory Comput* **2019**, *15* (7), 4122-4139. DOI: 10.1021/acs.jctc.9b00261.

(16) Zhang, C.; Lu, C.; Jing, Z.; Wu, C.; Piquemal, J.-P.; Ponder, J. W.; Ren, P. AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids. *J Chem Theory Comput* **2018**, *14* (4), 2084-2108. DOI: 10.1021/acs.jctc.7b01169.

(17) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. *J Chem Theory Comput* **2013**, *9* (9), 4046-4063. DOI: 10.1021/ct4003702.

(18) Nan, Y.; MacKerell, A. D., Jr. Balancing Group I Monatomic Ion-Polar Compound Interactions for Condensed Phase Simulation in the Polarizable Drude Force Field. *J Chem Theory Comput* **2024**, *20* (8), 3242-3257. DOI: 10.1021/acs.jctc.3c01380.

(19) Yu, Y.; Venable, R. M.; Thirman, J.; Chatterjee, P.; Kumar, A.; Pastor, R. W.; Roux, B.; MacKerell, A. D., Jr.; Klauda, J. B. Drude Polarizable Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Saturated and Monounsaturated Zwitterionic Lipids. *J Chem Theory Comput* **2023**, *19* (9), 2590-2605. DOI: 10.1021/acs.jctc.3c00203.

(20) Lin, F. Y.; Huang, J.; Pandey, P.; Rupakheti, C.; Li, J.; Roux, B. T.; MacKerell, A. D., Jr. Further Optimization and Validation of the Classical Drude Polarizable Protein Force Field. *J Chem Theory Comput* **2020**, *16* (5), 3221-3239. DOI: 10.1021/acs.jctc.0c00057.

(21) Pandey, P.; Aytenfisu, A. H.; MacKerell, A. D., Jr.; Mallajosyula, S. S. Drude Polarizable Force Field Parametrization of Carboxylate and N-Acetyl Amine Carbohydrate Derivatives. *J Chem Theory Comput* **2019**, *15* (9), 4982-5000. DOI: 10.1021/acs.jctc.9b00327.

(22) Lemkul, J. A.; MacKerell, A. D., Jr. Polarizable force field for RNA based on the classical drude oscillator. *J Comput Chem* **2018**, *39*, 2624-2646. DOI: 10.1002/jcc.25709.

(23) Song, L.; Lee, T.-S.; Zhu, C.; York, D. M.; Merz Jr., K. M. Validation of AMBER/GAFF for Relative Free Energy Calculations. *ChemRxiv* [Preprint] **2019**. DOI: 10.26434/chemrxiv.7653434.v1.
(24) 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* (12), 3155-3168. DOI: 10.1021/ci3003649.

(25) 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* (12), 3144-3154. DOI: 10.1021/ci300363c.

(26) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. *J Comput Chem* **2010**, *31* (4), 671-690. DOI: 10.1002/jcc.21367.

(27) Halgren, T. A. MMFF VII. Characterization of MMFF94, MMFF94s, and Other Widely Available Force Fields for Conformational Energies and for Intermolecular-Interaction Energies and Geometries. *J Comput Chem* **1999**, *20*, 730-748.

(28) Halgren, T. A. Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. *J Comput Chem* **1996**, *17* (5-6), 490-519. DOI: https://doi.org/10.1002/(SICI)1096-987X(199604)17:5.

(29) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L.; et al. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. *J Chem Theory Comput* **2016**, *12* (1), 281-296. DOI: 10.1021/acs.jctc.5b00864.

(30) Boothroyd, S.; Behara, P. K.; Madin, O. C.; Hahn, D. F.; Jang, H.; Gapsys, V.; Wagner, J. R.; Horton, J. T.; Dotson, D. L.; Thompson, M. W.; et al. Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field. *J Chem Theory Comput* **2023**, *19* (11), 3251-3275. DOI: 10.1021/acs.jctc.3c00039.

(31) Qiu, Y.; Smith, D. G. A.; Boothroyd, S.; Jang, H.; Hahn, D. F.; Wagner, J.; Bannan, C. C.; Gokey, T.; Lim, V. T.; Stern, C. D.; et al. Development and Benchmarking of Open Force Field v1.0.0-the Parsley Small-Molecule Force Field. *J Chem Theory Comput* **2021**, *17* (10), 6262-6280. DOI: 10.1021/acs.jctc.1c00571.

(32) Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1, A data set of 20 million calculated off-equilibrium conformations for organic molecules. *Sci Data* **2017**, *4*, 170193, Data Descriptor. DOI: 10.1038/sdata.2017.193 https://www.nature.com/articles/sdata2017193#supplementary-information.

(33) Botu, V.; Batra, R.; Chapman, J.; Ramprasad, R. Machine Learning Force Fields: Construction, Validation, and Outlook. *J Phys Chem C* **2017**, *121* (1), 511-522. DOI: 10.1021/acs.jpcc.6b10908.

(34) Chmiela, S.; Vassilev-Galindo, V.; Unke, O. T.; Kabylda, A.; Sauceda, H. E.; Tkatchenko, A.; Muller, K. R. Accurate global machine learning force fields for molecules with hundreds of atoms. *Sci Adv* **2023**, *9* (2), eadf0873. DOI: 10.1126/sciadv.adf0873.

(35) Zubatyuk, R.; Smith, J. S.; Leszczynski, J.; Isayev, O. Accurate and transferable multitask prediction of chemical properties with an atoms-in-molecules neural network. *Sci Adv* **2019**, *5* (8), eaav6490. DOI: 10.1126/sciadv.eav6490.

(36) Anstine, D. M.; Isayev, O. Generative Models as an Emerging Paradigm in the Chemical Sciences. *J Am Chem Soc* **2023**, *145* (16), 8736-8750. DOI: 10.1021/jacs.2c13467.

(37) Gastegger, M.; Behler, J.; Marquetand, P. Machine learning molecular dynamics for the simulation of infrared spectra. *Chem Sci* **2017**, *8* (10), 6924-6935, 10.1039/C7SC02267K. DOI: 10.1039/C7SC02267K.

(38) Marchand, J. R.; Caflisch, A. In silico fragment-based drug design with SEED. *Eur J Med Chem* **2018**, *156*, 907-917. DOI: 10.1016/j.ejmech.2018.07.042.

(39) Goossens, K.; Wroblowski, B.; Langini, C.; van Vlijmen, H.; Caflisch, A.; De Winter, H. Assessment of the Fragment Docking Program SEED. *J Chem Inf Model* **2020**, *60* (10), 4881-4893. DOI: 10.1021/acs.jcim.0c00556.

(40) Guterres, H.; Im, W. CHARMM-GUI-Based Induced Fit Docking Workflow to Generate Reliable Protein-Ligand Binding Modes. *J Chem Inf Model* **2023**, *63* (15), 4772-4779. DOI: 10.1021/acs.jcim.3c00416.

(41) Das, S.; Shimshi, M.; Raz, K.; Nitoker Eliaz, N.; Mhashal, A. R.; Ansbacher, T.; Major, D. T. EnzyDock: Protein-Ligand Docking of Multiple Reactive States along a Reaction Coordinate in Enzymes. *J Chem Theory Comput* **2019**, *15* (9), 5116-5134. DOI: 10.1021/acs.jctc.9b00366.

(42) Suh, D.; Schwartz, R.; Gupta, P. K.; Zev, S.; Major, D. T.; Im, W. CHARMM-GUI EnzyDocker for Protein-Ligand Docking of Multiple Reactive States along a Reaction Coordinate in Enzymes. *J Chem Theory Comput* **2025**. DOI: 10.1021/acs.jctc.4c01691.

(43) Ahmed, F.; Brooks, C. L., 3rd. FASTDock: A Pipeline for Allosteric Drug Discovery. *J Chem Inf Model* **2023**, *63* (22), 7219-7227. DOI: 10.1021/acs.jcim.3c00895.

(44) Ding, X.; Wu, Y.; Wang, Y.; Vilseck, J. Z.; Brooks, C. L., 3rd. Accelerated CDOCKER with GPUs, Parallel Simulated Annealing, and Fast Fourier Transforms. *J Chem Theory Comput* **2020**, *16* (6), 3910-3919. DOI: 10.1021/acs.jctc.0c00145.

(45) Vilseck, J. Z.; Cervantes, L. F.; Hayes, R. L.; Brooks, C. L., 3rd. Optimizing Multisite lambda-Dynamics Throughput with Charge Renormalization. *J Chem Inf Model* **2022**, *62* (6), 1479-1488. DOI: 10.1021/acs.jcim.2c00047.

(46) Ustach, V. D.; Lakkaraju, W. K.; Jo, S.; Yu, W.; Jiang, W.; MacKerell, A. D. Optimization and Evaluation of Site-Identification by Ligand Competitive Saturation (SILCS) as a Tool for Target-Based Ligand Optimization. *J. Chem. Inf. Model.* **2019**, *59* (6), 3018-3035. DOI: 10.1021/acs.jcim.9b00210 (accessed 2023/02/15/14:24:16). From DOI.org (Crossref).

(47) Kim, S.; Lee, J.; Jo, S.; Brooks, C. L., 3rd; Lee, H. S.; Im, W. CHARMM-GUI ligand reader and modeler for CHARMM force field generation of small molecules. *J Comput Chem* **2017**, *38* (21), 1879-1886. DOI: 10.1002/jcc.24829.

(48) Spiliotopoulos, D.; Wamhoff, E. C.; Lolli, G.; Rademacher, C.; Caflisch, A. Discovery of BAZ2A bromodomain ligands. *Eur J Med Chem* **2017**, *139*, 564-572. DOI: 10.1016/j.ejmech.2017.08.028.

(49) Mussakhmetov, A.; Shumilin, I. A.; Nugmanova, R.; Shabalin, I. G.; Baizhumanov, T.; Toibazar, D.; Khassenov, B.; Minor, W.; Utepbergenov, D. A transient post-translational modification of active site cysteine alters binding properties of the parkinsonism protein DJ-1. *Biochemical and Biophysical Research Communications* **2018**, *504* (1), 328-333. DOI: 10.1016/j.bbrc.2018.08.190.

(50) Dalle Vedove, A.; Cazzanelli, G.; Corsi, J.; Sedykh, M.; D'Agostino, V. G.; Caflisch, A.; Lolli, G. Identification of a BAZ2A Bromodomain Hit Compound by Fragment Joining. *ACS Bio Med Chem Au* **2021**, *1* (1), 5-10. DOI: 10.1021/acsbmedchemau.1c00016.
(51) Dalle Vedove, A.; Cazzanelli, G.; Batiste, L.; Marchand, J. R.; Spiliotopoulos, D.; Corsi, J.; D'Agostino, V. G.; Caflisch, A.; Lolli, G. Identification of a BAZ2A-Bromodomain Hit Compound by Fragment Growing. *ACS Med Chem Lett* **2022**, *13* (9), 1434-1443. DOI: 10.1021/acsmedchemlett.2c00173.

(52) Nai, F.; Flores Espinoza, M. P.; Invernizzi, A.; Vargas-Rosales, P. A.; Bobileva, O.; Herok, M.; Caflisch, A. Small-Molecule Inhibitors of the m7G-RNA Writer METTL1. *ACS Bio Med Chem Au* **2024**, *4* (2), 100-110. DOI: 10.1021/acsbiomedchemau.3c00030.

(53) Li, Y.; Bedi, R. K.; Nai, F.; von Roten, V.; Dolbois, A.; Zálešák, F.; Nachawati, R.; Huang, D.; Caflisch, A. Structure-based design of ligands of the m6A-RNA reader YTHDC1. *European Journal of Medicinal Chemistry Reports* **2022**, *5*. DOI: 10.1016/j.ejmcr.2022.100057.

(54) Peiffer, A. L.; Garlick, J. M.; Wu, Y.; Wotring, J. W.; Arora, S.; Harmata, A. S.; Bochar, D. A.; Stephenson, C. J.; Soellner, M. B.; Sexton, J. Z.; et al. TMPRSS2 Inhibitor Discovery Facilitated through an In Silico and Biochemical Screening Platform. *ACS Med Chem Lett* **2023**, *14* (6), 860-866. DOI: 10.1021/acsmedchemlett.3c00035.

(55) Kumar, A.; MacKerell, A. D., Jr. FFParam-v2.0: A Comprehensive Tool for CHARMM Additive and Drude Polarizable Force-Field Parameter Optimization and Validation. *J Phys Chem B* **2024**, *128* (18), 4385-4395. DOI: 10.1021/acs.jpcb.4c01314.

(56) Vanommeslaeghe, K.; MacKerell, A. D., Jr. CHARMM additive and polarizable force fields for biophysics and computer-aided drug design. *Biochim Biophys Acta* **2015**, *1850* (5), 861-871. DOI: 10.1016/j.bbgagen.2014.08.004.

(57) Allen, M. P.; Tildesley, D. J. *Computer Simulation of Liquids: Second Edition*; Oxford University Press, 2017.

(58) Irwin, J. J.; Tang, K. G.; Young, J.; Dandarchuluun, C.; Wong, B. R.; Khurelbaatar, M.; Moroz, Y. S.; Mayfield, J.; Sayle, R. A. ZINC20—A Free Ultralarge-Scale Chemical Database for Ligand Discovery. *J. Chem. Inf. Model.* **2020**, *60* (12), 6065-6073. DOI: 10.1021/acs.jcim.0c00675 (accessed 2021/05/12/15:33:02). From ACS Publications.

(59) *MarvinSketch*; Chemaxon: Hungary, 2019. https://chemaxon.com/products/marvin (accessed 2021/02/22/).

(60) *RDKit: Open-source cheminformatics*; California, US, 2019. http://www.rdkit.org/ (accessed 2022/07/20/08:56:23).

(61) Hanwell, M. D.; Curtis, D. E.; Lonie, D. C.; Vandermeersch, T.; Zurek, E.; Hutchison, G. R. Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. *J Cheminform* **2012**, *4* (1), 17. DOI: 10.1186/1758-2946-4-17.

(62) O'Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An open chemical toolbox. *J Cheminform* **2011**, *3* (1), 33. DOI: 10.1186/1758-2946-3-33.

(63) Xu, Y.; Vanommeslaeghe, K.; Aleksandrov, A.; MacKerell, A. D., Jr.; Nilsson, L. Additive CHARMM force field for naturally occurring modified ribonucleotides. *J Comput Chem* **2016**, *37* (10), 896-912. DOI: 10.1002/jcc.24307.

(64) Croitoru, A.; Park, S.-J.; Kumar, A.; Lee, J.; Im, W.; MacKerell, A. D., Jr.; Aleksandrov, A. Additive CHARMM36 Force Field for Nonstandard Amino Acids. *J Chem Theory Comput* **2021**, *17* (6), 3554-3570. DOI: 10.1021/acs.jctc.1c00254.

(65) Pracht, P.; Bohle, F.; Grimme, S. Automated exploration of the low-energy chemical space with fast quantum chemical methods. *Phys Chem Chem Phys* **2020**, *22* (14), 7169-7192. DOI: 10.1039/c9cp06869d.

(66) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; et al. Psi4: an open-source ab initio electronic structure program. *WIREs Comput Mol Sci* **2012**, *2* (4), 556-565. DOI: 10.1002/wcm9.93.

(67) Foresman, J. B.; Frisch, A. *Exploring chemistry with electronic structure methods*; Gaussian, Inc, 2015.

(68) Head-Gordon, M.; Pople, J. A.; Frisch, M. J. MP2 energy evaluation by direct methods. *Chem Phys Lett* **1988**, *153* (6), 503-506. DOI: 10.1016/0009-2614(88)85250-3.

(69) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions. *J Chem Phys* **1980**, *72* (1), 650-654. DOI: 10.1063/1.438955.

(70) Huang, L.; Roux, B. Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data. *J Chem Theory Comput* **2013**, *9* (8), 3543-3556. DOI: 10.1021/ct4003477.

(71) *Gaussian*; Gaussian Inc, Wallingford CT: 2009. https://gaussian.com/ (accessed 2019/02/22).

(72) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. *J Phys Chem B* **1998**, *102* (18), 3586-3616. DOI: 10.1021/jp973084f.

(73) Jorgensen, W. L. Transferable Intermolecular Potential Functions for Waters, Alcohols, and Ethers. Application to Liquid Water. *J. Am. Chem. Soc.* **1981**, *103*, 335.

(74) Boys, S. F.; Bernardi, F. The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. *Mol Phys* **1970**, *19* (4), 553-566, Molecular Physics. DOI: 10.1080/00268977000101561.

(75) Croitoru, A.; Aleksandrov, A. Parametrization of Force Field Bonded Terms under Structural Inconsistency. *J Chem Inf Model* **2022**, *62* (19), 4771-4782. DOI: 10.1021/acs.jcim.2c00950.

(76) Scott, A. P.; Radom, L. Harmonic Vibrational Frequencies: An Evaluation of Hartree-Fock, Møller-Plesset, Quadratic Configuration Interaction, Density Functional Theory, and Semiempirical Scale Factors. *J Phys Chem* **1996**, *100* (41), 16502-16513. DOI: 10.1021/jp960976r.

(77) Brooks, B. R.; Brooks, C. L.; MacKerell Jr, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The biomolecular simulation program. *J Comput Chem* **2009**, *30* (10), 1545-1614. DOI: 10.1002/jcc.21287.

(78) Aleksandrov, A. A Molecular Mechanics Model for Flavins. *J Comput Chem* **2019**, *40* (32), 2834-2842. DOI: 10.1002/jcc.26061.

(79) Powell, M. J. D. An efficient method for finding the minimum of a function of several variables without calculating derivatives. *The Computer Journal* **1964**, *7* (2), 155-162. DOI: 10.1093/comjnl/7.2.155.

(80) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. *Journal of Computational Chemistry* **2004**, *24* (2), 187-217. DOI: 10.1002/jcc.540040211.

(81) Smith, D. G. A.; Burns, L. A.; Simmonett, A. C.; Parrish, R. M.; Schieber, M. C.; Galvelis, R.; Kraus, P.; Kruse, H.; Di Remigio, R.; Alenaizan,
A.; et al. Psi4 1.4: Open-source software for high-throughput quantum chemistry. *J Chem Phys* **2020**, *152* (18), 184108. DOI: 10.1063/5.0006002.

(82) Kumar, A.; Yoluk, O.; MacKerell, A. D., Jr. FFParam: Standalone package for CHARMM additive and Drude polarizable force field parametrization of small molecules. *J Comput Chem* **2020**, *41* (9), 958-970. DOI: 10.1002/jcc.26138.

(83) Hwang, W.; Austin, S. L.; Blondel, A.; Boittier, E. D.; Boresch, S.; Buck, M.; Buckner, J.; Caflisch, A.; Chang, H. T.; Cheng, X.; et al. CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed. *J Phys Chem B* **2024**, DOI: 10.1021/acs.jpcb.4c04100.

(84) Kumar, A.; Pandey, P.; Chatterjee, P.; MacKerell, A. D., Jr. Deep Neural Network Model to Predict the Electrostatic Parameters in the Polarizable Classical Drude Oscillator Force Field. *J Chem Theory Comput* **2022**, *18* (3), 1711-1725. DOI: 10.1021/acs.jctc.1c01166.

(85) Rupakheti, C. R.; MacKerell, A. D., Jr.; Roux, B. Global Optimization of the Lennard-Jones Parameters for the Drude Polarizable Force Field. *J Chem Theory Comput* **2021**, *17* (11), 7085-7095. DOI: 10.1021/acs.jctc.1c00664.

(86) Boulanger, E.; Huang, L.; Rupakheti, C.; MacKerell, A. D., Jr.; Roux, B. Optimized Lennard-Jones Parameters for Druglike Small Molecules. *J Chem Theory Comput* **2018**, *14* (6), 3121-3131. DOI: 10.1021/acs.jctc.8b00172.

(87) Chu, J.-W.; Trout, B. L.; Brooks, B. R. A super-linear minimization scheme for the nudged elastic band method. *J Chem Phys* **2003**, *119* (24), 12708-12717. DOI: 10.1063/1.1627754.

(88) Eastman, P.; Galvelis, R.; Peláez, R. P.; Abreu, C. R. A.; Farr, S. E.; Gallicchio, E.; Gorenko, A.; Henry, M. M.; Hu, F.; Huang, J.; et al. OpenMM 8: Molecular Dynamics Simulation with Machine Learning Potentials. *J Phys Chem B* **2024**, *128* (1), 109-116. DOI: 10.1021/acs.jpcb.3c06662.

(89) Izaguirre, J. A.; Sweet, C. R.; Pande, V. S. Multiscale dynamics of macromolecules using normal mode Langevin. *Pac Symp Biocomput* **2010**, 240-251. DOI: 10.1142/9789814295291_0026.

(90) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. *J Chem Phys* **1995**, *103* (19), 8577-8593. DOI: 10.1063/1.470117.

(91) Chow, K.-H.; Ferguson, D. M. Isothermal-isobaric molecular dynamics simulations with Monte Carlo volume sampling. *Comput Phys Commun* **1995**, *91* (1), 283-289. DOI: 10.1016/0010-4655(95)00059-0.

(92) Åqvist, J.; Wennerström, P.; Nervall, M.; Bjelic, S.; Brandsdal, B. O. Molecular dynamics simulations of water and biomolecules with a Monte Carlo constant pressure algorithm. *Chem Phys Lett* **2004**, *384* (4), 288-294. DOI: 10.1016/j.cplett.2003.12.039.

(93) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. *Chem Phys Lett* **1994**, *222* (6), 529-539. DOI: 10.1016/0009-2614(94)00397-1.

(94) Pham, T. T.; Shirts, M. R. Optimal pairwise and non-pairwise alchemical pathways for free energy calculations of molecular transformation in solution phase. *J Chem Phys* **2012**, *136* (12). DOI: 10.1063/1.3697833.

(95) *OpenMMTools*; US, 2023. https://github.com/choderalab/openmmtools (accessed 2024/11/21).

(96) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T. A.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. *J Chem Phys* **1995**, *103*, 8577-8593. DOI: 10.1063/1.470117.

(97) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. *J Comput Phys* **1977**, *23* (3), 327-341. DOI: 10.1016/0021-9991(77)90098-5.

(98) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. *J Chem Phys* **2008**, *129* (12). DOI: 10.1063/1.2978177.

(99) Bhatia, P. A.; Randolph, J. T.; Schrimpf, M. R.; Zhang, Q. I. Diazzepine Derivatives Useful as Calcium Channel Blockers in the Treatment of Pain. US151387105, 2014.

(100) Hiesinger, K.; Dar'in, D.; Proschak, E.; Krasavin, M. Spirocyclic Scaffolds in Medicinal Chemistry. *J Med Chem* **2021**, *64* (1), 150-183. DOI: 10.1021/acs.jmedchem.0c01473.

(101) Teli, S.; Teli, P.; Soni, S.; Sahiba, N.; Agarwal, S. Synthetic aspects of 1,4- and 1,5-benzodiazepines using o-phenylenediamine: a study of past quinquennial. *RSC Adv* **2023**, *13* (6), 3694-3714. DOI: 10.1039/d2ra06045k.
Table of Contents Graphic

1,300 Drug-like molecules

<figure><img src="image_14.png" /></figure>

Validation

- QM

- Experimental

Heat of Vaporization
Hydration Free Energy
Density

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

Parameter
optimization

CGenFF