QUBEKit: Automating the Derivation of Force Field Parameters from Quantum Mechanics

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.

System-Specific Parameter Optimization for Nonpolarizable and Polarizable Force Fields

The accuracy of classical force fields (FFs) has been shown to be limited for the simulation of cation-protein systems despite their importance in understanding the processes of life. Improvements can result from optimizing the parameters of classical FFs or by extending the FF formulation by terms describing charge transfer (CT) and polarization (POL) effects. In this work, we introduce our implementation of the CTPOL model in OpenMM, which extends the classical additive FF formula by adding CT and POL. Furthermore, we present an open-source parametrization tool, called FFAFFURR, that enables the (system-specific) parametrization of OPLS-AA and CTPOL models. The performance of our workflow was evaluated by its ability to reproduce quantum chemistry energies and by molecular dynamics simulations of a zinc-finger protein.

Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

AutoParams: An Automated Web-Based Tool To Generate Force Field Parameters for Molecular Dynamics Simulations

Many research questions benefit from molecular dynamics simulations to observe the motions and conformations of molecules over time, which rely on force fields that describe sets of common molecules by category. With the increase of importance for large data sets used in machine learning and growing computational efficiency, the ability to rapidly create large numbers of force field inputs is of high importance. Unusual molecules, such as nucleotide analogues, functionalized carbohydrates, and modified amino acids, are difficult to describe consistently using standard force fields, requiring the development of custom parameters for each unique molecule. While these parameters may be created by individual users, the process can become time-consuming or may introduce errors that may not be immediately apparent. We present an open-source automated parameter generation service, AutoParams, which requires minimal input from the user and creates useful Amber force field parameter sets for most molecules, particularly those that combine molecular types (e.g., a carbohydrate functionalized with a benzene). We include hierarchical atom-typing logic that makes it straightforward to expand with additional force fields and settings, and options for creating monomers in polymers, such as functionalized amino acids. It can be straightforwardly linked to any charge generation program and currently has interfaces to Psi4, PsiRESP, and TeraChem. It is open source and is available via GitHub. It includes error checking and testing protocols to ensure the parameters will be sufficient for subsequent molecular dynamics simulations and streamlines the creation of force field databases.

Evaluation of the MACE force field architecture: From medicinal chemistry to materials science

The MACE architecture represents the state of the art in the field of machine learning force fields for a variety of in-domain, extrapolation, and low-data regime tasks. In this paper, we further evaluate MACE by fitting models for published benchmark datasets. We show that MACE generally outperforms alternatives for a wide range of systems, from amorphous carbon, universal materials modeling, and general small molecule organic chemistry to large molecules and liquid water. We demonstrate the capabilities of the model on tasks ranging from constrained geometry optimization to molecular dynamics simulations and find excellent performance across all tested domains. We show that MACE is very data efficient and can reproduce experimental molecular vibrational spectra when trained on as few as 50 randomly selected reference configurations. We further demonstrate that the strictly local atom-centered model is sufficient for such tasks even in the case of large molecules and weakly interacting molecular assemblies.

Evidence That Less Can Be More for Transferable Force Fields

Graph-based parameter assignment has been the basis for developing transferable force fields for molecular dynamics simulations for decades. Nevertheless, transferable force fields vary in how specifically terms are defined with respect to the molecular graph and the procedures for generating parametrization data. More-specific force-field terms increase the complexity of the force field, theoretically increasing accuracy but also increasing training data requirements. In contrast, less-specific force fields can be reused across larger regions of chemical space, theoretically reducing accuracy but also reducing the number of parameters and training data requirements. Here, the tradeoffs between force-field specificity and accuracy are quantified by parametrizing three new sets of force fields with varying levels of graph specificity, using a shared procedure for generating training data. These force fields are benchmarked for their ability to reproduce the structural features and liquid properties of 87 organic molecules at 146 distinct state points. The overall accuracy for properties that were directly trained on rapidly saturates as the graph specificity of the force-field increases. From this, we conclude there is at best a marginal benefit of using less transferable and more complex force fields with common sources of quantum-chemically derived training data. When looking at properties unseen during training, there is some evidence that the more-complex force fields even perform slightly worse. These results are rationalized by the fortuitous regularization of force fields based on less-specific and more-transferable atom types. Both the saturation in the accuracy of training properties and the marginally worse performance on off-target properties fundamentally contradict the expectation that bespoke force fields are generally more accurate, given their larger number of parameters, and suggests that increasing force-field complexity should be carefully justified against performance gains and balanced against available training data.

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.

Regularized by Physics: Graph Neural Network Parametrized Potentials for the Description of Intermolecular Interactions

Simulations of molecular systems using electronic structure methods are still not feasible for many systems of biological importance. As a result, empirical methods such as force fields (FF) have become an established tool for the simulation of large and complex molecular systems. The parametrization of FF is, however, time-consuming and has traditionally been based on experimental data. Recent years have therefore seen increasing efforts to automatize FF parametrization or to replace FF with machine-learning (ML) based potentials. Here, we propose an alternative strategy to parametrize FF, which makes use of ML and gradient-descent based optimization while retaining a functional form founded in physics. Using a predefined functional form is shown to enable interpretability, robustness, and efficient simulations of large systems over long time scales. To demonstrate the strength of the proposed method, a fixed-charge and a polarizable model are trained on ab initio potential-energy surfaces. Given only information about the constituting elements, the molecular topology, and reference potential energies, the models successfully learn to assign atom types and corresponding FF parameters from scratch. The resulting models and parameters are validated on a wide range of experimentally and computationally derived properties of systems including dimers, pure liquids, and molecular crystals.

Treating Polarization Effects in Charged and Polar Bio-Molecules Through Variable Electrostatic Parameters

Polarization plays important roles in charged and hydrogen bonding containing systems. Much effort ranging from the construction of physics-based models to quantum mechanism (QM)-based and machine learning (ML)-assisted models have been devoted to incorporating the polarization effect into the conventional force fields at different levels, such as atomic and coarse grained (CG). The application of polarizable force fields or polarization models was limited by two aspects, namely, computational cost and transferability. Different from physics-based models, no predetermining parameters were required in the QM-based approaches. Taking advantage of both the accuracy of QM calculations and efficiency of molecular mechanism (MM) and ML, polarization effects could be treated more efficiently while maintaining the QM accuracy. The computational cost could be reduced with variable electrostatic parameters, such as the charge, dipole, and electronic dielectric constant with the help of linear scaling fragmentation-based QM calculations and ML models. Polarization and entropy effects on the prediction of partition coefficient of druglike molecules are demonstrated by using both explicit or implicit all-atom molecular dynamics simulations and machine learning-assisted models. Directions and challenges for future development are also envisioned.

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).

Accurate Quantum-Mechanically Derived Force-Fields through a Fragment-Based Approach: Balancing Specificity and Transferability in the Prediction of Self-Assembly in Soft Matter

The wide range of time/length scales covered by self-assembly in soft matter makes molecular dynamics (MD) the ideal candidate for simulating such a supramolecular phenomenon at an atomistic level. However, the reliability of MD outcomes heavily relies on the accuracy of the adopted force-field (FF). The spontaneous re-ordering in liquid crystalline materials stands as a clear example of such collective self-assembling processes, driven by a subtle and delicate balance between supramolecular interactions and single-molecule flexibility. General-purpose transferable FFs often dramatically fail to reproduce such complex phenomena, for example, the error on the transition temperatures being larger than 100 K. Conversely, quantum-mechanically derived force-fields (QMD-FFs), specifically tailored for the target system, were recently shown (J. Phys. Chem. Lett.2022,13, 243) to allow for the required accuracy as they not only well reproduced transition temperatures but also yielded a quantitative agreement with the experiment on a wealth of structural, dynamic, and thermodynamic properties. The main drawback of this strategy stands in the computational burden connected to the numerous quantum mechanical (QM) calculations usually required for a target-specific parameterization, which has undoubtedly hampered the routine application of QMD-FFs. In this work, we propose a fragment-based strategy to extend the applicability of QMD-FFs, in which the amount of QM calculations is significantly reduced, being a single-molecule-optimized geometry and its Hessian matrix the only QM information required. To validate this route, a new FF is assembled for a large mesogen, exploiting the parameters obtained for two smaller liquid crystalline molecules, in this and previous work. Lengthy MD simulations are carried out with the new transferred QMD-FF, observing again a spontaneous re-orientation in the correct range of temperatures, with good agreement with the available experimental measures. The present results strongly suggest that a partial transfer of QMD-FF parameters can be invoked without a significant loss of accuracy, thus paving the way to exploit the method's intrinsic predictive capabilities in the simulation of novel soft materials.

Efficient Crystal Structure Prediction for Structurally Related Molecules with Accurate and Transferable Tailor-Made Force Fields

Crystal structure prediction (CSP) his generally used to complement experimental solid form screening and applied to individual molecules in drug development. The fast development of algorithms and computing resources offers the opportunity to use CSP earlier and for a broader range of applications in the drug design cycle. This study presents a novel paradigm of CSP specifically designed for structurally related molecules, referred to as Quick-CSP. The approach prioritizes more accurate physics through robust and transferable tailor-made force fields (TMFFs), such that significant efficiency gains are achieved through the reduction of expensive ab initio calculations. The accuracy of the TMFF is increased by the introduction of electrostatic multipoles, and the fragment-based force field parameterization scheme is demonstrated to be transferable for a family of chemically related molecules. The protocol is benchmarked with structurally related compounds from the Bromodomain and Extraterminal (BET) domain inhibitors series. A new convergence criterion is introduced that aims at performing only as many ab initio optimizations of crystal structures as required to locate the bottom of the crystal energy landscape within a user-defined accuracy. The overall approach provides significant cost savings ranging from three- to eight-fold less than the full-CSP workflow. The reported advancements expand the scope and utility of the underlying CSP building blocks as well as their novel reassembly to other applications earlier in the drug design cycle to guide molecule design and selection.

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

The scale of the parameter optimisation problem in traditional molecular mechanics force field construction means that design of a new force field is a long process, and sub-optimal choices made in the early stages can persist for many generations. We hypothesise that careful use of quantum mechanics to inform molecular mechanics parameter derivation (QM-to-MM mapping) should be used to significantly reduce the number of parameters that require fitting to experiment and increase the pace of force field development. Here, we design and train a collection of 15 new protocols for small, organic molecule force field derivation, and test their accuracy against experimental liquid properties. Our best performing model has only seven fitting parameters, yet achieves mean unsigned errors of just 0.031 g cm^-3 and 0.69 kcal mol^-1 in liquid densities and heats of vaporisation, compared to experiment. The software required to derive the designed force fields is freely available at https://github.com/qubekit/QUBEKit.

Molecular Dynamics with Constrained Nuclear Electronic Orbital Density Functional Theory: Accurate Vibrational Spectra from Efficient Incorporation of Nuclear Quantum Effects

Nuclear quantum effects play a crucial role in many chemical and biological systems involving hydrogen atoms yet are difficult to include in practical molecular simulations. In this paper, we combine our recently developed methods of constrained nuclear-electronic orbital density functional theory (cNEO-DFT) and constrained minimized energy surface molecular dynamics (CMES-MD) to create a new method for accurately and efficiently describing nuclear quantum effects in molecular simulations. By use of this new method, dubbed cNEO-MD, the vibrational spectra of a set of small molecules are calculated and compared with those from conventional ab initio molecular dynamics (AIMD) as well as from experiments. With the same formal scaling, cNEO-MD greatly outperforms AIMD in describing the vibrational modes with significant hydrogen motion characters, demonstrating the promise of cNEO-MD for simulating chemical and biological systems with significant nuclear quantum effects.

Recent progress in general force fields of small molecules

Recent advances in computational hardware and free energy algorithms enable a broader application of molecular simulation of binding interactions between receptors and small-molecule ligands. The underlying molecular mechanics force fields (FFs) for small molecules have also achieved advancements in accuracy, user-friendliness, and speed during the past several years (2018-2020). Besides the expansion of chemical space coverage of ligand-like molecules among major popular classical additive FFs and polarizable FFs, new charge models have been proposed for better accuracy and transferability, new chemical perception of avoiding predefined atom types have been applied, and new automated parameterization toolkits, including machine learning approaches, have been developed for users' convenience.

Linear Atomic Cluster Expansion Force Fields for Organic Molecules: Beyond RMSE

We demonstrate that fast and accurate linear force fields can be built for molecules using the atomic cluster expansion (ACE) framework. The ACE models parametrize the potential energy surface in terms of body-ordered symmetric polynomials making the functional form reminiscent of traditional molecular mechanics force fields. We show that the four- or five-body ACE force fields improve on the accuracy of the empirical force fields by up to a factor of 10, reaching the accuracy typical of recently proposed machine-learning-based approaches. We not only show state of the art accuracy and speed on the widely used MD17 and ISO17 benchmark data sets, but we also go beyond RMSE by comparing a number of ML and empirical force fields to ACE on more important tasks such as normal-mode prediction, high-temperature molecular dynamics, dihedral torsional profile prediction, and even bond breaking. We also demonstrate the smoothness, transferability, and extrapolation capabilities of ACE on a new challenging benchmark data set comprised of a potential energy surface of a flexible druglike molecule.

The Influence of Electronic Polarization on Nonlinear Optical Spectroscopy

The environment surrounding a chromophore can dramatically affect the energy absorption and relaxation process, as manifested in optical spectra. Simulations of nonlinear optical spectroscopy, such as two-dimensional electronic spectroscopy (2DES) and transient absorption (TA), will be influenced by the computational model of the environment. We here compare a fixed point charge molecular mechanics model and a quantum mechanical (QM) model of the environment in computed 2DES and TA spectra of Nile red in water and the chromophore of photoactive yellow protein (PYP) in water and protein environments. In addition to simulating these nonlinear optical spectra, we directly juxtapose the computed excitation energy correlation function to the dynamic Stokes shift function often used to analyze environment dynamics. Overall, we find that for the three systems studied here the mutual electronic polarization provided by the QM environment manifests in broader 2DES signals, as well as a larger reorganization energy and a larger static Stokes shift due to stronger coupling between the chromophore and the environment.

Topology Automated Force-Field Interactions (TAFFI): A Framework for Developing Transferable Force Fields

Force-field development has undergone a revolution in the past decade with the proliferation of quantum chemistry based parametrizations and the introduction of machine learning approximations of the atomistic potential energy surface. Nevertheless, transferable force fields with broad coverage of organic chemical space remain necessary for applications in materials and chemical discovery where throughput, consistency, and computational cost are paramount. Here, we introduce a force-field development framework called Topology Automated Force-Field Interactions (TAFFI) for developing transferable force fields of varying complexity against an extensible database of quantum chemistry calculations. TAFFI formalizes the concept of atom typing and makes it the basis for generating systematic training data that maintains a one-to-one correspondence with force-field terms. This feature makes TAFFI arbitrarily extensible to new chemistries while maintaining internal consistency and transferability. As a demonstration of TAFFI, we have developed a fixed-charge force-field, TAFFI-gen, from scratch that includes coverage for common organic functional groups that is comparable to established transferable force fields. The performance of TAFFI-gen was benchmarked against OPLS and GAFF for reproducing several experimental properties of 87 organic liquids. The consistent performance of these force fields, despite their distinct origins, validates the TAFFI framework while also providing evidence of the representability limitations of fixed-charge force fields.

Influence of non-adiabatic effects on linear absorption spectra in the condensed phase: Methylene blue

Modeling linear absorption spectra of solvated chromophores is highly challenging as contributions are present both from coupling of the electronic states to nuclear vibrations and from solute-solvent interactions. In systems where excited states intersect in the Condon region, significant non-adiabatic contributions to absorption line shapes can also be observed. Here, we introduce a robust approach to model linear absorption spectra accounting for both environmental and non-adiabatic effects from first principles. This model parameterizes a linear vibronic coupling (LVC) Hamiltonian directly from energy gap fluctuations calculated along molecular dynamics (MD) trajectories of the chromophore in solution, accounting for both anharmonicity in the potential and direct solute-solvent interactions. The resulting system dynamics described by the LVC Hamiltonian are solved exactly using the thermalized time-evolving density operator with orthogonal polynomials algorithm (T-TEDOPA). The approach is applied to the linear absorption spectrum of methylene blue in water. We show that the strong shoulder in the experimental spectrum is caused by vibrationally driven population transfer between the bright S₁ and the dark S₂ states. The treatment of the solvent environment is one of many factors that strongly influence the population transfer and line shape; accurate modeling can only be achieved through the use of explicit quantum mechanical solvation. The efficiency of T-TEDOPA, combined with LVC Hamiltonian parameterizations from MD, leads to an attractive method for describing a large variety of systems in complex environments from first principles.

Development and Benchmarking of Open Force Field v1.0.0-the Parsley Small-Molecule Force Field

We present a methodology for defining and optimizing a general force field for classical molecular simulations, and we describe its use to derive the Open Force Field 1.0.0 small-molecule force field, codenamed Parsley. Rather than using traditional atom typing, our approach is built on the SMIRKS-native Open Force Field (SMIRNOFF) parameter assignment formalism, which handles increases in the diversity and specificity of the force field definition without needlessly increasing the complexity of the specification. Parameters are optimized with the ForceBalance tool, based on reference quantum chemical data that include torsion potential energy profiles, optimized gas-phase structures, and vibrational frequencies. These quantum reference data are computed and are maintained with QCArchive, an open-source and freely available distributed computing and database software ecosystem. In this initial application of the method, we present essentially a full optimization of all valence parameters and report tests of the resulting force field against compounds and data types outside the training set. These tests show improvements in optimized geometries and conformational energetics and demonstrate that Parsley's accuracy for liquid properties is similar to that of other general force fields, as is accuracy on binding free energies. We find that this initial Parsley force field affords accuracy similar to that of other general force fields when used to calculate relative binding free energies spanning 199 protein-ligand systems. Additionally, the resulting infrastructure allows us to rapidly optimize an entirely new force field with minimal human intervention.

Systematic optimization of a fragment-based force field against experimental pure-liquid properties considering large compound families: application to oxygen and nitrogen compounds

The CombiFF approach is a workflow for the automated refinement of force-field parameters against experimental condensed-phase data, considering entire classes of organic molecules constructed using a fragment library via combinatorial isomer enumeration. One peculiarity of this approach is that it relies on an electronegativity-equalization scheme to account for induction effects within molecules, with values of the atomic hardness and electronegativity as electrostatic parameters, rather than the partial charges themselves. In a previous article [M. P. Oliveira, M. Andrey, S. R. Rieder, L. Kern, D. F. Hahn, S. Riniker, B. A. C. Horta and P. H. Hünenberger, J. Chem. Theory. Comput. 2020, 16, 7525], CombiFF was introduced and applied to calibrate a GROMOS-compatible united-atom force field for the saturated acyclic (halo-)alkane family. Here, this scheme is employed for the construction of a corresponding force field for saturated acyclic compounds encompassing eight common chemical functional groups involving oxygen and/or nitrogen atoms, namely: ether, aldehyde, ketone, ester, alcohol, carboxylic acid, amine, and amide. Monofunctional as well as homo-polyfunctional compounds are considered. A total of 1712 experimental liquid densities ρliq and vaporization enthalpies ΔHvap concerning 1175 molecules are used for the calibration (339 molecules) and validation (836 molecules) of the 102 non-bonded interaction parameters of the force field. Using initial parameter values based on the GROMOS 2016H66 parameter set, convergence is reached after five iterations. Given access to one processor per simulated system, this operation only requires a few days of wall-clock computing time. After optimization, the root-mean-square deviations from experiment are 29.9 (22.4) kg m-3 for ρliq and 4.1 (5.5) kJ mol-1 for ΔHvap for the calibration (validation) set. Thus, a very good level of agreement with experiment is achieved in terms of these two properties, although the errors are inhomogeneously distributed across the different chemical functional groups.

Modeling Molecular Emitters in Organic Light-Emitting Diodes with the Quantum Mechanical Bespoke Force Field

Combined molecular dynamics (MD) and quantum mechanics (QM) simulation procedures have gained popularity in modeling the spectral properties of functional organic molecules. However, the potential energy surfaces used to propagate long-time scale dynamics in these simulations are typically described using general, transferable force fields designed for organic molecules in their electronic ground states. These force fields do not typically include spectroscopic data in their training, and importantly, there is no general protocol for including changes in geometry or intermolecular interactions with the environment that may occur upon electronic excitation. In this work, we show that parameters tailored for thermally activated delayed fluorescence (TADF) emitters used in organic light-emitting diodes (OLEDs), in both their ground and electronically excited states, can be readily derived from a small number of QM calculations using the QUBEKit (QUantum mechanical BEspoke toolKit) software and improve the overall accuracy of these simulations.

Q-Force: Quantum Mechanically Augmented Molecular Force Fields

The quality of molecular dynamics simulations strongly depends on the accuracy of the underlying force fields (FFs) that determine all intra- and intermolecular interactions of the system. Commonly, transferable FF parameters are determined based on a representative set of small molecules. However, such an approach sacrifices accuracy in favor of generality. In this work, an open-source and automated toolkit named Q-Force is presented, which augments these transferable FFs with molecule-specific bonded parameters and atomic charges that are derived from quantum mechanical (QM) calculations. The molecular fragmentation procedure allows treatment of large molecules (>200 atoms) with a low computational cost. The generated Q-Force FFs can be used at the same computational cost as transferable FFs, but with improved accuracy: We demonstrate this for the vibrational properties on a set of small molecules and for the potential energy surface on a complex molecule (186 atoms) with photovoltaic applications. Overall, the accuracy, user-friendliness, and minimal computational overhead of the Q-Force protocol make it widely applicable for atomistic molecular dynamics simulations.

Automated Parameterization of Quantum Mechanically Derived Force Fields for Soft Materials and Complex Fluids: Development and Validation

The reliability of molecular dynamics (MD) simulations in predicting macroscopic properties of complex fluids and soft materials, such as liquid crystals, colloidal suspensions, or polymers, relies on the accuracy of the adopted force field (FF). We present an automated protocol to derive specific and accurate FFs, fully based on ab initio quantum mechanical (QM) data. The integration of the Joyce and Picky procedures, recently proposed by our group to provide an accurate description of simple liquids, is here extended to larger molecules, capable of exhibiting more complex fluid phases. While the standard Joyce protocol is employed to parameterize the intramolecular FF term, a new automated procedure is here proposed to handle the computational cost of the QM calculations required for the parameterization of the intermolecular FF term. The latter is thus obtained by integrating the old Picky procedure with a fragmentation reconstruction method (FRM) that allows for a reliable, yet computationally feasible sampling of the intermolecular energy surface at the QM level. The whole FF parameterization protocol is tested on a benchmark liquid crystal, and the performances of the resulting quantum mechanically derived (QMD) FF were compared with those delivered by a general-purpose, transferable one, and by the third, "hybrid" FF, where only the bonded terms were refined against QM data. Lengthy atomistic MD simulations are carried out with each FF on extended 5CB systems in both isotropic and nematic phases, eventually validating the proposed protocol by comparing the resulting macroscopic properties with other computational models and with experiments. The QMD-FF yields the best performances, reproducing both phases in the correct range of temperatures and well describing their structure, dynamics, and thermodynamic properties, thus providing a clear protocol that may be explored to predict such properties on other complex fluids or soft materials.

Optimized Halogen Atomic Radii for PBSA Calculations Using Off-Center Point Charges

In force-field methods, the usage of off-center point charges, also called extra points (EPs), is a common strategy to tackle the anisotropy of the electrostatic potential of covalently bonded halogens (X), thus allowing the description of halogen bonds (XBs) at the molecular mechanics/molecular dynamics (MM/MD) level. Diverse EP implementations exist in the literature differing on the charge sets and/or the X-EP distances. Poisson-Boltzmann and surface area (PBSA) calculations can be used to obtain solvation free energies (ΔG_solv) of small molecules, often to compute binding free energies (ΔG_bind) at the MM-PBSA level. This method depends, among other parameters, on the empirical assignment of atomic radii (PB radii). Given the multiplicity of off-center point-charge models and the lack of specific PB radii for halogens compatible with such implementations, in this work, we assessed the performance of PBSA calculations for the estimation of ΔG_solv values in water (ΔG_hyd), also conducting an optimization of the halogen PB radii (Cl, Br, and I) for each EP model. We not only expand the usage of EP models in the scope of the general AMBER force field (GAFF) but also provide the first optimized halogen PB radii in the context of the CHARMM general force field (CGenFF), thus contributing to improving the description of halogenated compounds in PBSA calculations.

Molecularly Informed Field Theories from Bottom-up Coarse-Graining

Polymer formulations possessing mesostructures or phase coexistence are challenging to simulate using atomistic particle-explicit approaches due to the disparate time and length scales, while the predictive capability of field-based simulations is hampered by the need to specify interactions at a coarser scale (e.g., χ-parameters). To overcome the weaknesses of both, we introduce a bottom-up coarse-graining methodology that leverages all-atom molecular dynamics to molecularly inform coarser field-theoretic models. Specifically, we use relative-entropy coarse-graining to parametrize particle models that are directly and analytically transformable into statistical field theories. We demonstrate the predictive capability of this approach by reproducing experimental aqueous poly(ethylene oxide) (PEO) cloud-point curves with no parameters fit to experimental data. This synergistic approach to multiscale polymer simulations opens the door to de novo exploration of phase behavior across a wide variety of polymer solutions and melt formulations.

Partial Charges Optimized by Genetic Algorithms for Deep Eutectic Solvent Simulations

Deep eutectic solvents (DESs) are a class of solvents often composed of ammonium-based chloride salts and a neutral hydrogen bond donor (HBD) at specific ratios. These cost-effective and environmentally friendly solvents have seen significant growth in multiple fields, including organic synthesis, and in materials and extractions because of their desirable properties. In the present work, a new software called genetic algorithm machine learning (GAML) was developed that utilizes a genetic algorithm (GA) approach to facilitate the development of optimized potentials for liquid simulation (OPLS)-based force field (FF) parameters for eight unique DESs based on three ammonium-based salts and five HBDs at multiple salt:HBD ratios. As an initial test of GAML, partial charges were created for 86 conventional solvents based on neutral organic molecules that yielded excellent overall mean absolute deviations (MADs) of 0.021 g/cm³, 0.63 kcal/mol, and 0.20 kcal/mol compared to experimental densities, heats of vaporization (ΔH_vap), and free energies of hydration (ΔG_hyd), respectively. FFs for DESs constructed from ethylammonium, N,N-diethylethanolammonium, and N-ethyl-N,N-dimethylethanolammonium chloride salts were then parameterized using GAML with exceptional agreement achieved at multiple temperatures for experimental densities, surface tensions, and viscosities with MADs of 0.024 g/cm³, 4.2 mN/m, and 5.3 cP, respectively.

Testing the Limitations of MD-Based Local Electric Fields Using the Vibrational Stark Effect in Solution: Penicillin G as a Test Case

Noncovalent interactions underlie nearly all molecular processes in the condensed phase from solvation to catalysis. Their quantification within a physically consistent framework remains challenging. Experimental vibrational Stark effect (VSE)-based solvatochromism can be combined with molecular dynamics (MD) simulations to quantify the electrostatic forces in solute-solvent interactions for small rigid molecules and, by extension, when these solutes bind in enzyme active sites. While generalizing this approach toward more complex (bio)molecules, such as the conformationally flexible and charged penicillin G (PenG), we were surprised to observe inconsistencies in MD-based electric fields. Combining synthesis, VSE spectroscopy, and computational methods, we provide an intimate view on the origins of these discrepancies. We observe that the electric fields are correlated to conformation-dependent effects of the flexible PenG side chain, including both the local solvation structure and solute conformational sampling in MD. Additionally, we identified that MD-based electric fields are consistently overestimated in three-point water models in the vicinity of charged groups; this cannot be entirely ameliorated using polarizable force fields (AMOEBA) or advanced water models. This work demonstrates the value of the VSE as a direct method for experiment-guided refinements of MD force fields and establishes a general reductionist approach to calibrating vibrational probes for complex (bio)molecules.

Implementation of the QUBE Force Field in SOMD for High-Throughput Alchemical Free-Energy Calculations

The quantum mechanical bespoke (QUBE) force-field approach has been developed to facilitate the automated derivation of potential energy function parameters for modeling protein-ligand binding. To date, the approach has been validated in the context of Monte Carlo simulations of protein-ligand complexes. We describe here the implementation of the QUBE force field in the alchemical free-energy calculation molecular dynamics simulation package SOMD. The implementation is validated by demonstrating the reproducibility of absolute hydration free energies computed with the QUBE force field across the SOMD and GROMACS software packages. We further demonstrate, by way of a case study involving two series of non-nucleoside inhibitors of HIV-1 reverse transcriptase, that the availability of QUBE in a modern simulation package that makes efficient use of graphics processing unit acceleration will facilitate high-throughput alchemical free-energy calculations.

Vibronic and Environmental Effects in Simulations of Optical Spectroscopy

Including both environmental and vibronic effects is important for accurate simulation of optical spectra, but combining these effects remains computationally challenging. We outline two approaches that consider both the explicit atomistic environment and the vibronic transitions. Both phenomena are responsible for spectral shapes in linear spectroscopy and the electronic evolution measured in nonlinear spectroscopy. The first approach utilizes snapshots of chromophore-environment configurations for which chromophore normal modes are determined. We outline various approximations for this static approach that assumes harmonic potentials and ignores dynamic system-environment coupling. The second approach obtains excitation energies for a series of time-correlated snapshots. This dynamic approach relies on the accurate truncation of the cumulant expansion but treats the dynamics of the chromophore and the environment on equal footing. Both approaches show significant potential for making strides toward more accurate optical spectroscopy simulations of complex condensed phase systems.

Semi-automated Optimization of the CHARMM36 Lipid Force Field to Include Explicit Treatment of Long-Range Dispersion

The development of the CHARMM lipid force field (FF) can be traced back to the early 1990s with its current version denoted CHARMM36 (C36). The parametrization of C36 utilized high-level quantum mechanical data and free energy calculations of model compounds before parameters were manually adjusted to yield agreement with experimental properties of lipid bilayers. While such manual fine-tuning of FF parameters is based on intuition and trial-and-error, automated methods can identify beneficial modifications of the parameters via their sensitivities and thereby guide the optimization process. This work introduces a semi-automated approach to reparametrize the CHARMM lipid FF with consistent inclusion of long-range dispersion through the Lennard-Jones particle-mesh Ewald (LJ-PME) approach. The optimization method is based on thermodynamic reweighting with regularization with respect to the C36 set. Two independent optimizations with different topology restrictions are presented. Targets of the optimizations are primarily liquid crystalline phase properties of lipid bilayers and the compression isotherm of monolayers. Pair correlation functions between water and lipid functional groups in aqueous solution are also included to address headgroup hydration. While the physics of the reweighting strategy itself is well-understood, applying it to heterogeneous, complex anisotropic systems poses additional challenges. These were overcome through careful selection of target properties and reweighting settings allowing for the successful incorporation of the explicit treatment of long-range dispersion, and we denote the newly optimized lipid force field as C36/LJ-PME. The current implementation of the optimization protocol will facilitate the future development of the CHARMM and related lipid force fields.

New atoms-in-molecules dispersion models for use in ab initio derived force fields

Recently, substantial research efforts have gone into bridging the accuracy-efficiency gap between parameterized force field models and quantum chemical calculations by extracting molecule-specific force fields directly from ab initio data in a robust and automated manner. One of the challenging aspects is deriving localized atomic polarizabilities for pairwise distributed dispersion models. The Tkatchenko-Scheffler model is based upon correcting free-atom C₆ coefficients according to the square of the ratio of the atom-in-molecule volume to the free-atom volume. However, it has recently been shown that a more accurate relationship can be found if static atomic polarizabilities are also taken into account. Using this relationship, we develop two modified Tkatchenko-Scheffler dispersion models and benchmark their performance against SAPT2+3 reference data and other commonly used dispersion models.

Automated parameterization of quantum-mechanically derived force-fields including explicit sigma holes: A pathway to energetic and structural features of halogen bonds in gas and condensed phase

In classical molecular dynamics, general purpose atomistic force-fields (FFs) often deliver inaccurate results when dealing with halogen bonds (XBs), notwithstanding their crucial role in many fields of science, ranging from material design to drug development. Given the large dimensions of the systems of interest, it would be therefore desirable to increase the FF accuracy maintaining the simplicity of the standard Lennard-Jones (LJ) plus point charge description to avoid an excessive computational cost. A simple yet effective strategy consists in introducing a number of virtual sites able to mimic the so-called "explicit σ-hole." In this work, we present an automated FF parameterization strategy based on a global optimization of both LJ and charge parameters with respect to accurate quantum mechanical data, purposely computed for the system under investigation. As a test case, we report on two homologue series, characterized either by weak or strong XBs, namely, the di-halogenated methanes and the mono-, di-, and tri-substituted acetonitriles, taking into consideration Cl, Br, and I substituents. The resulting quantum mechanically derived FFs are validated for each compound in the gas and in the condensed phase by comparing them to general purpose and specific FFs without virtual sites and to highly accurate reference quantum mechanical data. The results strongly support the adoption of the specific FFs with virtual sites, which overcome the other investigated models in representing both gas phase energetics and the structural patterns of the liquid phase structure related to the presence of XBs.

Charge Anisotropy of Nitrogen: Where Chemical Intuition Fails

For more than half a century computer simulations were developed and employed to study ensemble properties of a wide variety of atomic and molecular systems with tremendous success. Nowadays, a selection of force-fields is available that describe the interactions in such systems. A key feature of force-fields is an adequate description of the electrostatic potential (ESP). Several force-fields model the ESP via point charges positioned at the atom centers. A major shortcoming of this approach, its inability to model anisotropies in the ESP, can be mitigated using additional charge sites. It has been shown that nitrogen is the most problematic element abundant in many polymers as well as large molecules of biological origin. To tackle this issue, small organic molecules containing a single nitrogen atom were studied. In performing rigorous scans of the surroundings of these nitrogen atoms, positions where a single extra charge can enhance the ESP description the most were identified. Significant improvements are found for ammonia, amines, and amides. Interestingly, the optimal location for the extra charge does not correlate with the chemically intuitive position of the nitrogen lone pair. In fact, the placement of an extra charge in the lone-pair location does not lead to significant improvements in most cases.

Iron's Wake: The Performance of Quantum Mechanical-Derived Versus General-Purpose Force Fields Tested on a Luminescent Iron Complex

Recently synthetized iron complexes have achieved long-lived excited states and stabilities which are comparable, or even superior, to their ruthenium analogues, thus representing an eco-friendly and cheaper alternative to those materials based on rare metals. Most of computational tools which could help unravel the origin of this large efficiency rely on ab-initio methods which are not able, however, to capture the nanosecond time scale underlying these photophysical processes and the influence of their realistic environment. Therefore, it exists an urgent need of developing new low-cost, but still accurate enough, computational methodologies capable to deal with the steady-state and transient spectroscopy of transition metal complexes in solution. Following this idea, here we focus on the comparison between general-purpose transferable force-fields (FFs), directly available from existing databases, and specific quantum mechanical derived FFs (QMD-FFs), obtained in this work through the Joyce procedure. We have chosen a recently reported Fe^III complex with nanosecond excited-state lifetime as a representative case. Our molecular dynamics (MD) simulations demonstrated that the QMD-FF nicely reproduces the structure and the dynamics of the complex and its chemical environment within the same precision as higher cost QM methods, whereas general-purpose FFs failed in this purpose. Although in this particular case the chemical environment plays a minor role on the photo physics of this system, these results highlight the potential of QMD-FFs to rationalize photophysical phenomena provided an accurate QM method to derive its parameters is chosen.

The ONETEP linear-scaling density functional theory program

We present an overview of the onetep program for linear-scaling density functional theory (DFT) calculations with large basis set (plane-wave) accuracy on parallel computers. The DFT energy is computed from the density matrix, which is constructed from spatially localized orbitals we call Non-orthogonal Generalized Wannier Functions (NGWFs), expressed in terms of periodic sinc (psinc) functions. During the calculation, both the density matrix and the NGWFs are optimized with localization constraints. By taking advantage of localization, onetep is able to perform calculations including thousands of atoms with computational effort, which scales linearly with the number or atoms. The code has a large and diverse range of capabilities, explored in this paper, including different boundary conditions, various exchange-correlation functionals (with and without exact exchange), finite electronic temperature methods for metallic systems, methods for strongly correlated systems, molecular dynamics, vibrational calculations, time-dependent DFT, electronic transport, core loss spectroscopy, implicit solvation, quantum mechanical (QM)/molecular mechanical and QM-in-QM embedding, density of states calculations, distributed multipole analysis, and methods for partitioning charges and interactions between fragments. Calculations with onetep provide unique insights into large and complex systems that require an accurate atomic-level description, ranging from biomolecular to chemical, to materials, and to physical problems, as we show with a small selection of illustrative examples. onetep has always aimed to be at the cutting edge of method and software developments, and it serves as a platform for developing new methods of electronic structure simulation. We therefore conclude by describing some of the challenges and directions for its future developments and applications.

Deriving a Polarizable Force Field for Biomolecular Building Blocks with Minimal Empirical Calibration

Force field parametrization involves a complex set of linked optimization problems, with the goal of describing complex molecular interactions by using simple classical potential-energy functions that model Coulomb interactions, dispersion, and exchange repulsion. These functions comprise a set of atomic (and molecular) parameters and together with the bonded terms they constitute the molecular mechanics force field. Traditionally, many of these parameters have been fitted in a calibration approach in which experimental measures for thermodynamic and other relevant properties of small-molecule compounds are used for fitting and validation. As these approaches are laborious and time-consuming and represent an underdetermined optimization problem, we study methods to fit and derive an increasing number of parameters directly from electronic structure calculations, in order to greatly reduce possible parameter space for the remaining free parameters. In the current work we investigate a polarizable model with a higher order dispersion term for use in biomolecular simulation. Results for 49 biochemically relevant molecules are presented including updated parameters for hydrocarbon side chains. We show that our recently presented set of QM/MM derived atomic polarizabilities can be used in direct conjunction with partial charges and a higher order dispersion model that are quantum-mechanically determined, to freeze nearly all (i.e., 132 out of 138) nonbonded parameters to their quantum determined values.

Data-Driven Mapping of Gas-Phase Quantum Calculations to General Force Field Lennard-Jones Parameters

Molecular dynamics simulations are helpful tools for a range of applications, ranging from drug discovery to protein structure determination. The successful use of this technology largely depends on the potential function, or force field, used to determine the potential energy at each configuration of the system. Most force fields encode all of the relevant parameters to be used in distinct atom types, each associated with parameters for all parts of the force field, typically bond stretches, angle bends, torsions, and nonbonded terms accounting for van der Waals and electrostatic interactions. Much attention has been paid to the nonbonded parameters and their derivation, which are important in particular due to their governance of noncovalent interactions, such as protein-ligand binding. Parametrization involves adjusting the nonbonded parameters to minimize the error between simulation results and experimental properties, such as heats of vaporization and densities of neat liquids. In this setting, determining the best set of atom types is far from trivial, and the large number of parameters to be fit for the atom types in a typical force field can make it difficult to approach a true optimum. Here, we utilize a previously described Minimal Basis Iterative Stockholder (MBIS) method to carry out an atoms-in-molecules partitioning of electron densities. Information from these atomic densities is then mapped to Lennard-Jones parameters using a set of mapping parameters much smaller than the typical number of atom types in a force field. This approach is advantageous in two ways: it eliminates atom types by allowing each atom to have unique Lennard-Jones parameters, and it greatly reduces the number of parameters to be optimized. We show that this approach yields results comparable to those obtained with the typed GAFF 1.7 force field, even when trained on a relatively small amount of experimental data.

A machine learning based intramolecular potential for a flexible organic molecule

Here, we employ the kernel regression machine learning technique to construct an analytical potential that reproduces the quantum mechanical potential energy surface of a small, flexible, drug-like molecule, 3-(benzyloxy)pyridin-2-amine.

Identifying misbonded atoms in the 2019 CoRE metal–organic framework database

Databases of experimentally-derived metal–organic framework (MOF) crystal structures are useful for large-scale computational screening to identify which MOFs are best-suited for particular applications. However, these crystal structures must be cleaned to identify and/or correct various artifacts. The recently published 2019 CoRE MOF database (Chung et al., J. Chem. Eng. Data, 2019, 64, 5985–5998) reported thousands of experimentally-derived crystal structures that were partially cleaned to remove solvent molecules, to identify hundreds of disordered structures (approximately thirty of those were corrected), and to manually correct approximately 100 structures (e.g., adding missing hydrogen atoms). Herein, further cleaning of the 2019 CoRE MOF database is performed to identify structures with misbonded or isolated atoms: (i) structures containing an isolated atom, (ii) structures containing atoms too close together (i.e., overlapping atoms), (iii) structures containing a misplaced hydrogen atom, (iv) structures containing an under-bonded carbon atom (which might be caused by missing hydrogen atoms), and (v) structures containing an over-bonded carbon atom. This study should not be viewed as the final cleaning of this database, but rather as progress along the way towards the goal of someday achieving a completely cleaned set of experimentally-derived MOF crystal structures. We performed atom typing for all of the accepted structures to identify those structures that can be parameterized by previously reported forcefield precursors (Chen and Manz, RSC Adv., 2019, 9, 36492–36507). We report several forcefield precursors (e.g., net atomic charges, atom-in-material polarizabilities, atom-in-material dispersion coefficients, electron cloud parameters, etc.) for more than five thousand MOFs in the 2019 CoRE MOF database.

The 2019 CoRE MOF database was cleaned by identifying structures containing isolated atoms, overlapping atoms, misplaced hydrogens, and under/over-bonded carbons.

Modelling flexible protein–ligand binding in p38α MAP kinase using the QUBE force field

The accuracy of quantum mechanical bespoke (QUBE) force fields for protein–ligand binding free energy calculations are benchmarked against experiment.

Torsional Energy Barriers of Biaryls Could Be Predicted by Electron Richness/Deficiency of Aromatic Rings; Advancement of Molecular Mechanics toward Atom-Type Independence

Biaryl molecules are ubiquitous pharmacophores found in natural products and pharmaceuticals. In spite of this, existing molecular mechanics force fields are unable to accurately reproduce their torsional energy profiles, except for a few well-parametrized cases. This effectively limits the ability of structure-based drug design methods to correctly identify hits involving biaryls with confidence (e.g., during virtual screening, employing docking and/or molecular dynamics simulations). Continuing in our endeavor to quantify organic chemistry principles, we showed that the torsional energy profile of biaryl compounds could be computed on-the-fly based on the electron richness/deficiency of the aromatic rings. This method, called H-TEQ 4.0, was developed using a set of 131 biaryls. It was subsequently validated on a separate set of 100 diverse biaryls, including multisubstituted, bicyclic and tricyclic druglike molecules, and produced an average root-mean-square error (RMSE) of 0.95 kcal·mol^-1. For comparison, GAFF2 produced an RMSE of 3.88 kcal·mol^-1, owing to problems associated with the transferability of torsion parameters. The success of H-TEQ 4.0 provided further evidence that force fields could transition to become atom-type independent, providing that the correct chemical principles are used. Overall, this method solved the problem of transferability of biaryl torsion parameters, while simultaneously improving the overall accuracy of the force field.

A collection of forcefield precursors for metal-organic frameworks

A host of important performance properties for metal-organic frameworks (MOFs) and other complex materials can be calculated by modeling statistical ensembles. The principle challenge is to develop accurate and computationally efficient interaction models for these simulations. Two major approaches are (i) ab initio molecular dynamics in which the interaction model is provided by an exchange-correlation theory (e.g., DFT + dispersion functional) and (ii) molecular mechanics in which the interaction model is a parameterized classical force field. The first approach requires further development to improve computational speed. The second approach requires further development to automate accurate forcefield parameterization. Because of the extreme chemical diversity across thousands of MOF structures, this problem is still mostly unsolved today. For example, here we show structures in the 2014 CoRE MOF database contain more than 8 thousand different atom types based on first and second neighbors. Our results showed that atom types based on both first and second neighbors adequately capture the chemical environment, but atom types based on only first neighbors do not. For 3056 MOFs, we used density functional theory (DFT) followed by DDEC6 atomic population analysis to extract a host of important forcefield precursors: partial atomic charges; atom-in-material (AIM) C6, C8, and C10 dispersion coefficients; AIM dipole and quadrupole moments; various AIM polarizabilities; quantum Drude oscillator parameters; AIM electron cloud parameters; etc. Electrostatic parameters were validated through comparisons to the DFT-computed electrostatic potential. These forcefield precursors should find widespread applications to developing MOF force fields.