An efficient method for computing the QTAIM topology of a scalar field: The electron density case

An efficient method for computing the quantum theory of atoms in molecules (QTAIM) topology of the electron density (or other scalar field) is presented. A modified Newton–Raphson algorithm was implemented for finding the critical points (CP) of the electron density. Bond paths were constructed with the second‐order Runge–Kutta method. Vectorization of the present algorithm makes it to scale linearly with the system size. The parallel efficiency decreases with the number of processors (from 70% to 50%) with an average of 54%. The accuracy and performance of the method are demonstrated by computing the QTAIM topology of the electron density of a series of representative molecules. Our results show that our algorithm might allow to apply QTAIM analysis to large systems (carbon nanotubes, polymers, fullerenes) considered unreachable until now. © 2012 Wiley Periodicals, Inc.     

An algorithm to delineate and integrate topological basins in a three-dimensional quantum mechanical density function

The growing activity in the area of Quantum Chemical Topology warrants a new algorithm to delineate topological basins in 3D scalar fields other than the electron density. A method based on the "octal tree search algorithm" of computer graphics is proposed to reach this goal. We illustrate the algorithm on the L(r) function, which is the negative of the Laplacian of the electron density. Because of its complicated topology, even in a simple test molecule such as water, it benefits from the octal tree algorithm as a robust, compact, and general technique to find the boundaries of topological basins. For the first time, we are able to compute the population and volume of the core and valence (bonding and nonbonding, i.e., lone pair) basins given by L(r)'s topology.     

Shared and closed-shell O-O interactions in silicates

Bond paths of maximum electron density spanning O-O edges shared between equivalent or quasiequivalent MOn (n > 4) coordination polyhedra are not uncommon electron density features displayed by silicates. On the basis of the positive values for the local electronic energy density, H(rc), at the bond critical points, rc, they qualify as weak "closed-shell" interactions. As observed for M-O bonded interactions (M = first and second row metal atoms), the electron density, rho(rc), and the Laplacian of the electron density increase in a regular way as the separation between the O atoms, R(O-O), decreases. A simple model, based on R(O-O) and the distances of the Si atoms from the midpoint between adjacent pairs of O atoms, partitions the O-O bond paths in the high-pressure silica polymorph coesite into two largely disjoint domains, one with and one without bond paths. The occurrence of O-O bond paths shared in common between equivalent coordination polyhedra suggests that they may be grounded in some cases on factors other than bonded interactions, particularly since they are often displayed by inert procrystal representations of the electron density. In these cases, it can be argued that the accumulation of the electron density along the paths has its origin, at least in part, in the superposition of the peripheral electron density distributions of the metal M atoms occupying the edge-sharing polyhedra. On the other hand, the accumulation of electron density along the paths may stabilize a structure by shielding the adjacent M atoms in the edge-sharing polyhedra. For closed-shell Li-O, Na-O, and Mg-O interactions, H(rc) is positive and increases as the value of rho(rc) increases, unlike the "shared" Be-O, B-O, C-O, Al-O, Si-O, P-O, and S-O interactions, where H(rc) is negative and decreases as rho(rc) increases. The H(rc) values for the weak closed-shell O-O interactions also increase as rho(rc) increases, as observed for the closed-shell M-O interactions. On the basis of the bond critical point properties and the negative H(rc) value, the O-O interaction comprising the O2 molecule in silica III qualifies as a shared interaction.     

New insights in quantum chemical topology studies using numerical grid-based analyses

New insights in Quantum Chemical Topology of one-electron density functions have been proposed here by using a recent grid-based algorithm (Tang et al., J Phys Condens Matter 2009, 21, 084204), initially designed for the decomposition of the electron density. Beyond the charge analysis, we show that this algorithm is suitable for different scalar functions showing a more complex topology, that is, the Laplacian of the electron density, the electron localization function (ELF), and the molecular electrostatic potential (MEP). This algorithm makes use of a robust methodology enabling to numerically assign the data points of three-dimensional grids to basin volumes, and it has the advantage of requiring only the values of the scalar function without details on the wave function used to build the grid. Our implementation is briefly outlined (program named TopChem), its capabilities are examined, and technical aspects in terms of CPU requirement and accuracy of the results are discussed. Illustrative examples for individual molecules and crystalline solids obtained with gaussian and plane-wave-based density functional theory calculations are presented. Special attention was given to the MEP because its topological analysis is complex and scarce.     

On the similarity between molecular electron densities, electrostatic potentials and bare nuclear potentials

The similarities among the molecular contours of three scalar fields, viz. electron density (ED), electrostatic potential (ESP) and bare nuclear potential (BNP) have been investigated. The topological resemblance between ESP and ED contour diagrams (as prompted by the Thomas-Fermi model) is more pronounced than that for BNP and ED contour diagrams (as predicted by the local density functional model of Parr, Gadre and Bartolotti) with three-membered ring systems as test cases. An analysis of critical points of these distributions has also been included. Thus it may be conjectured that ED maps may prove useful in predicting reactive sites in molecules.     

Rapid topography mapping of scalar fields: Large molecular clusters

An efficient and rapid algorithm for topography mapping of scalar fields, molecular electron density (MED) and molecular electrostatic potential (MESP) is presented. The highlight of the work is the use of fast function evaluation by Deformed-atoms-in-molecules (DAM) method. The DAM method provides very rapid as well as sufficiently accurate function and gradient evaluation. For mapping the topography of large systems, the molecular tailoring approach (MTA) is invoked. This new code is tested out for mapping the MED and MESP critical points (CP's) of small systems. It is further applied to large molecular clusters viz. (H(2)O)(25), (C(6)H(6))(8) and also to a unit cell of valine crystal at MP2∕6-31+G(d) level of theory. The completeness of the topography is checked by extensive search as well as applying the Poincare?-Hopf relation. The results obtained show that the DAM method in combination with MTA provides a rapid and efficient route for mapping the topography of large molecular systems.     

Electron density distributions calculated for the nickel sulfides millerite, vaesite, and heazlewoodite and nickel metal: a case for the importance of ni-ni bond paths for electron transport

Bond paths and the bond critical point properties (the electron density (rho) and the Hessian of rho at the bond critical points (bcp's)) have been calculated for the bonded interactions comprising the nickel sulfide minerals millerite, NiS, vaesite, NiS(2), and heazlewoodite, Ni(3)S(2), and Ni metal. The experimental Ni-S bond lengths decrease linearly as the magnitudes of the properties each increases in value. Bond paths exist between the Ni atoms in heazlewoodite and millerite for the Ni-Ni separations that match the shortest separation in Ni metal, an indicator that the Ni atoms are bonded. The bcp properties of the bonded interactions in Ni metal are virtually the same as those in heazlewoodite and millerite. Ni-Ni bond paths are absent in vaesite where the Ni-Ni separations are 60% greater than those in Ni metal. The bcp properties for the Ni-Ni bonded interactions scatter along protractions of the Ni-S bond length-bcp property trends, suggesting that the two bonded interactions have similar characteristics. Ni-Ni bond paths radiate throughout Ni metal and the metallic heazlewoodite structures as continuous networks whereas the Ni-Ni paths in millerite, a p,d-metal displaying ionic and covalent features, are restricted to isolated Ni(3) rings. Electron transport in Ni metal and heazlewoodite is pictured as occurring along the bond paths, which behave as networks of atomic size wires that radiate in a contiguous circuit throughout the two structures. Unlike heazlewoodite, the electron transport in millerite is pictured as involving a cooperative hopping of the d-orbital electrons from the Ni(3) rings comprising Ni(3)S(9) clusters to Ni(3) rings in adjacent clusters via the p-orbitals on the interconnecting S atoms. Vaesite, an insulator at low temperatures and a doped semiconductor at higher temperatures, lacks Ni-Ni bond paths. The net charges conferred on the Ni and S atoms are about a quarter of their nominal charges for the atoms in millerite and vaesite with the net charge on Ni increasing with increasing Ni-S bond length. Reduced net charges are observed on the Ni atoms in heazlewoodite and are related to its Ni-Ni metal bonded interactions and to the greater covalent character of its bonds. Local energy density and bond critical point properties of the electron density distributions indicate that the Ni-S and Ni-Ni bonded interactions are intermediate in character between ionic and covalent.     

A recursive equation method for the determination of density and heat capacity: Comparison between isentropic and isothermal integration paths

Among different ways to obtain the thermodynamic properties of a fluid in the liquid phase, there is the possibility to solve, numerically, a system of two differential equations where density and specific heat capacity are the unknown variables expressed as functions of temperature and pressure. By means of general methods, like Predictor–Corrector or Runge–Kutta, this system can be solved integrating along isothermal paths, once the initial conditions are known along an isobar and the speed of sound values are measured over the $p''–''T$ region of interest. In this work a new perturbative method, called recursive equation method (REM), based on the analytic recursive determination of the coefficients describing density and specific heat capacity, is proposed. The main advantage of REM is the possibility to also get the value of the uncertainty associated with the solutions, obtained by means of standard methods of error analysis. Another improvement introduced by this algorithm is the possibility to solve the system of equations following arbitrary integration paths. In particular, a comparison between the results of density and heat capacity, obtained by integrating along isentropic and isothermal paths, is presented, both for water and for acetone.     

GPUs as boosters to analyze scalar and vector fields in quantum chemistry

The analysis of scalar and vector fields in quantum chemistry is an essential task for the computational chemistry community, where such quantities must be evaluated rapidly to perform a particular study. For example, the atoms in molecules approach proposed by Bader has become popular; however, this method demands significant computational resources to compute the involved tasks in short times. In this article, we discuss the importance of graphics processing units (GPU) to analyze electron density, and related fields, implementing several scalar, and vector fields within the graphics processing units for atoms and molecules (GPUAM) code developed by a group of the Universidad Autónoma Metropolitana in México City. With this application, the quantum chemistry community can perform demanding computational tasks on a desktop, where CPUs and GPUs are used to their maximum capabilities. The performance of GPUAM is tested in several systems and over different GPUs, where a GPU installed in a workstation converts it to a robust high-performance computing system.     

The vorticity of the current density tensor and 3D-aromaticity

The induced current density ( J ( r )) vector field has been extensively used as a criterion of aromaticity and electron delocalization. However, the selection of the direction of the perturbing magnetic field ( B ) is arbitrary and in the case of three-dimensional electron delocalization/aromaticity the selection could be ambiguous. The J ( r ) has also recently received some criticism as an aromaticity index. We propose the Trace of the Vorticity of the Current Density tensor (TVCD) scalar field as a more suitable quantity for the evaluation of electron delocalization of three-dimensional systems. It does not depend on the orientation of B and contrary to other related scalar fields like the anisotropy of the induced current density, it does not lose the information of the direction of the currents. We show that not only the currents parallel to the molecular plane are important in the evaluation of the aromaticity, as is largely believed, but also the perpendicular ones, which information is included in the TVCD. The TVCD is very useful to study planar and 3D delocalized molecules (eg, fullerenes). Moreover, the integration of the TVCD over an internal surface of the 3D-cages serves as index for 3D-aromaticity.     

Characterizing the oxygen-oxygen interaction in the dinitramide anion

Atomic interactions between oxygen atoms have been analyzed in terms of the Theory of Atoms in Molecules for the biguanidium dinitramide and biguanidium bis-dinitramide crystals. The electron density has been derived from X-ray diffraction data obtained at 90 K, and the potential energy density has been calculated using the density functional approach. Bond critical points have been found on the O(1)...O(4) interatomic line in both the electron density and potential energy density gradient fields. The bond path and its associated virial path have been obtained. The interaction has been identified as a bonding closed-shell type interaction.     

Electron probability distribution in AIM and ELF basins

The probabilities of finding a certain number of electrons enclosed in a given volume is calculated and discussed for a series of molecules. Two different methodologies to do the partition of the molecular space in separate volumes are investigated: the Atoms in the Molecules, AIM, topologic analysis of the density, and the topologic analysis of the Electron Localization Function (ELF). The formulas to calculate the probability distribution are reviewed and the way to implement them shortly explained. For a series of molecules, we present how the probability distribution complement the chemical information about the localization of the electrons in certain regions of the space. The calculations show that the probability of finding Z electrons in the AIM atomic basin associated to an atom of atomic number Z is, in general, low, even when the average number of electrons is close to Z. The probability distribution on the ELF basins associated to bonds yields new insight about the nature of the respective bond.     

DAMQT 2.1.0: A new version of the DAMQT package enabled with the topographical analysis of electron density and electrostatic potential in molecules

DAMQT‐2.1.0 is a new version of DAMQT package which includes topographical analysis of molecular electron density (MED) and molecular electrostatic potential (MESP), such as mapping of critical points (CPs), creating molecular graphs, and atomic basins. Mapping of CPs is assisted with algorithmic determination of Euler characteristic in order to provide a necessary condition for locating all possible CPs. Apart from the mapping of CPs and determination of molecular graphs, the construction of MESP‐based atomic basin is a new and exclusive feature introduced in DAMQT‐2.1.0. The GUI in DAMQT provides a user‐friendly interface to run the code and visualize the final outputs. MPI libraries have been implemented for all the tasks to develop the parallel version of the software. Almost linear scaling of computational time is achieved with the increasing number of processors while performing various aspects of topography. A brief discussion of molecular graph and atomic basin is provided in the current article highlighting their chemical importance. Appropriate example sets have been presented for demonstrating the functions and efficiency of the code. © 2015 Wiley Periodicals, Inc.     

The effect of electron correlation on the topological and atomic properties of the electron density distributions of molecules

The topological properties of the electron density and the properties of an atom in a molecule are calculated by means of second-order Møller-Plesset perturbation theory (MP2) and compared with the results of configuration interaction calculations (C12) which include all single and double substitutions from the Hartree-Fock reference configuration. A software package for analyzing the effects of electron correlation on the topological properties of the electron density of molecules is described. H₂CO is used to provide a numerical example and to indicate that the number of bond critical points is unaffected by the inclusion of electron correlation. Correlation leads to only a small shift in the positions of bond critical points and a small change in the electron density at bond critical points. It is further shown that the energy of an atom in a molecule can be calculated to an accuracy of 1 kcal/mol and the electron population of an atom to about 0.001e. A statistical method is used to show that the deviation of the MP2 correlation correction relative to the CI2 correlation correction for a variety of atomic properties is about 25%.     

A nonadditive methanol force field: bulk liquid and liquid-vapor interfacial properties via molecular dynamics simulations using a fluctuating charge model

We study the bulk and interfacial properties of methanol via molecular dynamics simulations using a CHARMM (Chemistry at HARvard Molecular Mechanics) fluctuating charge force field. We discuss the parametrization of the electrostatic model as part of the ongoing CHARMM development for polarizable protein force fields. The bulk liquid properties are in agreement with available experimental data and competitive with existing fixed-charge and polarizable force fields. The liquid density and vaporization enthalpy are determined to be 0.809 g/cm3 and 8.9 kcal/mol compared to the experimental values of 0.787 g/cm3 and 8.94 kcal/mol, respectively. The liquid structure as indicated by radial distribution functions is in keeping with the most recent neutron diffraction results; the force field shows a slightly more ordered liquid, necessarily arising from the enhanced condensed phase electrostatics (as evidenced by an induced liquid phase dipole moment of 0.7 D), although the average coordination with two neighboring molecules is consistent with the experimental diffraction study as well as with recent density functional molecular dynamics calculations. The predicted surface tension of 19.66+/-1.03 dyn/cm is slightly lower than the experimental value of 22.6 dyn/cm, but still competitive with classical force fields. The interface demonstrates the preferential molecular orientation of molecules as observed via nonlinear optical spectroscopic methods. Finally, via canonical molecular dynamics simulations, we assess the model's ability to reproduce the vapor-liquid equilibrium from 298 to 423 K, the simulation data then used to obtain estimates of the model's critical temperature and density. The model predicts a critical temperature of 470.1 K and critical density of 0.312 g/cm3 compared to the experimental values of 512.65 K and 0.279 g/cm3, respectively. The model underestimates the critical temperature by 8% and overestimates the critical density by 10%, and in this sense is roughly equivalent to the underlying fixed-charge CHARMM22 force field.     

PathOpt—A global transition state search approach: Outline of algorithm

We propose a new algorithm to determine reaction paths and test its capability for Ar₁₂ and Ar₁₃ clusters. Its main ingredient is a search for the local minima on a (n?1) dimensional hyperplane (n = dimension of the complete system in Cartesian coordinates) lying perpendicular to the straight line connection between initial and final states. These minima are part of possible reaction paths and are, hence, used as starting points for an uphill search to the next transition state. First, path fragments are obtained from subsequent relaxations starting from these transition states. They can be combined with information from the straight line connection procedure to obtain complete paths. Our test computations for Ar₁₂ and Ar₁₃ clusters prove that PathOpt delivers several reaction paths in one round. © 2013 Wiley Periodicals, Inc.     

Fast generation of molecular surfaces from 3D data fields with an enhanced “marching cube” algorithm

An improved version of the “marching cubes” algorithm [W. Lorensen and H. Cline, Comp. Graph. 21 , (1987)] for the generation of isosurfaces from 3D data fields is presented and applied to molecular surfaces. The new algorithm avoids inconsistent pattern definitions of the original one, which lead to artificial gaps. The advantage of a logarithmic interpolation procedure, in particular for data fields typically occurring in molecular science, is demonstrated. An example is the generation of molecular surfaces based upon electron density data. © 1993 John Wiley & Sons, Inc.     

Critical point representations of electron density maps for the comparison of benzodiazepine-type ligands

A procedure for the comparison of three-dimensional electron density distributions is proposed for similarity searches between pharmacological ligands at various levels of crystallographic resolution. First, a graph representation of molecular electron density distributions is generated using a critical point analysis approach. Pairwise as well as multiple comparisons between the obtained graphs of critical points are then carried out using a Monte Carlo/simulated annealing technique, and results are compared with genetic algorithm solutions.     

Topology and Chemistry

The determinants of chemical bonding are the chemical properties of the atoms and the constraints of three-dimensional (3-D) space into which the atoms must fit, but topology provides a convenient way of describing the resultant structure. This paper explores the topologies of various scalar fields associated with atoms in molecules and crystals and what they can tell us about chemical bonding. The scalar fields examined are the electron density, the electrostatic potential, and two simplified electrostatic potentials in which the contributions of the electron cores have been removed, namely the Madelung and the covalent field. Not all of the information contained in these fields is present in the topology but, since the topology is insensitive to the details of the field, it can often be determined using simplified calculations. Although the same topological model is used to explore all four fields, each field has its own distinctive topology and each provides different information about the nature of chemical bonding and structure. The analysis of these topologies, when combined with simple electrostatic theory and a few empirical observations, leads to a quantitative model of localized chemical bonding. In the process, the analysis provides insights into the nature of chemical bonding.     

Dimensional scaling treatment with relativistic corrections for stable multiply charged atomic ions in high-frequency super-intense laser fields

We present a theoretical framework which describes multiply charged atomic ions, their stability within super-intense laser fields, and also lay corrections to the systems due to relativistic effects. Dimensional scaling calculations with relativistic corrections for systems: H, H(-), H(2 -), He, He(-), He(2 -), He(3 -) within super-intense laser fields were completed. Also completed were three-dimensional self consistent field calculations to verify the dimensionally scaled quantities. With the aforementioned methods the system's ability to stably bind "additional" electrons through the development of multiple isolated regions of high potential energy leading to nodes of high electron density is shown. These nodes are spaced far enough from each other to minimize the electronic repulsion of the electrons, while still providing adequate enough attraction so as to bind the excess electrons into orbitals. We have found that even with relativistic considerations these species are stably bound within the field. It was also found that performing the dimensional scaling calculations for systems within the confines of laser fields to be a much simpler and more cost-effective method than the supporting D = 3 SCF method. The dimensional scaling method is general and can be extended to include relativistic corrections to describe the stability of simple molecular systems in super-intense laser fields.