As an extension of the Outlier FLOODing (OFLOOD) method [Harada et al., J. Comput. Chem. 2015, 36, 763], the sparsity of the outliers defined by a hierarchical clustering algorithm, FlexDice, was considered to achieve an efficient conformational search as sparsity-weighted “OFLOOD.” In OFLOOD, FlexDice detects areas of sparse distribution as outliers. The outliers are regarded as candidates that have high potential to promote conformational transitions and are employed as initial structures for conformational resampling by restarting molecular dynamics simulations. When detecting outliers, FlexDice defines a rank in the hierarchy for each outlier, which relates to sparsity in the distribution. In this study, we define a lower rank (first ranked), a medium rank (second ranked), and the highest rank (third ranked) outliers, respectively. For instance, the first-ranked outliers are located in a given conformational space away from the clusters (highly sparse distribution), whereas those with the third-ranked outliers are nearby the clusters (a moderately sparse distribution). To achieve the conformational search efficiently, resampling from the outliers with a given rank is performed. As demonstrations, this method was applied to several model systems: Alanine dipeptide, Met-enkephalin, Trp-cage, T4 lysozyme, and glutamine binding protein. In each demonstration, the present method successfully reproduced transitions among metastable states. In particular, the first-ranked OFLOOD highly accelerated the exploration of conformational space by expanding the edges. In contrast, the third-ranked OFLOOD reproduced local transitions among neighboring metastable states intensively. For quantitatively evaluations of sampled snapshots, free energy calculations were performed with a combination of umbrella samplings, providing rigorous landscapes of the biomolecules. © 2015 Wiley Periodicals, Inc.

Biologically rare events play important roles in understanding functions. To computationally reproduce them, Outlier FLOODing (OFLOOD) method is powerful, in which sparse distributions of biological states are detected as outliers and intensively resampled by MD simulations. As an extension, sparsity-weighted OFLOOD method is newly proposed, in which a hierarchical clustering defines ranks of outliers. Accordingly to the ranks, the confirmational resampling from outliers is performed, accerlarating the conformational sampling of bio-molecules.

The UV-induced photochemistry of HCFC-132b (CF2ClCH2Cl) was investigated by computing excited-state properties with time-dependent density functional theory (TDDFT), multiconfigurational second-order perturbation theory (CASPT2), and coupled cluster with singles, doubles, and perturbative triples (CCSD(T)). Excited states calculated with TDDFT show good agreement with CASPT2 and CCSD(T) results, correctly predicting the main excited-states properties. Simulations of ultrafast nonadiabatic dynamics in the gas phase were performed, taking into account 25 electronic states at TDDFT level starting in two different spectral windows (8.5 ± 0.25 and 10.0 ± 0.25 eV). Experimental data measured at 123.6 nm (10 eV) is in very good agreement with our simulations. The excited-state lifetimes are 106 and 191 fs for the 8.5 and 10.0 eV spectral windows, respectively. Internal conversion to the ground state occurred through several different reaction pathways with different products, where 2Cl, C-Cl bond breakage, and HCl are the main photochemical pathways in the low-excitation region, representing 95% of all processes. On the other hand, HCl, HF, and C-Cl bond breakage are the main reaction pathways in the higher excitation region, with 77% of the total yield. © 2015 Wiley Periodicals, Inc.

HCFC-132b is an important industrial compound, with a strong impact on health and environment. Upon UV irradiation, it decomposes into dozens of different photoproducts. In this article, nonadiabatic dynamics simulation is used to explain how photo-decomposition takes place through the competition between diverse reaction pathways in the subpicosecond time scale.

On page 78 (DOI: 10.1002/jcc.24021), Ramon Carbó-Dorca discusses aromaticity, quantummultimolecular polyhedral, and the quantumQSPR fundamental equation. First, a concise description of the Kekulé's historical origin of aromaticity and the actual state of the question is given. After this, it is argued that still room is left to the discussion about the quantummechanical foundation existence of aromaticity. In order to perform that, quantum multimolecular polyhedra (QMP) are defined: they are based onmolecular density functions sets attached to QMP vertices. Fromthere, collective QMP distances, QSPR fundamental equation and aromaticity descriptors are proposed as away to construct an equation, able to estimate aromaticity via expectation values of Hermitian operators.

On page 18 (DOI: 10.1002/jcc.23914), the Gibbs energies of association between primary alkyl ammonium ions and crown ethers in solution are measured and calculated by Andreas J. Achazi, Larissa K. S. von Krbek, Christoph A. Schalley, and Beate Paulus.Measurementswere carried out by isothermal titration calorimetry. Calculations were done as accurate as possible for the gas phasewith DFT-D3(BJ). The Gibbs energies to transfer the educts in the gas phase and the products back in solution were calculated with the solvation model COSMO-RS in order to get the Gibbs energies of association in solution. Calculated andmeasured Gibbs energies of association in solution agreewell and reveal a strong solvent-dependent ion pair effect.

The cover image shows the largestof a family of fullerenes used for extrapolating to the graphene limit, as presented by Lukas N. Wirz, Ralf Tonner, Andreas Hermann, Rebecca Sure, and Peter Schwerdtfeger on page 10 (DOI: 10.1002/jcc.23894). The structures were obtained from a newly developed force field treated subsequently by density functional theory. Our results confirm Paul von Ragué Schleyer's hypothesis that C60 is not especially stable – 60 is not amagic number – comparedwith other fullerenes.

We present new algorithms to improve the performance of ENUF method (F. Hedman, A. Laaksonen, Chem. Phys. Lett. 425, 2006, 142) which is essentially Ewald summation using Non-Uniform FFT (NFFT) technique. A NearDistance algorithm is developed to extensively reduce the neighbor list size in real-space computation. In reciprocal-space computation, a new algorithm is developed for NFFT for the evaluations of electrostatic interaction energies and forces. Both real-space and reciprocal-space computations are further accelerated by using graphical processing units (GPU) with CUDA technology. Especially, the use of CUNFFT (NFFT based on CUDA) very much reduces the reciprocal-space computation. In order to reach the best performance of this method, we propose a procedure for the selection of optimal parameters with controlled accuracies. With the choice of suitable parameters, we show that our method is a good alternative to the standard Ewald method with the same computational precision but a dramatically higher computational efficiency. © 2015 Wiley Periodicals, Inc.

By using graphical processing units and developing new algorithms both for real-space and reciprocal-space computations, the computational efficiency of the ENUF method used to evaluate electrostatic interactions in particle-based simulations is extraordinarily improved. A procedure for the selection of the optimal parameters with controlled accuracies is proposed. With this method, the same computational precision as the standard Ewald method is achieved, but with dramatically higher computational efficiency.

We respond to the two questions posed by Weinhold, Schleyer, and McKee (WSM) in their study of cis-2-butene (Weinhold et al., J Comput Chem 2014, 35, 1499), in which they solicit explanations for the relative conformational energies of this molecule in terms of the Quantum Theory of Atoms in Molecules (QTAIM). WSM requested answers to the questions: (1) why is cis-2-butene less stable than trans-2-butene despite the presence of a hydrogen-hydrogen (H⋯H) bond path in the former but not in the latter if the H⋯H bond path is stabilizing? (2) Why is the potential well of the conformational global minimum of cis-2-butene only 0.8 kcal/mol deep when the H⋯H bonding is stabilizing by 5 kcal/mol? Both questions raised by WSM are answered by considering the changes in the energies of all atoms as a function of the rotation of one of the two methyl groups from the minimum-energy structure, which exhibits the H⋯H bond path, to the transition state, which is devoid of this bond path. It is found that the stability gained by the H⋯H bonding interaction is cancelled by the destabilization of one of the ethylenic carbon atoms which, alone, destabilizes the system by as much as 5 kcal/mol in the global minimum conformation. Further, it is found that the 1.1 kcal/mol stability of trans-2-butene with respect to the cis-isomer is driven by the considerable destabilization of the ethylenic carbons by 11 kcal/mol, while the changes in the atomic energies of the other corresponding atoms in the two isomers account for the observed different stabilities. The error introduced into QTAIM atomic energies by neglecting the virials of the forces on the nuclei for partially optimized structures is discussed. © 2015 Wiley Periodicals, Inc.

Atomic origin of the locally stabilizing H⋯H contact in cis-2-butene from virial QTAIM atomic energies along the potential energy surface of methyl rotation, i.e., in terms of atomic sub-potential energy surfaces.

We develop a new method for calculating the hydration free energy (HFE) of a protein with any net charge. The polar part of the energetic component in the HFE is expressed as a linear combination of four geometric measures (GMs) of the protein structure and the generalized Born (GB) energy plus a constant. The other constituents in the HFE are expressed as linear combinations of the four GMs. The coefficients (including the constant) in the linear combinations are determined using the three-dimensional reference interaction site model (3D-RISM) theory applied to sufficiently many protein structures. Once the coefficients are determined, the HFE and its constituents of any other protein structure are obtained simply by calculating the four GMs and GB energy. Our method and the 3D-RISM theory give perfectly correlated results. Nevertheless, the computation time required in our method is over four orders of magnitude shorter.

Although the hydration free energy (HFE) is one of the most important factors in studies on the structural stability of a protein, its calculation is significantly difficult in computational cost and accuracy. We develop a new method for calculating the HFE by combining the generalized Born model and the morphometric approach. Our method gives almost the same result as that from the three-dimensional reference interaction site model (3D-RISM) theory with drastic reduction of computational cost.

The Shift-and-invert parallel spectral transformations (SIPs), a computational approach to solve sparse eigenvalue problems, is developed for massively parallel architectures with exceptional parallel scalability and robustness. The capabilities of SIPs are demonstrated by diagonalization of density-functional based tight-binding (DFTB) Hamiltonian and overlap matrices for single-wall metallic carbon nanotubes, diamond nanowires, and bulk diamond crystals. The largest (smallest) example studied is a 128,000 (2000) atom nanotube for which ∼330,000 (∼5600) eigenvalues and eigenfunctions are obtained in ∼190 (∼5) seconds when parallelized over 266,144 (16,384) Blue Gene/Q cores. Weak scaling and strong scaling of SIPs are analyzed and the performance of SIPs is compared with other novel methods. Different matrix ordering methods are investigated to reduce the cost of the factorization step, which dominates the time-to-solution at the strong scaling limit. A parallel implementation of assembling the density matrix from the distributed eigenvectors is demonstrated. © 2015 Wiley Periodicals, Inc.

Massively parallel supercomputers can extend the size of systems that we can study with quantum chemistry methods. However, the eigenvalue problem remains the bottleneck of scalability for density-functional based tight-binding (DFTB) or semi-empirical molecular orbital methods. We present a sparse eigensolver that enables DFTB calculations for systems with more than 100,000 atoms utilizing more than 200,000 CPU cores.

A number of model Diels-Alder (D-A) cycloaddition reactions (H2CCH2 + cyclopentadiene and H2CCHX + 1,3-butadiene, with X = H, F, CH3, OH, CN, NH2, and NO) were studied by static (transition state - TS and IRC) and dynamics (quasiclassical trajectories) approaches to establish the (a)synchronous character of the concerted mechanism. The use of static criteria, such as the asymmetry of the TS geometry, for classifying and quantifying the (a)synchronicity of the concerted D-A reaction mechanism is shown to be severely limited and to provide contradictory results and conclusions when compared to the dynamics approach. The time elapsed between the events is shown to be a more reliable and unbiased criterion and all the studied D-A reactions, except for the case of H2CCHNO, are classified as synchronous, despite the gradual and quite distinct degrees of (a)symmetry of the TS structures. © 2015 Wiley Periodicals, Inc.

Model Diels-Alder cycloaddition reactions were studied by static and dynamics approaches to establish the (a)synchronous character of the concerted mechanism. The use of static criteria, such as the asymmetry of the TS geometry, for classifying and quantifying the (a)synchronicity of the concerted reaction mechanism provides contradictory results and conclusions when compared to the dynamics approach.

Fragment-based linear scaling quantum chemistry methods are a promising tool for the accurate simulation of chemical and biomolecular systems. Because of the coupled inter-fragment electrostatic interactions, a dual-layer iterative scheme is often employed to compute the fragment electronic structure and the total energy. In the dual-layer scheme, the self-consistent field (SCF) of the electronic structure of a fragment must be solved first, then followed by the updating of the inter-fragment electrostatic interactions. The two steps are sequentially carried out and repeated; as such a significant total number of fragment SCF iterations is required to converge the total energy and becomes the computational bottleneck in many fragment quantum chemistry methods. To reduce the number of fragment SCF iterations and speed up the convergence of the total energy, we develop here a new SCF scheme in which the inter-fragment interactions can be updated concurrently without converging the fragment electronic structure. By constructing the global, block-wise Fock matrix and density matrix, we prove that the commutation between the two global matrices guarantees the commutation of the corresponding matrices in each fragment. Therefore, many highly efficient numerical techniques such as the direct inversion of the iterative subspace method can be employed to converge simultaneously the electronic structure of all fragments, reducing significantly the computational cost. Numerical examples for water clusters of different sizes suggest that the method shall be very useful in improving the scalability of fragment quantum chemistry methods. © 2015 Wiley Periodicals, Inc.

Fragment-based, linear-scaling, quantum chemistry methods hold great potential for application in accurate simulations of complex molecular systems. Because of the coupled inter-fragment interactions, however, conventional fragment methods often employ a dual-layer SCF scheme to converge both the intra- and inter-fragment interactions. The total number of fragment SCF iterations displays a size-dependence to the target molecule. A new global SCF scheme developed here shows excellent performance over a broad range of molecular sizes and is applicable to other fragment methods.

We show that the empirical 18-electron rule of transition metal chemistry corresponds to an intrinsic saturation limit for the 3c/4e hyperbonding interactions that are a ubiquitous feature of D-block aggregation phenomena. Such a “rule” therefore requires no “p-orbital participation,” “d2sp3 hybridization,” “valence shell expansion,” or other p-type intrusions into the Aufbau orbital filling sequence. Instead, 18e stability corresponds to the natural terminus of post-Lewis 3c/4e resonance-type stabilizations that lead to successive ligand additions (and formal increments of electron count) at a transition metal center, all within the normal (s + 5d) confines of D-block valence space. © 2015 Wiley Periodicals, Inc.

Electronic iconography of the 18e rule, illustrating how a normal-valent (12e) transition metal MX3 species (using s and d orbitals at M to accommodate three MX bonds, blue dots, and three lone pairs, red dots) achieves the maximally stabilizing eM = 18 “count” (all dots) by adding three L donor ligands (green dots) that participate in all possible 3c/4e hyperbonding interactions (nLσ*MX) with MX bonds.

In this report, we summarize and describe the recent unique updates and additions to the Molcas quantum chemistry program suite as contained in release version 8. These updates include natural and spin orbitals for studies of magnetic properties, local and linear scaling methods for the Douglas–Kroll–Hess transformation, the generalized active space concept in MCSCF methods, a combination of multiconfigurational wave functions with density functional theory in the MC-PDFT method, additional methods for computation of magnetic properties, methods for diabatization, analytical gradients of state average complete active space SCF in association with density fitting, methods for constrained fragment optimization, large-scale parallel multireference configuration interaction including analytic gradients via the interface to the Columbus package, and approximations of the CASPT2 method to be used for computations of large systems. In addition, the report includes the description of a computational machinery for nonlinear optical spectroscopy through an interface to the QM/MM package Cobramm. Further, a module to run molecular dynamics simulations is added, two surface hopping algorithms are included to enable nonadiabatic calculations, and the DQ method for diabatization is added. Finally, we report on the subject of improvements with respects to alternative file options and parallelization. © 2015 Wiley Periodicals, Inc.

The Molcas quantum chemistry program package has a long history, and with the release of Molcas 8 in 2014, it continues to offer state-of-the-art tools for computational chemistry. This article summarizes some of the most significant additions and improvements included in the package in the last 6 years. There are sections on electron correlation methods, relativistic features, molecular dynamics, gradients and optimizations, and technical features.

The conformational dynamics of a macromolecule can be modulated by a number of factors, including changes in environment, ligand binding, and interactions with other macromolecules, among others. We present a method that quantifies the differences in macromolecular conformational dynamics and automatically extracts the structural features responsible for these changes. Given a set of molecular dynamics (MD) simulations of a macromolecule, the norms of the differences in covariance matrices are calculated for each pair of trajectories. A matrix of these norms thus quantifies the differences in conformational dynamics across the set of simulations. For each pair of trajectories, covariance difference matrices are parsed to extract structural elements that undergo changes in conformational properties. As a demonstration of its applicability to biomacromolecular systems, the method, referred to as DIRECT-ID, was used to identify relevant ligand-modulated structural variations in the β2-adrenergic (β2AR) G-protein coupled receptor. Micro-second MD simulations of the β2AR in an explicit lipid bilayer were run in the apo state and complexed with the ligands: BI-167107 (agonist), epinephrine (agonist), salbutamol (long-acting partial agonist), or carazolol (inverse agonist). Each ligand modulated the conformational dynamics of β2AR differently and DIRECT-ID analysis of the inverse-agonist vs. agonist-modulated β2AR identified residues known through previous studies to selectively propagate deactivation/activation information, along with some previously unidentified ligand-specific microswitches across the GPCR. This study demonstrates the utility of DIRECT-ID to rapidly extract functionally relevant conformational dynamics information from extended MD simulations of large and complex macromolecular systems. © 2015 Wiley Periodicals, Inc.

Characterizing changes in conformational dynamics of a macromolecule is important for understanding its function. To rapidly quantify the change and identify structural features that contribute most significantly to the conformational changes, a covariance matrix-based method, called DIRECT-ID, is presented. Ligand-modulated conformational changes in the β2-adrenergic GPCR are quantified using DIRECT-ID and both previously known and new ligand-specific microswitches are identified using the method.

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.

We compare changes in solvation free energy upon folding provided by several implicit solvent models and the TIP3P explicit solvent model. Inconsistencies of an unexpected magnitude were found across the models, which could only be corrected by using settings that were nonphysical or system-specific.

We present a low rank moment expansion of the linear density-density response function. The general interacting (fully nonlocal) density-density response function is calculated by means of its spectral decomposition via an iterative Lanczos diagonalization technique within linear density functional perturbation theory. We derive a unitary transformation in the space of the eigenfunctions yielding subspaces with well-defined moments. This transformation generates the irreducible representations of the density-density response function with respect to rotations within SO(3). This allows to separate the contributions to the electronic response density from different multipole moments of the perturbation. Our representation maximally condenses the physically relevant information of the density-density response function required for intermolecular interactions, yielding a considerable reduction in dimensionality. We illustrate the performance and accuracy of our scheme by computing the electronic response density of a water molecule to a complex interaction potential. © 2015 Wiley Periodicals, Inc.

Intermolecular interactions lead to changes of the electronic charge density. We provide an efficient scheme to express these changes for an arbitrary interaction with an universal yet low-rank tensor. To that extend, we transform the linear density–density response function from its spectral decomposition to a more condensed representation, separating the contributions to the electronic response density from different multipole moments of the perturbation.