Volume 140, Issue 21, 07 June 2014

Markov state models (MSMs) have been demonstrated to be a powerful method for computationally studying intramolecular processes such as protein folding and macromolecular conformational changes. In this article, we present a new approach to construct MSMs that is applicable to modeling a broad class of multimolecular assembly reactions. Distinct structures formed during assembly are distinguished by their undirected graphs, which are defined by strong subunit interactions. Spatial inhomogeneities of free subunits are accounted for using a recently developed Gaussianbased signature. Simplifications to this state identification are also investigated. The feasibility of this approach is demonstrated on two different coarsegrained models for virus selfassembly. We find good agreement between the dynamics predicted by the MSMs and long, unbiased simulations, and that the MSMs can reduce overall simulation time by orders of magnitude.
 COMMUNICATIONS


Communication: Kinetic and pairing contributions in the dielectric spectra of electrolyte solutions
View Description Hide DescriptionIn the late 1970s, Hubbard and Onsager predicted that adding salt to a polar solution would result in a reduced dielectric permittivity that arises from the unexpected tendency of solvent dipoles to align opposite to the applied field. Here we develop a novel nonequilibrium molecular dynamics simulation approach to determine this decrement accurately. Using a thermodynamic consistent allatom force field we show that for an aqueous solution containing sodium chloride around 4.8 mol/l, this effect accounts for 12% of the total dielectric permittivity. The dielectric decrement can be strikingly different if a less accurate force field for the ions is used. Using the widespread GROMOS parameters, we observe in fact an increment of the dielectric permittivity rather than a decrement, caused by ion pairing and introduced by a too low dispersion force.

Communication: Structure characterization of hard sphere packings in amorphous and crystalline states
View Description Hide DescriptionThe channel size distribution in hard sphere systems, based on the local neighbor correlation of four particle positions, is investigated for all volume fractions up to jamming. For each particle, all three particle combinations of neighbors define channels, which are relevant for the concept of caging. The analysis of the channel size distribution is shown to be very useful in distinguishing between gaseous, liquid, partially and fully crystallized, and glassy (random) jammed states. A common microstructural feature of four coplanar particles is observed in crystalline and glassy jammed states, suggesting the presence of “hidden” twodimensional order in threedimensional random close packings.

Communication: Integral equation theory for pair correlation functions in a crystal
View Description Hide DescriptionA method for calculating pair correlation functions in a crystal is developed. The method is based on separating the one and twoparticle correlation functions into the symmetry conserving and the symmetry broken parts. The conserving parts are calculated using the integral equation theory of homogeneous fluids. The symmetry broken part of the direct pair correlation function is calculated from a series written in powers of order parameters and that of the total pair correlation function from the OrnsteinZernike equation. The results found for a twodimensional hexagonal lattice show that the method provides accurate and detailed informations about the pair correlation functions in a crystal.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

Using Markov state models to study selfassembly
View Description Hide DescriptionMarkov state models (MSMs) have been demonstrated to be a powerful method for computationally studying intramolecular processes such as protein folding and macromolecular conformational changes. In this article, we present a new approach to construct MSMs that is applicable to modeling a broad class of multimolecular assembly reactions. Distinct structures formed during assembly are distinguished by their undirected graphs, which are defined by strong subunit interactions. Spatial inhomogeneities of free subunits are accounted for using a recently developed Gaussianbased signature. Simplifications to this state identification are also investigated. The feasibility of this approach is demonstrated on two different coarsegrained models for virus selfassembly. We find good agreement between the dynamics predicted by the MSMs and long, unbiased simulations, and that the MSMs can reduce overall simulation time by orders of magnitude.

Minima hopping guided path search: An efficient method for finding complex chemical reaction pathways
View Description Hide DescriptionThe Minima Hopping global optimization method uses physically realizable molecular dynamics moves in combination with an energy feedback that guarantees the escape from any potential energy funnel. For the purpose of finding reaction pathways, we argue that Minima Hopping is particularly suitable as a guide through the potential energy landscape and as a generator for pairs of minima that can be used as input structures for methods capable of finding transition states between two minima. For LennardJones benchmark systems we compared this Minima Hopping guided path search method to a known approach for the exploration of potential energy landscapes that is based on deterministic modefollowing. Although we used a stabilized modefollowing technique that reliably allows to follow distinct directions when escaping from a local minimum, we observed that Minima Hopping guided path search is far superior in finding lowestbarrier reaction pathways. We, therefore, suggest that Minima Hopping guided path search can be used as a simple and efficient way to identify energetically lowlying chemical reaction pathways. Finally, we applied the Minima Hopping guided path search approach to 75atom and 102atom LennardJones systems. For the 75atom system we found pathways whose highest energies are significantly lower than the highest energy along the previously published lowestbarrier pathway. Furthermore, many of these pathways contain a smaller number of intermediate transition states than the previously publish lowestbarrier pathway. In case of the 102atom system Minima Hopping guided path search found a previously unknown and energetically lowlying funnel.

General coalescence conditions for the exact wave functions. II. Higherorder relations for manyparticle systems
View Description Hide DescriptionWe derived the necessary conditions that must be satisfied by the nonrelativistic timeindependent exact wave functions for manyparticle systems at a twoparticle coalescence (or cusp) point. Some simple conditions are known to be Kato's cusp condition (CC) and Rassolov and Chipman's CC. In a previous study, we derived an infinite number of necessary conditions that twoparticle wave functions must satisfy at a coalescence point. In the present study, we extend these conditions to manyparticle systems. They are called general coalescence conditions (GCCs), and Kato's CC and Rassolov and Chipman's CC are included as special conditions. GCCs can be applied not only to Coulombic systems but also to any system in which the interaction between two particles is represented in a power series of interparticle distances. We confirmed the correctness of our derivation of the GCCs by applying the exact wave function of a harmonium in electronelectron and electronnucleus coalescence situations. In addition, we applied the free complement (FC) wave functions of a helium atom to the GCCs to examine the accuracy of the FC wave function in the context of a coalescence situation.

Collecting highorder interactions in an effective pairwise intermolecular potential using the hydrated ion concept: The hydration of Cf^{3+}
View Description Hide DescriptionThis work proposes a new methodology to build interaction potentials between a highly charged metal cation and water molecules. These potentials, which can be used in classical computer simulations, have been fitted to reproduce quantum mechanical interaction energies (MP2 and BP86) for a wide range of [M(H2O) n ]^{ m+}(H2O)ℓ clusters (n going from 6 to 10 and ℓ from 0 to 18). A flexible and polarizable water shell model (Mobile Charge Density of Harmonic Oscillator) has been coupled to the cationwater potential. The simultaneous consideration of polyhydrated clusters and the polarizability of the interacting particles allows the inclusion of the most important manybody effects in the new polarizable potential. Applications have been centered on the californium, Cf(III) the heaviest actinoid experimentally studied in solution. Two different strategies to select a set of about 2000 structures which are used for the potential building were checked. Monte Carlo simulations of Cf(III)+500 H2O for three of the intermolecular potentials predict an aquaion structure with coordination number close to 8 and average in the range 2.43–2.48 Å, whereas the fourth one is closer to 9 with = 2.54 Å. Simulated EXAFS spectra derived from the structural Monte Carlo distribution compares fairly well with the available experimental spectrum for the simulations bearing 8 water molecules. An angular distribution similar to that of a square antiprism is found for the octacoordination.

The density matrix functional approach to electron correlation: Dynamic and nondynamic correlation along the full dissociation coordinate
View Description Hide DescriptionFor chemistry an accurate description of bond weakening and breaking is vital. The great advantage of density matrix functionals, as opposed to density functionals, is their ability to describe such processes since they naturally cover both nondynamical and dynamical correlation. This is obvious in the LöwdinShull functional, the exact natural orbital functional for twoelectron systems. We present in this paper extensions of this functional for the breaking of a single electron pair bond in Nelectron molecules, using LiH, BeH^{+}, and Li2 molecules as prototypes. Attention is given to the proper formulation of the functional in terms of not just J and K integrals but also the twoelectron L integrals (K integrals with a different distribution of the complex conjugation of the orbitals), which is crucial for the calculation of response functions. Accurate energy curves are obtained with extended LöwdinShull functionals along the complete dissociation coordinate using full CI calculations as benchmark.

Improved initial guess for minimum energy path calculations
View Description Hide DescriptionA method is presented for generating a good initial guess of a transition path between given initial and final states of a system without evaluation of the energy. An objective function surface is constructed using an interpolation of pairwise distances at each discretization point along the path and the nudged elastic band method then used to find an optimal path on this image dependent pair potential (IDPP) surface. This provides an initial path for the more computationally intensive calculations of a minimum energy path on an energy surface obtained, for example, by ab initio or density functional theory. The optimal path on the IDPP surface is significantly closer to a minimum energy path than a linear interpolation of the Cartesian coordinates and, therefore, reduces the number of iterations needed to reach convergence and averts divergence in the electronic structure calculations when atoms are brought too close to each other in the initial path. The method is illustrated with three examples: (1) rotation of a methyl group in an ethane molecule, (2) an exchange of atoms in an island on a crystal surface, and (3) an exchange of two Siatoms in amorphous silicon. In all three cases, the computational effort in finding the minimum energy path with DFT was reduced by a factor ranging from 50% to an order of magnitude by using an IDPP path as the initial path. The time required for parallel computations was reduced even more because of load imbalance when linear interpolation of Cartesian coordinates was used.

Energy relaxation and separation of a hot electronhole pair in organic aggregates from a timedependent wavepacket diffusion method
View Description Hide DescriptionThe timedependent wavepacket diffusive method [X. Zhong and Y. Zhao, J. Chem. Phys.138, 014111 (2013)] is extended to investigate the energy relaxation and separation of a hot electronhole pair in organic aggregates with incorporation of Coulomb interaction and electronphonon coupling. The pair initial condition generated by laser pulse is represented by a Gaussian wavepacket with a central momentum. The results reveal that the hot electron energy relaxation is very well described by two rate processes with the fast rate much larger than the slow one, consistent with experimental observations, and an efficient electronhole separation is accomplished accompanying the fast energy relaxation. Furthermore, although the extra energy indeed helps the separation by overcoming the Coulomb interaction, the width of initial wavepacket is much sensitive to the separation efficiency and the narrower wavepacket generates the more separated charges. This behavior may be useful to understand the experimental controversy of the hot carrier effect on charge separation.

Fluctuationinduced transport of two coupled particles: Effect of the interparticle interaction
View Description Hide DescriptionWe consider a system of two coupled particles fluctuating between two states, with different interparticle interaction potentials and particle friction coefficients. An external action drives the interstate transitions that induces reciprocating motion along the internal coordinate x (the interparticle distance). The system moves unidirectionally due to rectification of the internal motion by asymmetric friction fluctuations and thus operates as a dimeric motor that converts input energy into net movement. We focus on how the law of interaction between the particles affects the dimer transport and, in particular, the role of thermal noise in the motion inducing mechanism. It is argued that if the interaction potential behaves at large distances as x ^{α}, depending on the value of the exponent α, the thermal noise plays a constructive (α > 2), neutral (α = 2), or destructive (α < 2) role. In the case of α = 1, corresponding piecewise linear potential profiles, an exact solution is obtained and discussed in detail.

Orderparameteraided temperatureaccelerated sampling for the exploration of crystal polymorphism and solidliquid phase transitions
View Description Hide DescriptionThe problem of predicting polymorphism in atomic and molecular crystals constitutes a significant challenge both experimentally and theoretically. From the theoretical viewpoint, polymorphism prediction falls into the general class of problems characterized by an underlying rough energy landscape, and consequently, free energy based enhanced sampling approaches can be brought to bear on the problem. In this paper, we build on a scheme previously introduced by two of the authors in which the lengths and angles of the supercell are targeted for enhanced sampling via temperature accelerated adiabatic free energy dynamics[T. Q. Yu and M. E. Tuckerman, Phys. Rev. Lett.107, 015701 (2011)]. Here, that framework is expanded to include general order parameters that distinguish different crystalline arrangements as target collective variables for enhanced sampling. The resulting free energy surface, being of quite high dimension, is nontrivial to reconstruct, and we discuss one particular strategy for performing the free energy analysis. The method is applied to the study of polymorphism in xenon crystals at high pressure and temperature using the Steinhardt order parameters without and with the supercell included in the set of collective variables. The expected fcc and bcc structures are obtained, and when the supercell parameters are included as collective variables, we also find several new structures, including fcc states with hcp stacking faults. We also apply the new method to the solidliquid phase transition in copper at 1300 K using the same Steinhardt order parameters. Our method is able to melt and refreeze the system repeatedly, and the free energy profile can be obtained with high efficiency.

Molecular dynamics saddle search adaptive kinetic Monte Carlo
View Description Hide DescriptionA method for accelerating molecular dynamics simulations in rare event systems is described. From each new state visited, high temperature molecular dynamics trajectories are used to discover the set of escape mechanisms and rates. This event table is provided to the adaptive kinetic Monte Carlo algorithm to model the evolution of the system from state to state. Importantly, an estimator for the completeness of the calculated rate table in each state is derived. The method is applied to three model systems: adatom diffusion on Al(100); island diffusion on Pt(111); and vacancy cluster ripening in bulk Fe. Connections to the closely related temperature accelerated dynamics method of Voter and coworkers is discussed.

Adiabatic state preparation study of methylene
View Description Hide DescriptionQuantum computers attract much attention as they promise to outperform their classical counterparts in solving certain type of problems. One of them with practical applications in quantum chemistry is simulation of complex quantum systems. An essential ingredient of efficient quantum simulation algorithms are initial guesses of the exact wave functions with high enough fidelity. As was proposed in AspuruGuzik et al. [Science309, 1704 (2005)], the exact ground states can in principle be prepared by the adiabatic state preparation method. Here, we apply this approach to preparation of the lowest lying multireference singlet electronic state of methylene and numerically investigate preparation of this state at different molecular geometries. We then propose modifications that lead to speeding up the preparation process. Finally, we decompose the minimal adiabatic state preparation employing the direct mapping in terms of twoqubit interactions.

Block diagonalization of the equationofmotion coupled cluster effective Hamiltonian: Treatment of diabatic potential constants and triple excitations
View Description Hide DescriptionWe present a diabatization method applicable to spectroscopic studies based on EquationofMotion Coupled Cluster (EOMCC) energies and biorthogonal wavefunctions that uses the Block Diagonalization (BD) approaches of Cederbaum et al.[L. S. Cederbaum, J. Schirmer, and H. D. Meyer, J. Phys. A: Math. Gen.22, 2427 (1989)] and Domcke et al.[W. Domcke and C. Woywod, Chem. Phys. Lett.216, 362 (1993); W. Domcke, C. Woywod, and M. Stengle, Chem. Phys. Lett.226, 257 (1994)]. The method gives excellent agreement with coupling constants calculated using the analytic gradient approach of Ichino et al.[T. Ichino, J. Gauss, and J. F. Stanton, J. Chem. Phys.130, 174105 (2009)]. While the BD method is a finite difference approach, it can be applied at any geometry, can generate (pointwise) diabatic potential energy surfaces, and can be used with EOM wavefunctions that include triple (or higher) excitations. The method is applied to several model systems and its sensitivity to orbital choice, excitation space, and projection space is explored.

Seniority zero pair coupled cluster doubles theory
View Description Hide DescriptionCoupled cluster theory with single and double excitations accurately describes weak electron correlation but is known to fail in cases of strong static correlation. Fascinatingly, however, pair coupled cluster doubles (pCCD), a simplified version of the theory limited to pair excitations that preserve the seniority of the reference determinant (i.e., the number of unpaired electrons), has mean field computational cost and is an excellent approximation to the full configuration interaction (FCI) of the paired space provided that the orbital basis defining the pairing scheme is adequately optimized. In previous work, we have shown that optimization of the pairing scheme in the seniority zero FCI leads to a very accurate description of static correlation. The same conclusion extends to pCCD if the orbitals are optimized to make the pCCD energy stationary. We here demonstrate these results with numerous examples. We also explore the contributions of different seniority sectors to the coupled cluster doubles (CCD) correlation energy using different orbital bases. We consider both HartreeFock and Brueckner orbitals, and the role of orbital localization. We show how one can pair the orbitals so that the role of the Brueckner orbitals at the CCD level is retained at the pCCD level. Moreover, we explore ways of extending CCD to accurately describe strongly correlated systems.

Projected senioritytwo orbital optimization of the antisymmetric product of onereference orbital geminal
View Description Hide DescriptionWe present a new, nonvariational orbitaloptimization scheme for the antisymmetric product of onereference orbital geminal wave function. Our approach is motivated by the observation that an orbitaloptimized seniorityzero configuration interaction (CI) expansion yields similar results to an orbitaloptimized seniorityzeroplustwo CI expansion [L. Bytautas, T. M. Henderson, C. A. JimenezHoyos, J. K. Ellis, and G. E. Scuseria, J. Chem. Phys.135, 044119 (2011)]. A numerical analysis is performed for the C2 and LiF molecules, for the CH2 singlet diradical as well as for the symmetric stretching of hypothetical (linear) hydrogen chains. For these test cases, the proposed orbitaloptimization protocol yields similar results to its variational orbital optimization counterpart, but prevents symmetrybreaking of molecular orbitals in most cases.

A Brownian dynamics algorithm for colloids in curved manifolds
View Description Hide DescriptionThe manyparticle Langevin equation, written in local coordinates, is used to derive a Brownian dynamics simulation algorithm to study the dynamics of colloids moving on curved manifolds. The predictions of the resulting algorithm for the particular case of free particles diffusing along a circle and on a sphere are tested against analytical results, as well as with simulation data obtained by means of the standard Brownian dynamics algorithm developed by Ermak and McCammon [J. Chem. Phys.69, 1352 (1978)] using explicitly a confining external field. The latter method allows constraining the particles to move in regions very tightly, emulating the diffusion on the manifold. Additionally, the proposed algorithm is applied to strong correlated systems, namely, paramagnetic colloids along a circle and soft colloids on a sphere, to illustrate its applicability to systems made up of interacting particles.

When do we need to account for the geometric phase in excited state dynamics?
View Description Hide DescriptionWe investigate the role of the geometric phase (GP) in an internal conversion process when the system changes its electronic state by passing through a conical intersection (CI). Local analysis of a twodimensional linear vibronic coupling (LVC) model Hamiltonian near the CI shows that the role of the GP is twofold. First, it compensates for a repulsion created by the socalled diagonal Born–Oppenheimer correction. Second, the GP enhances the nonadiabatic transition probability for a wavepacket part that experiences a central collision with the CI. To assess the significance of both GP contributions we propose two indicators that can be computed from parameters of electronic surfaces and initial conditions. To generalize our analysis to Ndimensional systems we introduce a reduction of a general Ndimensional LVC model to an effective 2D LVC model using a mode transformation that preserves shorttime dynamics of the original Ndimensional model. Using examples of the bis(methylene) adamantyl and butatriene cations, and the pyrazine molecule we have demonstrated that their effective 2D models reproduce the shorttime dynamics of the corresponding full dimensional models, and the introduced indicators are very reliable in assessing GP effects.

Improved parameterization of interatomic potentials for rare gas dimers with densitybased energy decomposition analysis
View Description Hide DescriptionWe examine interatomic interactions for rare gas dimers using the densitybased energy decomposition analysis (DEDA) in conjunction with computational results from CCSD(T) at the complete basis set (CBS) limit. The unique DEDA capability of separating frozen density interactions from density relaxation contributions is employed to yield clean interaction components, and the results are found to be consistent with the typical physical picture that density relaxations play a very minimal role in rare gas interactions. Equipped with each interaction component as reference, we develop a new threeterm molecular mechanical force field to describe rare gas dimers: a smeared charge multipole model for electrostatics with charge penetration effects, a B3LYPD3 dispersion term for asymptotically correct longrange attractions that is screened at shortrange, and a BornMayer exponential function for the repulsion. The resulted force field not only reproduces rare gas interaction energies calculated at the CCSD(T)/CBS level, but also yields each interaction component (electrostatic or van der Waals) which agrees very well with its corresponding reference value.