Alma Carolina Escobosa

Correction for ‘An automated protocol to construct flexibility parameters for classical forcefields: applications to metal–organic frameworks’ by Reza Ghanavati et al., RSC Adv., 2024, 14, 22714–22762, https://doi.org/10.1039/D4RA01859A.

Atom-centered point-charge models are computationally efficient and commonly used electrostatic models to build forcefields used in classical simulations of materials. To assess their performance, we evaluated various atomic charge assignment methods across the following material types: (a) organic molecules, (b) inorganic molecules, (c) heterodiatomic molecules, (d) transition metal complexes, and (e) nanoporous crystals. We compared 12 atomic charge assignment methods for molecular systems and 6 for nanoporous crystals. In this article, we introduce a computationally efficient quadrupole-dipole-resorption (QDR) method that improves the accuracy of stockholder-partitioning models (e.g., DDEC6) for approximately reproducing the electrostatic potential and molecular dipole and quadrupole moments. For each electrostatic model, we computed the electrostatic potential's root-mean-squared error (RMSE) and relative RMSE (RRMSE) using the material's QM-computed electrostatic potential as a reference. The electrostatic RMSE and RRMSE were computed for a training dataset containing 21 geometric conformations per material and a validation dataset containing 20 new geometric conformations per material. Raincloud plots were prepared to visualize the resulting data distributions. For each charge assignment method in the nonperiodic materials, we also computed and compared the root-mean-squared charge transfer magnitude, correlations to other charge assignment methods, conformational sensitivity, etc. Among point-charge models, multiframe ESP methods provided the best overall accuracy for reproducing the electrostatic potential across different system conformations, but they require a training set containing many geometric conformations. The QDR-DDEC6 and CM5 methods provided good conformational transferability and electrostatic potential accuracy across various conformations even when trained only on a single optimized ground-state geometry. A Pareto plot was prepared illustrating the tradeoff between conformational sensitivity and accuracy for reproducing the electrostatic potential of individual conformations. Including atomic dipoles (e.g., QDR-DDEC6_ad, DDEC6_ad, and MBIS_ad) greatly improved the electrostatic model for individual conformations, and QDR-DDEC6_ad outperformed all atom-centered point-charge models. We recommend that more computationally efficient methods be developed to use electrostatic models containing atom-centered point charges plus atomic dipoles in flexible forcefields. Finally, some electron-density partitioning approaches have the key advantage of providing accurate results even when applied to dense solids under high pressures, and we demonstrated this using 11 sodium chloride crystals having various stoichiometries.

In this work, forcefield flexibility parameters were constructed and validated for more than 100 metal–organic frameworks (MOFs). We used atom typing to identify bond types, angle types, and dihedral types associated with bond stretches, angle bends, dihedral torsions, and other flexibility interactions. Our work used Manz's angle-bending and dihedral-torsion model potentials. For a crystal structure containing N_atoms in its unit cell, the number of independent flexibility interactions is 3(N_atoms − 1). Because the number of bonds, angles, and dihedrals is normally much larger than 3(N_atoms − 1), these internal coordinates are redundant. To reduce (but not eliminate) this redundancy, our protocol prunes dihedral types in a way that preserves symmetry equivalency. Next, each dihedral type is classified as non-rotatable, hindered, rotatable, or linear. We introduce a smart selection method that identifies which particular torsion modes are important for each rotatable dihedral type. Then, we computed the force constants for all flexibility interactions together via LASSO regression (i.e., regularized linear least-squares fitting) of the training dataset. LASSO automatically identifies and removes unimportant forcefield interactions. For each MOF, the reference dataset was quantum-mechanically-computed in VASP via DFT with dispersion and included: (i) finite-displacement calculations along every independent atom translation mode, (ii) geometries randomly sampled via ab initio molecular dynamics (AIMD), (iii) the optimized ground-state geometry using experimental lattice parameters, and (iv) rigid torsion scans for each rotatable dihedral type. After training, the flexibility model was validated across geometries that were not part of the training dataset. For each MOF, we computed the goodness of fit (R-squared value) and the root-mean-squared error (RMSE) separately for the training and validation datasets. We compared flexibility models with and without bond–bond cross terms. Even without cross terms, the model yielded R-squared values of 0.910 (avg across all MOFs) ± 0.018 (st. dev.) for atom-in-material forces in the validation datasets. Our SAVESTEPS protocol should find widespread applications to parameterize flexible forcefields for material datasets. We performed molecular dynamics simulations using these flexibility parameters to compute heat capacities and thermal expansion coefficients for two MOFs.