Low‐energy structures of benzene clusters with a novel accurate potential surface

The benzene‐benzene (Bz‐Bz) interaction is present in several chemical systems and it is known to be crucial in understanding the specificity of important biological phenomena. In this work, we propose a novel Bz‐Bz analytical potential energy surface which is fine‐tuned on accurate ab initio calculations in order to improve its reliability. Once the Bz‐Bz interaction is modeled, an analytical function for the energy of the clusters may be obtained by summing up over all pair potentials. We apply an evolutionary algorithm (EA) to discover the lowest‐energy structures of clusters (for ), and the results are compared with previous global optimization studies where different potential functions were employed. Besides the global minimum, the EA also gives the structures of other low‐lying isomers ranked by the corresponding energy. Additional ab initio calculations are carried out for the low‐lying isomers of and clusters, and the global minimum is confirmed as the most stable structure for both sizes. Finally, a detailed analysis of the low‐energy isomers of the n = 13 and 19 magic‐number clusters is performed. The two lowest‐energy isomers show S₆ and C₃ symmetry, respectively, which is compatible with the experimental results available in the literature. The structures reported here are all non‐symmetric, showing two central Bz molecules surrounded by 12 nearest‐neighbor monomers in the case of the five lowest‐energy structures. © 2015 Wiley Periodicals, Inc.     

Computational study of basis set and electron correlation effects on anapole magnetizabilities of chiral molecules

In the presence of a static, nonhomogeneous magnetic field, represented by the axial vector at the origin of the coordinate system and by the polar vector , assumed to be spatially uniform, the chiral molecules investigated in this paper carry an orbital electronic anapole, described by the polar vector . The electronic interaction energy of these molecules in nonordered media is a cross term, coupling and via , one third of the trace of the anapole magnetizability a_αβ tensor, that is, . Both and W^BC have opposite sign in the two enantiomeric forms, a fact quite remarkable from the conceptual point of view. The magnitude of predicted in the present computational investigation for five chiral molecules is very small and significantly biased by electron correlation contributions, estimated at the density functional level via three different functionals. © 2016 Wiley Periodicals, Inc.     

Structure,stability, and nature of bonding in carbon monoxide bound complexes (E = group 14 element; X = H,F, Cl,Br, I)

A density functional theory study is performed to predict the structures and stability of carbon monoxide (CO) bound (E = C, Si, Ge, Sn, Pb; X = H, F, Cl, Br, I) complexes. The possibility of bonding through both C‐ and O‐sides of CO is considered. Thermochemical analysis reveals that all the dissociation processes producing CO and are endothermic in nature whereas most of the dissociation reactions are endergonic in nature at room temperature. The nature of bonding in E? C/O bonds is analyzed via Wiberg bond index, natural population analysis, electron density, and energy decomposition analyses in conjunction with natural orbitals for chemical valence scheme. In comparison to C? O stretching frequency ( ) in free CO, while a red shift is noted in O‐side binding, the C‐side binding results in a blue shift in . The relative change in values in CO bound complexes on changing either E or X is scrutinized and possible explanation is provided in terms of polarization in the σ‐ and π‐orbitals and the relative strength of C→E or O→E σ‐donation and E→C or E→O π‐back‐donation. © 2016 Wiley Periodicals, Inc.     

FEREBUS: Highly parallelized engine for kriging training

FFLUX is a novel force field based on quantum topological atoms, combining multipolar electrostatics with IQA intraatomic and interatomic energy terms. The program FEREBUS calculates the hyperparameters of models produced by the machine learning method kriging. Calculation of kriging hyperparameters ( θ and p ) requires the optimization of the concentrated log‐likelihood . FEREBUS uses Particle Swarm Optimization (PSO) and Differential Evolution (DE) algorithms to find the maximum of . PSO and DE are two heuristic algorithms that each use a set of particles or vectors to explore the space in which is defined, searching for the maximum. The log‐likelihood is a computationally expensive function, which needs to be calculated several times during each optimization iteration. The cost scales quickly with the problem dimension and speed becomes critical in model generation. We present the strategy used to parallelize FEREBUS, and the optimization of through PSO and DE. The code is parallelized in two ways. MPI parallelization distributes the particles or vectors among the different processes, whereas the OpenMP implementation takes care of the calculation of , which involves the calculation and inversion of a particular matrix, whose size increases quickly with the dimension of the problem. The run time shows a speed‐up of 61 times going from single core to 90 cores with a saving, in one case, of ～98% of the single core time. In fact, the parallelization scheme presented reduces computational time from 2871 s for a single core calculation, to 41 s for 90 cores calculation. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Quantum Monte Carlo with very large multideterminant wavefunctions

An algorithm to compute efficiently the first two derivatives of (very) large multideterminant wavefunctions for quantum Monte Carlo calculations is presented. The calculation of determinants and their derivatives is performed using the Sherman–Morrison formula for updating the inverse Slater matrix. An improved implementation based on the reduction of the number of column substitutions and on a very efficient implementation of the calculation of the scalar products involved is presented. It is emphasized that multideterminant expansions contain in general a large number of identical spin‐specific determinants: for typical configuration interaction‐type wavefunctions the number of unique spin‐specific determinants ( ) with a non‐negligible weight in the expansion is of order . We show that a careful implementation of the calculation of the N_det ‐dependent contributions can make this step negligible enough so that in practice the algorithm scales as the total number of unique spin‐specific determinants, , over a wide range of total number of determinants (here, N_det up to about one million), thus greatly reducing the total computational cost. Finally, a new truncation scheme for the multideterminant expansion is proposed so that larger expansions can be considered without increasing the computational time. The algorithm is illustrated with all‐electron fixed‐node diffusion Monte Carlo calculations of the total energy of the chlorine atom. Calculations using a trial wavefunction including about 750,000 determinants with a computational increase of ～400 compared to a single‐determinant calculation are shown to be feasible. © 2016 Wiley Periodicals, Inc.     

Topological analysis of the electron delocalization range

The electron delocalization range function EDR( ) (Janesko et al., J. Chem. Phys. 2014, 141, 144104) quantifies the extent to which an electron at point in a calculated wavefunction delocalizes over distance d. This work shows how topological analysis distills chemically useful information out of the EDR. Local maxima (attractors) in the EDR occur in regions such as atomic cores, covalent bonds, and lone pairs where the wavefunction is dominated by a single orbital lobe. The EDR characterizes each attractor in terms of a delocalization length D and a normalization , which are qualitatively consistent with the size of the orbital lobe and the number of lobes in the orbital. Attractors identify the progressively more delocalized atomic shells in heavy atoms, the interplay of delocalization and strong (nondynamical) correlation in stretched and dissociating covalent bonds, the locations of valence and weakly bound electrons in anionic water clusters, and the chemistry of different reactive sites on metal clusters. Application to ammonia dissociation over silicon illustrates how this density‐matrix‐based analysis can give insight into realistic systems. © 2016 Wiley Periodicals, Inc.     

The performance of density functional and wavefunction‐based methods for 2D and 3D structures of Au10

The transition from 2D to 3D structures in small gold clusters occurs around 10 atoms. Density functional theory predicts a planar structure for in contrast to recent second‐order Møller–Plesset perturbation theory calculations, which predict a 3D arrangement. The validity of the use of single‐reference second‐order Møller–Plesset theory for near metallic systems remains, however, questionable. On the other hand, it is less than clear how well density functional approximations perform for such clusters. We, therefore, decided to carry out quantum chemical calculations for using a variety of different density functionals as well as wavefunction‐based methods including coupled cluster theory to compare the different energetically low lying 2D and 3D cluster isomers. The results are perhaps not encouraging showing that most computational methods do not predict correctly the energetic sequence of isomers compared to coupled cluster theory. As perturbative triple corrections in the coupled cluster treatment change the order in cluster stability, the onset of 2D to 3D transition in these gold clusters remains elusive. As expected, second‐order Møller–Plesset theory is not suitable for correctly describing such systems.     

Ab initio potential energy surface and vibration–rotation energy levels of germanium dicarbide,GeC2

The accurate ground‐state potential energy surface of germanium dicarbide, GeC₂, has been determined from ab initio calculations using the coupled‐cluster approach. The core–electron correlation, higher‐order valence‐electron correlation, and scalar relativistic effects were taken into account. The potential energy surface of GeC₂ was shown to be extraordinarily flat near the T‐shaped equilibrium configuration. The potential energy barrier to the linear CCGe configuration was predicted to be 1218 cm⁻¹. The vibration–rotation energy levels of some GeC₂ isotopologues were calculated using a variational method. The vibrational bending mode ν₃ was found to be highly anharmonic, with the fundamental wavenumber being only 58 cm⁻¹. Vibrational progressions due to this mode were predicted for the , , and states of GeC₂. © 2018 Wiley Periodicals, Inc.     

Energy decomposition analysis of cation–π, metal ion–lone pair,hydrogen bonded,charge‐assisted hydrogen bonded,and π–π interactions

This study probes the nature of noncovalent interactions, such as cation–π, metal ion–lone pair (M–LP), hydrogen bonding (HB), charge‐assisted hydrogen bonding (CAHB), and π–π interactions, using energy decomposition schemes—density functional theory (DFT)–symmetry‐adapted perturbation theory and reduced variational space. Among cation–π complexes, the polarization and electrostatic components are the major contributors to the interaction energy (IE) for metal ion–π complexes, while for onium ion–π complexes ( , , , and ) the dispersion component is prominent. For M–LP complexes, the electrostatic component contributes more to the IE except the dicationic metal ion complexes with H₂S and PH₃ where the polarization component dominates. Although electrostatic component dominates for the HB and CAHB complexes, dispersion is predominant in π–π complexes.Copyright © 2015 Wiley Periodicals, Inc.     

Horizontal vectorization of electron repulsion integrals

We present an efficient implementation of the Obara–Saika algorithm for the computation of electron repulsion integrals that utilizes vector intrinsics to calculate several primitive integrals concurrently in a SIMD vector. Initial benchmarks display a 2–4 times speedup with AVX instructions over comparable scalar code, depending on the basis set. Speedup over scalar code is found to be sensitive to the level of contraction of the basis set, and is best for quartets when l_D = 0 or , which makes such a vectorization scheme particularly suitable for density fitting. The basic Obara–Saika algorithm, how it is vectorized, and the performance bottlenecks are analyzed and discussed. © 2016 Wiley Periodicals, Inc.     

Aqueous acidities of primary benzenesulfonamides: Quantum chemical predictions based on density functional theory and SMD

Aqueous of selected primary benzenesulfonamides are predicted in a systematic manner using density functional theory methods and the SMD solvent model together with direct and proton exchange thermodynamic cycles. Some test calculations were also performed using high‐level composite CBS‐QB3 approach. The direct scheme generally does not yield a satisfactory agreement between calculated and measured acidities due to a severe overestimation of the Gibbs free energy changes of the gas‐phase deprotonation reaction by the used exchange‐correlation functionals. The relative values calculated using proton exchange method compare to experimental data very well in both qualitative and quantitative terms, with a mean absolute error of about 0.4 units. To achieve this accuracy, we find it mandatory to perform geometry optimization of the neutral and anionic species in the gas and solution phases separately, because different conformations are stabilized in these two cases. We have attempted to evaluate the effect of the conformer‐averaged free energies in the predictions, and the general conclusion is that this procedure is highly too costly as compared with the very small improvement we have gained. © 2015 Wiley Periodicals, Inc.     

An anisotropic coarse‐grained model based on Gay–Berne and electric multipole potentials and its application to simulate a DMPC bilayer in an implicit solvent model

In this work, we aim at optimizing the performance of the anisotropic GBEMP model, which adopts a framework by combining a Gay–Berne (GB) anisotropic potential with an electric multipole (EMP) potential, in simulating a DMPC lipid bilayer in an implicit solvent model. First, the Gay–Berne parameters were initially obtained by fitting to atomistic profiles of van der Waals interactions between homodimers of molecular fragments while EMP parameters was directly derived from the expansion of point multipoles at predefined EMP sites. Second, the GB and EMP parameters for DMPC molecule were carefully optimized to be comparable to AMBER atomistic model in the calculations of the dipole moments of DMPC monomers adopting different conformations as well as the nonbonded interactions between two DMPC molecules adopting different conformations and separated at various distances. Finally, the GB parameters for DMPC were slightly adjusted in simulating a 72 DMPC bilayer system so that our GBEMP model would be able to reproduce a few important structural properties, namely, thickness (), area per lipid ( ) and volume per lipid ( ). Meanwhile, the atomistic and experimental results for electron density profiles and order parameters were reproduced reasonably well by the GBEMP model, demonstrating the promising feature of GBEMP model in modeling lipid systems. Finally, we have shown that current GBEMP model is more efficient by a factor of about 25 than AMBER atomistic point charge model. © 2015 Wiley Periodicals, Inc.     

Accurate double many‐body expansion potential energy surface by extrapolation to the complete basis set limit and dynamics calculations for ground state of NH2

An accurate single‐sheeted double many‐body expansion potential energy surface is reported for the title system. A switching function formalism has been used to warrant the correct behavior at the and dissociation channels involving nitrogen in the ground and first excited states. The topographical features of the novel global potential energy surface are examined in detail, and found to be in good agreement with those calculated directly from the raw ab initio energies, as well as previous calculations available in the literature. The novel surface can be using to treat well the Renner–Teller degeneracy of the and states of . Such a work can both be recommended for dynamics studies of the reaction and as building blocks for constructing the double many‐body expansion potential energy surface of larger nitrogen/hydrogen‐containing systems. In turn, a test theoretical study of the reaction has been carried out with the method of quantum wave packet on the new potential energy surface. Reaction probabilities, integral cross sections, and differential cross sections have been calculated. Threshold exists because of the energy barrier (68.5 meV) along the minimum energy path. On the curve of reaction probability for total angular momentum J = 0, there are two sharp peaks just above threshold. The value of integral cross section increases quickly from zero to maximum with the increase of collision energy, and then stays stable with small oscillations. The differential cross section result shows that the reaction is a typical forward and backward scatter in agreement with experimental measurement result. © 2013 Wiley Periodicals, Inc.     

Distinguishing and quantifying the torquoselectivity in competitive ring‐opening reactions using the stress tensor and QTAIM

Currently the theories to explain and predict the classification of the electronic reorganization due to the torquoselectivity of a ring‐opening reaction cannot accommodate the directional character of the reaction pathway; the torquoselectivity is a type of stereoselectivity and therefore is dependent on the pathway. Therefore, in this investigation we introduced new measures from quantum theory of atoms in molecules and the stress tensor to clearly distinguish and quantify the transition states of the inward (TSIC) and outward (TSOC) conrotations of competitive ring‐opening reactions of 3‐(trifluoromethyl)cyclobut‐1‐ene and 1‐cyano‐1‐methylcyclobutene. We find the metallicity ξ( r _b) of the ring‐opening bond does not occur exactly at the transition state in agreement with transition state theory. The vector‐based stress tensor response β_σ was used to distinguish the effect of the CN, CH₃, and CF₃ groups on the TSIC and TSOC paths that was consistent with the ellipticity ε, the total local energy density H( r_b ) and the stress tensor stiffness S_σ. We determine the directional properties of the TSIC and TSOC ring‐opening reactions by constructing a stress tensor space with trajectories (s) with length l in real space, longer l correlated with the lowest density functional theory‐evaluated total energy barrier and hence will be more thermodynamically favored. © 2016 Wiley Periodicals, Inc.     

Structure and electrochemical properties for complexes of nitrocompounds with inorganic ions: A theoretical approach

Reduction and oxidation (redox) reactions are widely used for removal of nitrocompounds from contaminated soil and water. Structures and redox properties for complexes of nitrocompounds, such as 2,4,6‐trinitrotoluene (TNT), 2,4‐dinitrotoluene (DNT), 2,4‐dinitroanisole (DNAN), and 5‐nitro‐2,4‐dihydro‐3H?1,2,4‐triazol‐3‐one (NTO), with common inorganic ions (Na⁺, Cl^?, ) were investigated at the SMD(Pauling)/PCM(Pauling)/MPWB1K/TZVP level of theory. Atoms in molecules (AIM) theory was applied to analyze the topological properties of the bond critical points involved in the interactions between the nitrocompounds and the ions. Topological analyses show that intermolecular interactions of the types O(N)…Na⁺, C‐H…Cl^?( ), and C…Cl^?( ) may be discussed as noncovalent closed‐shell interactions, while N‐H···Cl^?( ) hydrogen bonds are partially covalent in nature. Complexation causes significant decrease of redox activity of the nitrocompounds. Analysis of the reduction potentials of the complexes obtained through application of the Pourbaix diagram of an iron/water system revealed that sodium complexes of NTO might be reduced by metallic iron. © 2016 Wiley Periodicals, Inc.     

Fast search algorithms for computational protein design

One of the main challenges in computational protein design (CPD) is the huge size of the protein sequence and conformational space that has to be computationally explored. Recently, we showed that state‐of‐the‐art combinatorial optimization technologies based on Cost Function Network (CFN) processing allow speeding up provable rigid backbone protein design methods by several orders of magnitudes. Building up on this, we improved and injected CFN technology into the well‐established CPD package Osprey to allow all Osprey CPD algorithms to benefit from associated speedups. Because Osprey fundamentally relies on the ability of to produce conformations in increasing order of energy, we defined new strategies combining CFN lower bounds, with new side‐chain positioning‐based branching scheme. Beyond the speedups obtained in the new ‐CFN combination, this novel branching scheme enables a much faster enumeration of suboptimal sequences, far beyond what is reachable without it. Together with the immediate and important speedups provided by CFN technology, these developments directly benefit to all the algorithms that previously relied on the DEE/ combination inside Osprey* and make it possible to solve larger CPD problems with provable algorithms. © 2016 Wiley Periodicals, Inc.     

Polygonal current models for polycyclic aromatic hydrocarbons and graphene sheets of various shapes

Assuming that graphene is an “infinite alternant” polycyclic aromatic hydrocarbon resulting from tessellation of a surface by only six‐membered carbon rings, planar fragments of various size and shape (hexagon, triangle, rectangle, and rhombus) have been considered to investigate their response to a magnetic field applied perpendicularly. Allowing for simple polygonal current models, the diatropicity of a series of polycyclic textures has been reliably determined by comparing quantitative indicators, the π‐electron contribution to I^B, the magnetic field‐induced current susceptibility of the peripheral circuit, to and to , respectively the out‐of‐plane components of the magnetizability tensor and of the magnetic shielding tensor at the center of mass. Extended numerical tests and the analysis based on the polygonal model demonstrate that (i) and yield inadequate and sometimes erroneous measures of diatropicity, as they are heavily flawed by spurious geometrical factors, (ii) I^B values computed by simple polygonal models are valid quantitative indicators of aromaticity on the magnetic criterion, preferable to others presently available, whenever current susceptibility cannot be calculated ab initio as a flux integral, (iii) the hexagonal shape is the most effective to maximize the strength of π‐electron currents over the molecular perimeter, (iv) the edge current strength of triangular and rhombic graphene fragments is usually much smaller than that of hexagonal ones, (v) doping by boron and nitrogen nuclei can regulate and even inhibit peripheral ring currents, (vi) only for very large rectangular fragments can substantial current strengths be expected. © 2017 Wiley Periodicals, Inc.     

Free energies of solvation in the context of protein folding: Implications for implicit and explicit solvent models

Implicit solvent models for biomolecular simulations have been developed to use in place of more expensive explicit models; however, these models make many assumptions and approximations that are likely to affect accuracy. Here, the changes in free energies of solvation upon folding of several fast folding proteins are calculated from previously run μs–ms simulations with a number of implicit solvent models and compared to the values needed to be consistent with the explicit solvent model used in the simulations. In the majority of cases, there is a significant and substantial difference between the values calculated from the two approaches that is robust to the details of the calculations. These differences could only be remedied by selecting values for the model parameters—the internal dielectric constant for the polar term and the surface tension coefficient for the nonpolar term—that were system‐specific or physically unrealistic. We discuss the potential implications of our findings for both implicit and explicit solvent simulations. © 2015 Wiley Periodicals, Inc.     

Effects of zero point vibration on the reaction dynamics of water dimer cations following ionization

Reactions of water dimer cation following ionization have been investigated by means of a direct ab initio molecular dynamics method. In particular, the effects of zero point vibration and zero point energy (ZPE) on the reaction mechanism were considered in this work. Trajectories were run on two electronic potential energy surfaces (PESs) of : ground state (²A″‐like state) and the first excited state (²A′ ‐ like state). All trajectories on the ground‐state PES lead to the proton‐transferred product: H₂O⁺(Wd)‐H₂O(Wa) → OH(Wd)‐H₃O⁺(Wa), where Wd and Wa refer to the proton donor and acceptor water molecules, respectively. Time of proton transfer (PT) varied widely from 15 to 40 fs (average time of PT = 30.9 fs). The trajectories on the excited‐state PES gave two products: an intermediate complex with a face‐to‐face structure (H₂O‐OH₂)⁺ and a PT product. However, the proton was transferred to the opposite direction, and the reverse PT was found on the excited‐state PES: H₂O(Wd)‐H₂O⁺ (Wa) → H₃O⁺(Wd)‐OH(Wa). This difference occurred because the ionizing water molecule in the dimer switched between the ground and excited states. The reaction mechanism of and the effects of ZPE are discussed on the basis of the results. © 2017 Wiley Periodicals, Inc.