Eleven NanoHUB Simulation Tools Using RASPA Software To Demonstrate Classical Atomistic Simulations of Fluids and Nanoporous Materials.

Eleven interactive simulation tools were created on nanoHUB to help users learn how to perform classical atomistic simulations. These tools enable users to perform classical Monte Carlo and molecular dynamics simulations using RASPA software. These tools use comparatively small numbers of production cycles to keep the runtimes short, so that users will not be discouraged by long wait times to see results. Here, we show that these tools produce results of sufficient accuracy and reproducibility for learning purposes. The 11 tools developed were as follows: (1) calculation of the self-diffusion constant of gas molecules in metal-organic frameworks (MOFs), (2) gas adsorption in MOFs using the grand canonical ensemble, (3) Henry's coefficient calculator for gas molecules in MOFs and a zeolite, (4) adsorption of a gas mixture in a MOF, (5) self-diffusion of a gas mixture in a MOF, (6) void fraction calculation for several MOFs and zeolites, (7) surface area calculation for several MOFs and zeolites, (8) calculation of radial distribution function and self-diffusion constant for several pure gases, (9) energy distribution of adsorption sites using a probe molecule in MOFs, (10) molecular dynamics simulation of pure fluids in the NPT ensemble, and (11) gas adsorption in MOFs using the Gibbs ensemble.

Apples to apples comparison of standardized to unstandardized principal component analysis of methods that assign partial atomic charges in molecules.

Articles by Cho et al. (ChemPhysChem, 2020, 21, 688-696) and Manz (RSC Adv., 2020, 10, 44121-44148) performed unstandardized and standardized, respectively, principal component analysis (PCA) to study atomic charge assignment methods for molecular systems. Both articles used subsets of atomic charges computed by Cho et al.; however, the data subsets employed were not strictly identical. Herein, an element by element analysis of this dataset is first performed to compare the spread of charge values across individual chemical elements and charge assignment methods. This reveals an underlying problem with the reported Becke partial atomic charges in this dataset. Due to their unphysical values, these Becke charges were not included in the subsequent PCA. Standardized and unstandardized PCA are performed across two datasets: (i) 19 charge assignment methods having a complete basis set limit and (ii) all 25 charge assignment methods (excluding Becke) for which Cho et al. computed atomic charges. The dataset contained ∼2000 molecules having a total of 29 907 atoms in materials. The following five methods (listed here in alphabetical order) showed the greatest correlation to the first principal component in standardized and unstandardized PCA: DDEC6, Hirshfeld-I, ISA, MBIS, and MBSBickelhaupt (note: MBSBickelhaupt does not appear in the 19 methods dataset). For standardized PCA, the DDEC6 method ranked first followed closely by MBIS. For unstandardized PCA, Hirshfeld-I (19 methods) or MBSBickelhaupt (25 methods) ranked first followed by DDEC6 in second place (both 19 and 25 methods).

Correction: Introducing DDEC6 atomic population analysis: part 1. Charge partitioning theory and methodology.

[This corrects the article DOI: 10.1039/C6RA04656H.].

Parametrization of Nonbonded Force Field Terms for Metal-Organic Frameworks Using Machine Learning Approach.

The enormous structural and chemical diversity of metal-organic frameworks (MOFs) forces researchers to actively use simulation techniques as often as experiments. MOFs are widely known for their outstanding adsorption properties, so a precise description of the host-guest interactions is essential for high-throughput screening aimed at ranking the most promising candidates. However, highly accurate ab initio calculations cannot be routinely applied to model thousands of structures due to the demanding computational costs. Furthermore, methods based on force field (FF) parametrization suffer from low transferability. To resolve this accuracy-efficiency dilemma, we applied a machine learning (ML) approach: extreme gradient boosting. The trained models reproduced the atom-in-material quantities, including partial charges, polarizabilities, dispersion coefficients, quantum Drude oscillator, and electron cloud parameters, with accuracy similar to the reference data set. The aforementioned FF precursors make it possible to thoroughly describe noncovalent interactions typical for MOF-adsorbate systems: electrostatic, dispersion, polarization, and short-range repulsion. The presented approach can also readily facilitate hybrid atomistic simulation/ML workflows.

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.

Seven confluence principles: a case study of standardized statistical analysis for 26 methods that assign net atomic charges in molecules.

This article studies two kinds of information extracted from statistical correlations between methods for assigning net atomic charges (NACs) in molecules. First, relative charge transfer magnitudes are quantified by performing instant least squares fitting (ILSF) on the NACs reported by Cho et al. (ChemPhysChem, 2020, 21, 688-696) across 26 methods applied to ∼2000 molecules. The Hirshfeld and Voronoi deformation density (VDD) methods had the smallest charge transfer magnitudes, while the quantum theory of atoms in molecules (QTAIM) method had the largest charge transfer magnitude. Methods optimized to reproduce the molecular dipole moment (e.g., ACP, ADCH, CM5) have smaller charge transfer magnitudes than methods optimized to reproduce the molecular electrostatic potential (e.g., CHELPG, HLY, MK, RESP). Several methods had charge transfer magnitudes even larger than the electrostatic potential fitting group. Second, confluence between different charge assignment methods is quantified to identify which charge assignment method produces the best NAC values for predicting via linear correlations the results of 20 charge assignment methods having a complete basis set limit across the dataset of ∼2000 molecules. The DDEC6 NACs were the best such predictor of the entire dataset. Seven confluence principles are introduced explaining why confluent quantitative descriptors offer predictive advantages for modeling a broad range of physical properties and target applications. These confluence principles can be applied in various fields of scientific inquiry. A theory is derived showing confluence is better revealed by standardized statistical analysis (e.g., principal components analysis of the correlation matrix and standardized reversible linear regression) than by unstandardized statistical analysis. These confluence principles were used together with other key principles and the scientific method to make assigning atom-in-material properties non-arbitrary. The N@C60 system provides an unambiguous and non-arbitrary falsifiable test of atomic population analysis methods. The HLY, ISA, MK, and RESP methods failed for this material.

New scaling relations to compute atom-in-material polarizabilities and dispersion coefficients: part 1. Theory and accuracy.

Polarizabilities and London dispersion forces are important to many chemical processes. Force fields for classical atomistic simulations can be constructed using atom-in-material polarizabilities and C n (n = 6, 8, 9, 10…) dispersion coefficients. This article addresses the key question of how to efficiently assign these parameters to constituent atoms in a material so that properties of the whole material are better reproduced. We develop a new set of scaling laws and computational algorithms (called MCLF) to do this in an accurate and computationally efficient manner across diverse material types. We introduce a conduction limit upper bound and m-scaling to describe the different behaviors of surface and buried atoms. We validate MCLF by comparing results to high-level benchmarks for isolated neutral and charged atoms, diverse diatomic molecules, various polyatomic molecules (e.g., polyacenes, fullerenes, and small organic and inorganic molecules), and dense solids (including metallic, covalent, and ionic). We also present results for the HIV reverse transcriptase enzyme complexed with an inhibitor molecule. MCLF provides the non-directionally screened polarizabilities required to construct force fields, the directionally-screened static polarizability tensor components and eigenvalues, and environmentally screened C6 coefficients. Overall, MCLF has improved accuracy compared to the TS-SCS method. For TS-SCS, we compared charge partitioning methods and show DDEC6 partitioning yields more accurate results than Hirshfeld partitioning. MCLF also gives approximations for C8, C9, and C10 dispersion coefficients and quantum Drude oscillator parameters. This method should find widespread applications to parameterize classical force fields and density functional theory (DFT) + dispersion methods.

Bond orders of the diatomic molecules.

Bond order quantifies the number of electrons dressed-exchanged between two atoms in a material and is important for understanding many chemical properties. Diatomic molecules are the smallest molecules possessing chemical bonds and play key roles in atmospheric chemistry, biochemistry, lab chemistry, and chemical manufacturing. Here we quantum-mechanically calculate bond orders for 288 diatomic molecules and ions. For homodiatomics, we show bond orders correlate to bond energies for elements within the same chemical group. We quantify and discuss how semicore electrons weaken bond orders for elements having diffuse semicore electrons. Lots of chemistry is effected by this. We introduce a first-principles method to represent orbital-independent bond order as a sum of orbital-dependent bond order components. This bond order component analysis (BOCA) applies to any spin-orbitals that are unitary transformations of the natural spin-orbitals, with or without periodic boundary conditions, and to non-magnetic and (collinear or non-collinear) magnetic materials. We use this BOCA to study all period 2 homodiatomics plus Mo2, Cr2, ClO, ClO-, and Mo2(acetate)4. Using Manz's bond order equation with DDEC6 partitioning, the Mo-Mo bond order was 4.12 in Mo2 and 1.46 in Mo2(acetate)4 with a sum of bond orders for each Mo atom of ∼4. Our study informs both chemistry research and education. As a learning aid, we introduce an analogy between bond orders in materials and message transmission in computer networks. We also introduce the first working quantitative heuristic model for all period 2 homodiatomic bond orders. This heuristic model incorporates s-p mixing to give heuristic bond orders of ¾ (Be2), 1¾ (B2), 2¾ (C2), and whole number bond orders for the remaining period 2 homodiatomics.

New scaling relations to compute atom-in-material polarizabilities and dispersion coefficients: part 2. Linear-scaling computational algorithms and parallelization.

We present two algorithms to compute system-specific polarizabilities and dispersion coefficients such that required memory and computational time scale linearly with increasing number of atoms in the unit cell for large systems. The first algorithm computes the atom-in-material (AIM) static polarizability tensors, force-field polarizabilities, and C 6, C 8, C 9, C 10 dispersion coefficients using the MCLF method. The second algorithm computes the AIM polarizability tensors and C 6 coefficients using the TS-SCS method. Linear-scaling computational cost is achieved using a dipole interaction cutoff length function combined with iterative methods that avoid large dense matrix multiplies and large matrix inversions. For MCLF, Richardson extrapolation of the screening increments is used. For TS-SCS, a failproof conjugate residual (FCR) algorithm is introduced that solves any linear equation system having Hermitian coefficients matrix. These algorithms have mathematically provable stable convergence that resists round-off errors. We parallelized these methods to provide rapid computation on multi-core computers. Excellent parallelization efficiencies were obtained, and adding parallel processors does not significantly increase memory requirements. This enables system-specific polarizabilities and dispersion coefficients to be readily computed for materials containing millions of atoms in the unit cell. The largest example studied herein is an ice crystal containing >2 million atoms in the unit cell. For this material, the FCR algorithm solved a linear equation system containing >6 million rows, 7.57 billion interacting atom pairs, 45.4 billion stored non-negligible matrix components used in each large matrix-vector multiplication, and ∼19 million unknowns per frequency point (>300 million total unknowns).

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.

Synthesis and Characterization of N,N'-Bismesityl Phenanthrene-9,10-diimine and Imine-Nitrone.

The sterically bulky compounds N,N'-bismesityl phenanthrene-9,10-diimine [1] and imine-nitrone [2] were synthesized. To the best of our knowledge, this is the first report of the synthesis of a bulky steric imine-nitrone accessed from the secondary ketimine using urea hydrogen peroxide over methyltrioxorhenium catalyst. Purified compounds were characterized using 1H and 13C NMR, high-resolution mass spectrometry, and infrared spectrometry. We report the first crystal structure of compound 1. Detailed IR bands of compounds 1 and 2 were assigned by comparing experimentally measured spectra to individually animated modes of quantum mechanically computed spectra. We believe these compounds may be of use as bidentate ligands in the synthesis of novel organometallic compounds. The asymmetric N and O coordination sites of compound 2 might impart interesting electronic effects to organometallic compounds compared to the symmetric N,N'-coordination sites of compound 1.

Introducing DDEC6 atomic population analysis: part 4. Efficient parallel computation of net atomic charges, atomic spin moments, bond orders, and more.

The DDEC6 method is one of the most accurate and broadly applicable atomic population analysis methods. It works for a broad range of periodic and non-periodic materials with no magnetism, collinear magnetism, and non-collinear magnetism irrespective of the basis set type. First, we show DDEC6 charge partitioning to assign net atomic charges corresponds to solving a series of 14 Lagrangians in order. Then, we provide flow diagrams for overall DDEC6 analysis, spin partitioning, and bond order calculations. We wrote an OpenMP parallelized Fortran code to provide efficient computations. We show that by storing large arrays as shared variables in cache line friendly order, memory requirements are independent of the number of parallel computing cores and false sharing is minimized. We show that both total memory required and the computational time scale linearly with increasing numbers of atoms in the unit cell. Using the presently chosen uniform grids, computational times of ∼9 to 94 seconds per atom were required to perform DDEC6 analysis on a single computing core in an Intel Xeon E5 multi-processor unit. Parallelization efficiencies were usually >50% for computations performed on 2 to 16 cores of a cache coherent node. As examples we study a B-DNA decamer, nickel metal, supercells of hexagonal ice crystals, six X@C60 endohedral fullerene complexes, a water dimer, a Mn12-acetate single molecule magnet exhibiting collinear magnetism, a Fe4O12N4C40H52 single molecule magnet exhibiting non-collinear magnetism, and several spin states of an ozone molecule. Efficient parallel computation was achieved for systems containing as few as one and as many as >8000 atoms in a unit cell. We varied many calculation factors (e.g., grid spacing, code design, thread arrangement, etc.) and report their effects on calculation speed and precision. We make recommendations for excellent performance.

Expanding the Scope of Density Derived Electrostatic and Chemical Charge Partitioning to Thousands of Atoms.

The density derived electrostatic and chemical (DDEC/c3) method is implemented into the onetep program to compute net atomic charges (NACs), as well as higher-order atomic multipole moments, of molecules, dense solids, nanoclusters, liquids, and biomolecules using linear-scaling density functional theory (DFT) in a distributed memory parallel computing environment. For a >1000 atom model of the oxygenated myoglobin protein, the DDEC/c3 net charge of the adsorbed oxygen molecule is approximately -1e (in agreement with the Weiss model) using a dynamical mean field theory treatment of the iron atom, but much smaller in magnitude when using the generalized gradient approximation. For GaAs semiconducting nanorods, the system dipole moment using the DDEC/c3 NACs is about 5% higher in magnitude than the dipole computed directly from the quantum mechanical electron density distribution, and the DDEC/c3 NACs reproduce the electrostatic potential to within approximately 0.1 V on the nanorod's solvent-accessible surface. As examples of conducting materials, we study (i) a 55-atom Pt cluster with an adsorbed CO molecule and (ii) the dense solids Mo2C and Pd3V. Our results for solid Mo2C and Pd3V confirm the necessity of a constraint enforcing exponentially decaying electron density in the tails of buried atoms.

Comment on "Extending Hirshfeld-I to bulk and periodic materials".

In recent years, several methods have been developed that partition the electron density among atoms using spherically symmetric atomic weights. D. E. P. Vanpoucke, P. Bultinck, and I. Van Driessche (J. Comput. Chem. 2012, doi: 10.1002/jcc.23088) recently reported a periodic implementation of the Hirshfeld-I method that uses a combination of Becke-style and uniform integration grids and modified atomic reference densities to compute net atomic charges in periodic materials. Herein, this method is discussed in the context of earlier periodic implementations of the Hirshfeld-I method, the Iterated Stockholder Atoms method, and the density derived electrostatic and chemical method.

Improved Atoms-in-Molecule Charge Partitioning Functional for Simultaneously Reproducing the Electrostatic Potential and Chemical States in Periodic and Nonperiodic Materials.

We develop a nonempirical atoms-in-molecules (AIM) method for computing net atomic charges that simultaneously reproduce chemical states of atoms in a material and the electrostatic potential V(r) outside its electron distribution. This method gives accurate results for a variety of periodic and nonperiodic materials including molecular systems, solid surfaces, porous solids, and nonporous solids. This method, called DDEC/c3, improves upon our previously published DDEC/c2 method (Manz, T. A.; Sholl, D. S. J. Chem. Theory Comput. 2010, 6, 2455-2468) by accurately treating nonporous solids with short bond lengths. Starting with the theory all AIM charge partitioning functionals with spherically symmetric atomic weights must satisfy, the form of the DDEC/c3 functional is derived from first principles. The method is designed to converge robustly by avoiding conditions that lead to nearly flat optimization landscapes. In addition to net atomic charges, the method can also compute atomic multipoles and atomic spin moments. Calculations performed on a variety of systems demonstrate the method's accuracy, computational efficiency, and good agreement with available experimental data. Comparisons to a variety of other charge assignment methods (Bader, natural population analysis, electrostatic potential fitting, Hirshfeld, iterative Hirshfeld, and iterative stockholder atoms) show that the DDEC/c3 net atomic charges are well-suited for constructing flexible force-fields for atomistic simulations.

Methods for Computing Accurate Atomic Spin Moments for Collinear and Noncollinear Magnetism in Periodic and Nonperiodic Materials.

The partitioning of electron spin density among atoms in a material gives atomic spin moments (ASMs), which are important for understanding magnetic properties. We compare ASMs computed using different population analysis methods and introduce a method for computing density derived electrostatic and chemical (DDEC) ASMs. Bader and DDEC ASMs can be computed for periodic and nonperiodic materials with either collinear or noncollinear magnetism, while natural population analysis (NPA) ASMs can be computed for nonperiodic materials with collinear magnetism. Our results show Bader, DDEC, and (where applicable) NPA methods give similar ASMs, but different net atomic charges. Because they are optimized to reproduce both the magnetic field and the chemical states of atoms in a material, DDEC ASMs are especially suitable for constructing interaction potentials for atomistic simulations. We describe the computation of accurate ASMs for (a) a variety of systems using collinear and noncollinear spin DFT, (b) highly correlated materials (e.g., magnetite) using DFT+U, and (c) various spin states of ozone using coupled cluster expansions. The computed ASMs are in good agreement with available experimental results for a variety of periodic and nonperiodic materials. Examples considered include the antiferromagnetic metal organic framework Cu3(BTC)2, several ozone spin states, mono- and binuclear transition metal complexes, ferri- and ferro-magnetic solids (e.g., Fe3O4, Fe3Si), and simple molecular systems. We briefly discuss the theory of exchange-correlation functionals for studying noncollinear magnetism. A method for finding the ground state of systems with highly noncollinear magnetism is introduced. We use these methods to study the spin-orbit coupling potential energy surface of the single molecule magnet Fe4C40H52N4O12, which has highly noncollinear magnetism, and find that it contains unusual features that give a new interpretation to experimental data.

A dimensionless reaction coordinate for quantifying the lateness of transition states.

The Hammond-Leffler postulate asserts that transition states of exothermic reactions are reactant-like (early), whereas transition states of endothermic reactions are product-like (late). Related postulates have been proposed to describe the sensitivity of activation barriers for reactions occurring on catalytic surfaces to the catalyst structure. To evaluate the validity of these postulates for different chemical reactions, a general method for classifying transition states as either early or late is needed. One can envision a dimensionless reaction coordinate that changes continuously and monotonically from 0 to 1 along a minimum energy reaction pathway. The value of the dimensionless reaction coordinate for the transition state (W(TS)) classifies transition states as (a) early when W(TS) < 0.5, (b) late when W(TS) > 0.5, and (c) equidistant between reactants and products when W(TS) = 0.5. In this article, we derive such a dimensionless reaction coordinate and illustrate its usefulness for several different chemical reactions.

Chemically Meaningful Atomic Charges That Reproduce the Electrostatic Potential in Periodic and Nonperiodic Materials.

Net atomic charges (NACs) can be used both to understand the chemical states of atoms in a material as well as to represent the electrostatic potential, V, of the material outside its electron distribution. However, many existing definitions of NACs have limitations that prevent them from adequately fulfilling this dual purpose. Some charge methods are not applicable to periodic materials or are inaccurate for systems containing buried atoms, while others work for both periodic and nonperiodic materials containing buried atoms but give NACS that do not accurately reproduce V. We present a new approach, density derived electrostatic and chemical (DDEC) charges, that overcomes these limitations by simultaneously optimizing the NACs to be chemically meaningful and to reproduce V outside the electron distribution. This atoms-in-molecule method partitions the total electron density among atoms and uses a distributed multipole expansion to formally reproduce V exactly outside the electron distribution. We compare different methods for computing NACs for a broad range of materials that are periodic in zero, one, two, and three dimensions. The DDEC method consistently performs well for systems with and without buried atoms, including molecules, nonporous solids, solid surfaces, and porous solids like metal organic frameworks.

The nature of aryloxide and arylsulfide ligand bonding in dimethyltitanium complexes containing cyclopentadienyl ligation.

A series of dimethyltitanium compounds [CpTi(EAr)Me2](E = O, S) ligated by one cyclopentadienyl (Cp) and one aryloxide (OAr) or arylsulfide (SAr) have been structurally characterized in order to gain a better understanding of aryloxide and arylsulfide bonding in these systems. Experimental structures were compared to those predicted by density functional theory (DFT). Bonding in the arylsulfide systems was found to be significantly different from bonding in the aryloxide systems. The aryloxide ligands exhibited wide Ti-O-Ar angles > or = 150 degrees) with the Ar group oriented proximal to the Cp group. DFT computations revealed two conformers for the arylsulfide systems. Arylsulfides with the Ar group proximal to the Cp group had a predicted Ti-S-Ar angle of approximately 120 degrees while those with the Ar group distal to the Cp group had a measured and predicted Ti-S-Ar angle of approximately 100 degrees. Molecular and natural bond orbital (NBO) analyses were employed to explain the nature of ligand bonding in these systems.