Exploration and validation of force field design protocols through QM-to-MM mapping

Excitonic Configuration Interaction: Going Beyond the Frenkel Exciton Model

We present the excitonic configuration interaction (ECI) method─a fragment-based analogue of the CI method for electronic structure calculations of multichromophoric systems. It can also be viewed as a generalization of the exciton approach, with the following properties: (i) It constructs the effective Hamiltonian exclusively from monomer calculations. (ii) It employs the strong orthogonality assumption and is exact within McWeeny's group function theory, thus requiring only one-electron density matrices of the monomer states. (iii) It is agnostic of the monomer electronic structure method, allowing us to use/combine different methods. (iv) It includes an embedding point charge scheme (called excitonic Hartree-Fock, EHF) to improve the accuracy of the monomer states, but such that the effective full-system Hamiltonian is not explicitly dependent on the embedding. (v) It is systematically improvable, by expanding the set of monomer states and by including configurations where two or more monomers are excited (defining the ECIS, ECISD, etc., methods). The performance of ECI is assessed by computing the absorption spectrum of two exemplary multichromophoric systems, using CIS as the monomer electronic structure method. The accuracy of ECI significantly depends on the chosen embedding charges and the ECI expansion. The most accurate assessed combinations─ECIS or ECISD with EHF embedding─yield spectra that agree qualitatively and quantitatively with full-system direct calculations, with deviations of the excitation energies below 0.1 eV. We also show that ECISD based on CIS monomer calculations can predict states where two monomers are excited simultaneously (e.g., triplet-triplet double-local excitations) that are inaccessible in a full-system CIS calculation.

Axial H-Bonding Solvent Controls Inhomogeneous Spectral Broadening, While Peripheral H-Bonding Solvent Controls Vibronic Broadening: Cresyl Violet in Methanol

The dynamics of the nuclei of both a chromophore and its condensed-phase environment control many spectral features, including the vibronic and inhomogeneous broadening present in spectral line shapes. For the cresyl violet chromophore in methanol, we here analyze and isolate the effect of specific chromophore-solvent interactions on simulated spectral densities, reorganization energies, and linear absorption spectra. Employing both chromophore and its condensed-phase environment control many spectral features, including the vibronic and inhomogeneous broadening present in spectral line shapes. For the cresyl violet chromophore in methanol, we here analyze and isolate the effect of specific chromophore-solvent interactions on simulated spectral densities, reorganization energies, and linear absorption spectra. Employing both force field and ab initio molecular dynamics trajectories along with the inclusion of only certain solvent molecules in the excited-state calculations, we determine that the methanol molecules axial to the chromophore are responsible for the majority of inhomogeneous broadening, with a single methanol molecule that forms an axial hydrogen bond dominating the response. The strong peripheral hydrogen bonds do not contribute to spectral broadening, as they are very stable throughout the dynamics and do not lead to increased energy-gap fluctuations. We also find that treating the strong peripheral hydrogen bonds as molecular mechanical point charges during the molecular dynamics simulation underestimates the vibronic coupling. Including these peripheral hydrogen bonding methanol molecules in the quantum-mechanical region in a geometry optimization increases the vibronic coupling, suggesting that a more advanced treatment of these strongly interacting solvent molecules during the molecular dynamics trajectory may be necessary to capture the full vibronic spectral broadening.

Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field

We introduce the Open Force Field (OpenFF) 2.0.0 small molecule force field for drug-like molecules, code-named Sage, which builds upon our previous iteration, Parsley. OpenFF force fields are based on direct chemical perception, which generalizes easily to highly diverse sets of chemistries based on substructure queries. Like the previous OpenFF iterations, the Sage generation of OpenFF force fields was validated in protein-ligand simulations to be compatible with AMBER biopolymer force fields. In this work, we detail the methodology used to develop this force field, as well as the innovations and improvements introduced since the release of Parsley 1.0.0. One particularly significant feature of Sage is a set of improved Lennard-Jones (LJ) parameters retrained against condensed phase mixture data, the first refit of LJ parameters in the OpenFF small molecule force field line. Sage also includes valence parameters refit to a larger database of quantum chemical calculations than previous versions, as well as improvements in how this fitting is performed. Force field benchmarks show improvements in general metrics of performance against quantum chemistry reference data such as root-mean-square deviations (RMSD) of optimized conformer geometries, torsion fingerprint deviations (TFD), and improved relative conformer energetics (ΔΔE). We present a variety of benchmarks for these metrics against our previous force fields as well as in some cases other small molecule force fields. Sage also demonstrates improved performance in estimating physical properties, including comparison against experimental data from various thermodynamic databases for small molecule properties such as ΔHmix, ρ(x), ΔGsolv, and ΔGtrans. Additionally, we benchmarked against protein-ligand binding free energies (ΔGbind), where Sage yields results statistically similar to previous force fields. All the data is made publicly available along with complete details on how to reproduce the training results at https://github.com/openforcefield/openff-sage.

What is the Optimal Dipole Moment for Nonpolarizable Models of Liquids?

In classical nonpolarizable models, electrostatic interactions are usually described by assigning fixed partial charges to interaction sites. Despite the multitude of methods and theories proposed over the years for partial charge assignment, a fundamental question remains─what is the correct degree of polarization that a fixed-charge model should possess to provide the best balance of interactions (including induction effects) and yield the best description of the potential energy surface of a liquid phase? We address this question by approaching it from two separate and independent viewpoints: the QUantum mechanical BEspoke (QUBE) approach, which assigns bespoke force field parameters for individual molecules from ab initio calculations with minimal empirical fitting, and the Polarization-Consistent Approach (PolCA) force field, based on empirical fitting of force field parameters with an emphasis on transferability by rigorously accounting for polarization effects in the parameterization process. We show that the two approaches yield consistent answers to the above question, namely, that the dipole moment of the model should be approximately halfway between those of the gas and the liquid phase. Crucially, however, the reference liquid-phase dipole needs to be estimated using methods that explicitly consider both mean-field and local contributions to polarization. In particular, continuum dielectric models are inadequate for this purpose because they cannot account for local effects and therefore significantly underestimate the degree of polarization of the molecule. These observations have profound consequences for the development, validation, and testing of nonpolarizable models.

Open Force Field BespokeFit: Automating Bespoke Torsion Parametrization at Scale

The development of accurate transferable force fields is key to realizing the full potential of atomistic modeling in the study of biological processes such as protein-ligand binding for drug discovery. State-of-the-art transferable force fields, such as those produced by the Open Force Field Initiative, use modern software engineering and automation techniques to yield accuracy improvements. However, force field torsion parameters, which must account for many stereoelectronic and steric effects, are considered to be less transferable than other force field parameters and are therefore often targets for bespoke parametrization. Here, we present the Open Force Field QCSubmit and BespokeFit software packages that, when combined, facilitate the fitting of torsion parameters to quantum mechanical reference data at scale. We demonstrate the use of QCSubmit for simplifying the process of creating and archiving large numbers of quantum chemical calculations, by generating a dataset of 671 torsion scans for druglike fragments. We use BespokeFit to derive individual torsion parameters for each of these molecules, thereby reducing the root-mean-square error in the potential energy surface from 1.1 kcal/mol, using the original transferable force field, to 0.4 kcal/mol using the bespoke version. Furthermore, we employ the bespoke force fields to compute the relative binding free energies of a congeneric series of inhibitors of the TYK2 protein, and demonstrate further improvements in accuracy, compared to the base force field (MUE reduced from 0.56_0.39^0.77 to 0.42_0.28^0.59 kcal/mol and R² correlation improved from 0.72_0.35^0.87 to 0.93_0.84^0.97).