Volume 138, Issue 5, 07 February 2013
Index of content:

A multiscale model for a colloidpolymer mixture is developed. The colloids are described as point particles interacting with each other and with the polymers with strongly repulsive potentials, while polymers interact with each other with a softer potential. The fluid in the suspension is taken into account by the multiparticle collision dynamics method (MPC). Considering a slit geometry where the suspension is confined between parallel repulsive walls, different possibilities for the hydrodynamic boundary conditions (b.c.) at the walls (slip versus stick) are treated. Quenching experiments are considered, where the system volume is suddenly reduced (keeping the density of the solvent fluid constant, while the colloid and polymer particle numbers are kept constant) and thus an initially homogeneous system is quenched deeply into the miscibility gap, where it is unstable. For various relative concentrations of colloids and polymers, the time evolution of the growing colloidrich and polymerrich domains are studied by molecular dynamics simulation, taking hydrodynamic effects mediated by the solvent into account via MPC. It is found that the domain size ℓ_{ d }(t) grows with time t as ℓ_{ d }(t) ∝ t ^{1/3} for stick and (at late stages) as ℓ_{ d }(t) ∝ t ^{2/3} for slip b.c., while breakup of percolating structures can cause a transient “arrest” of growth. While these findings apply for films that are 5–10 colloid diameters wide, for ultrathin films (1.5 colloid diameters wide) a regime with ℓ_{ d }(t) ∝ t ^{1/2} is also identified for rather shallow quenches.
 ARTICLES


Theoretical Methods and Algorithms

An accurate and linearscaling method for calculating chargetransfer excitation energies and diabatic couplings
View Description Hide DescriptionQuantum–mechanical methods that are both computationally fast and accurate are not yet available for electronic excitations having charge transfer character. In this work, we present a significant step forward towards this goal for those charge transfer excitations that take place between noncovalently bound molecules. In particular, we present a method that scales linearly with the number of noncovalently bound molecules in the system and is based on a twopronged approach: The molecular electronic structure of brokensymmetry chargelocalized states is obtained with the frozen density embedding formulation of subsystem densityfunctional theory; subsequently, in a postSCF calculation, the fullelectron Hamiltonian and overlap matrix elements among the chargelocalized states are evaluated with an algorithm which takes full advantage of the subsystem DFT density partitioning technique. The method is benchmarked against coupledcluster calculations and achieves chemical accuracy for the systems considered for intermolecular separations ranging from hydrogenbond distances to tens of Ångstroms. Numerical examples are provided for molecular clusters comprised of up to 56 noncovalently bound molecules.

Coupling of kinetic Monte Carlo simulations of surface reactions to transport in a fluid for heterogeneous catalytic reactor modeling
View Description Hide DescriptionWe have developed a method to couple kinetic Monte Carlo simulations of surface reactions at a molecular scale to transport equations at a macroscopic scale. This method is applicable to steady state reactors. We use a finite difference upwinding scheme and a gaptooth scheme to efficiently use a limited amount of kinetic Monte Carlo simulations. In general the stochastic kinetic Monte Carlo results do not obey mass conservation so that unphysical accumulation of mass could occur in the reactor. We have developed a method to perform mass balance corrections that is based on a stoichiometry matrix and a leastsquares problem that is reduced to a nonsingular set of linear equations that is applicable to any surface catalyzed reaction. The implementation of these methods is validated by comparing numerical results of a reactor simulation with a unimolecular reaction to an analytical solution. Furthermore, the method is applied to two reaction mechanisms. The first is the ZGB model for CO oxidation in which inevitable poisoning of the catalyst limits the performance of the reactor. The second is a model for the oxidation of NO on a Pt(111) surface, which becomes active due to lateral interaction at high coverages of oxygen. This reaction model is based on ab initio density functional theory calculations from literature.

Manybody dispersion interactions from the exchangehole dipole moment model
View Description Hide DescriptionIn this article, we present the extension of the exchangehole dipole moment model (XDM) of dispersion interactions to the calculation of twobody and threebody dispersion energy terms to any order, 2^{ l }pole oscillator strengths, and polarizabilities. By using the newlyformulated coefficients, we study the relative importance of the higherorder twobody and the leading nonadditive threebody (tripledipole) interactions in gasphase as well as in condensed systems. We show that the twobody terms up to R ^{−10}, but not the terms of higherorder, are essential in the correct description of the dispersion energy, while there are a number of difficulties related to the choice of the damping function, which precludes the use threebody tripledipole contributions in XDM. We conclude that further study is required before the threebody term can be used in production XDM densityfunctional calculations and point out the salient problems regarding its use.

Signatures of discrete breathers in coherent state quantum dynamics
View Description Hide DescriptionIn classical mechanics, discrete breathers (DBs) – a spatial timeperiodic localization of energy – are predicted in a large variety of nonlinear systems. Motivated by a conceptual bridging of the DB phenomena in classical and quantum mechanical representations, we study their signatures in the dynamics of a quantum equivalent of a classical mechanical point in phase space – a coherent state. In contrast to the classical point that exhibits either delocalized or localized motion, the coherent state shows signatures of both localized and delocalized behavior. The transition from normal to local modes have different characteristics in quantum and classical perspectives. Here, we get an insight into the connection between classical and quantum perspectives by analyzing the decomposition of the coherent state into system's eigenstates, and analyzing the spacial distribution of the wavefunction density within these eigenstates. We find that the delocalized and localized eigenvalue components of the coherent state are separated by a mixed region, where both kinds of behavior can be observed. Further analysis leads to the following observations. Considered as a function of coupling, energy eigenstates go through avoided crossings between tunneling and nontunneling modes. The dominance of tunneling modes in the high nonlinearity region is compromised by the appearance of new types of modes – high order tunneling modes – that are similar to the tunneling modes but have attributes of nontunneling modes. Certain types of excitations preferentially excite higher order tunneling modes, allowing one to study their properties. Since autocorrelation functions decrease quickly in highly nonlinear systems, shorttime dynamics are sufficient for modeling quantum DBs. This work provides a foundation for implementing modern semiclassical methods to model quantum DBs, bridging classical and quantum mechanical signatures of DBs, and understanding spectroscopic experiments that involve a coherent state.

A quadratically convergent VBSCF method
View Description Hide DescriptionA quadratically convergent valence bond selfconsistent field method is described where the simultaneous optimisation of orbitals and the coefficients of the configurations (VB structures) is based on a NewtonRaphson scheme. The applicability of the method is demonstrated in actual calculations. The convergence and efficiency are compared with the SuperCI method. A necessary condition to achieve convergence in the NewtonRaphson method is that the Hessian is positive definite. When this is not the case, a combination of the SuperCI and NewtonRaphson methods is shown to be an optimal choice instead of shifting the eigenvalues of the Hessian to make it positive definite. In the combined method, the first few iterations are performed with the SuperCI method and then the NewtonRaphson scheme is switched on based on an internal indicator. This approach is found computationally a more economical choice than using either the NewtonRaphson or SuperCI method alone to perform a full optimisation of the nonorthogonal orbitals.

Reaction coordinates, onedimensional Smoluchowski equations, and a test for dynamical selfconsistency
View Description Hide DescriptionWe propose a method for identifying accurate reaction coordinates among a set of trial coordinates. The method applies to special cases where motion along the reaction coordinate follows a onedimensional Smoluchowski equation. In these cases the reaction coordinate can predict its own shorttime dynamical evolution, i.e., the dynamics projected from multiple dimensions onto the reaction coordinate depend only on the reaction coordinate itself. To test whether this property holds, we project an ensemble of short trajectory swarms onto trial coordinates and compare projections of individual swarms to projections of the ensemble of swarms. The comparison, quantified by the KullbackLeibler divergence, is numerically performed for each isosurface of each trial coordinate. The ensemble of short dynamical trajectories is generated only once by sampling along an initial order parameter. The initial order parameter should separate the reactants and products with a free energy barrier, and distributions on isosurfaces of the initial parameter should be unimodal. The method is illustrated for three model free energy landscapes with anisotropic diffusion. Where exact coordinates can be obtained from KramersLangerBerezhkovskiiSzabo theory, results from the new method agree with the exact results. We also examine characteristics of systems where the proposed method fails. We show how dynamical selfconsistency is related (through the ChapmanKolmogorov equation) to the earlier isocommittor criterion, which is based on longer paths.

The Schrödinger equation with friction from the quantum trajectory perspective
View Description Hide DescriptionSimilarity of equations of motion for the classical and quantum trajectories is used to introduce a friction term dependent on the wavefunction phase into the timedependent Schrödinger equation. The term describes irreversible energy loss by the quantum system. The force of friction is proportional to the velocity of a quantum trajectory. The resulting Schrödinger equation is nonlinear, conserves wavefunction normalization, and evolves an arbitrary wavefunction into the ground state of the system (of appropriate symmetry if applicable). Decrease in energy is proportional to the average kinetic energy of the quantum trajectory ensemble. Dynamics in the high friction regime is suitable for simple models of reactions proceeding with energy transfer from the system to the environment. Examples of dynamics are given for single and symmetric and asymmetric double well potentials.

Structural fluctuation of protein in water around its native state: A new statistical mechanics formulation
View Description Hide DescriptionA new statistical mechanics formulation of characterizing the structural fluctuation of protein correlated with that of water is presented based on the generalized Langevin equation and the 3Dreference interaction site model (RISM)/RISM theory of molecular liquids. The displacement vector of atom positions, and their conjugated momentum, are chosen for the dynamic variables for protein, while the density fields of atoms and their momentum fields are chosen for water. Projection of other degrees of freedom onto those dynamic variables using the standard projection operator method produces essentially two equations, which describe the time evolution of fluctuation concerning the density field of solvent and the conformation of protein around an equilibrium state, which are coupled with each other. The equation concerning the protein dynamics is formally akin to that of the coupled Langevin oscillators, and is a generalization of the latter, to atomic level. The most intriguing feature of the new equation is that it contains the variancecovariance matrix as the “Hessian” term describing the “force” restoring an equilibrium conformation, which is the second moment of the fluctuation of atom positions. The “Hessian” matrix is naturally identified as the second derivative of the free energy surface around the equilibrium. A method to evaluate the Hessian matrix based on the 3DRISM/RISM theory is proposed. Proposed also is an application of the present formulation to the molecular recognition, in which the conformational fluctuation of protein around its native state becomes an important factor as exemplified by so called “induced fitting.”

The orbitalspecific virtual local triples correction: OSVL(T)
View Description Hide DescriptionA local method based on orbital specific virtuals (OSVs) for calculating the perturbative triples correction in local coupled cluster calculations is presented. In contrast to the previous approach based on projected atomic orbitals (PAOs), described by Schütz [J. Chem. Phys.113, 9986 (Year: 2000)]10.1063/1.1323265, the new scheme works without any ad hoc truncations of the virtual space to domains. A single threshold defines the pair and triple specific virtual spaces completely and automatically. It is demonstrated that the computational cost of the method scales linearly with molecular size. Employing the recommended threshold a similar fraction of the correlation energy is recovered as with the original PAO method at a somewhat lower cost. A benchmark for 52 reactions demonstrates that for reaction energies the intrinsic accuracy of the coupled cluster with singles and doubles excitations and a perturbative treatment of triples excitations method can be reached by OSVlocal coupled cluster theory with singles and doubles and perturbative triples, provided a MP2 correction is applied that accounts for basis set incompleteness errors as well as for remaining domain errors. As an application example the interaction energies of the guaninecytosine dimers in the WatsonCrick and stacked arrangements are investigated at the level of local coupled cluster theory with singles and doubles and perturbative triples. Based on these calculations we propose new completebasissetlimit estimates for these interaction energies at this level of theory.

Going beyond the frozen core approximation: Development of coordinatedependent pseudopotentials and application to
View Description Hide DescriptionMixed quantum/classical (MQC) simulations treat the majority of a system classically and reserve quantum mechanics only for a few degrees of freedom that actively participate in the chemical process(es) of interest. In MQC calculations, the quantum and classical degrees of freedom are coupled together using pseudopotentials. Although most pseudopotentials are developed empirically, there are methods for deriving pseudopotentials using the results of quantum chemistry calculations, which guarantee that the explicitlytreated valence electron wave functions remain orthogonal to the implicitlytreated core electron orbitals. Whether empirical or analytically derived in nature, to date all such pseudopotentials have been subject to the frozen core approximation (FCA) that ignores how changes in the nuclear coordinates alter the core orbitals, which in turn affects the wave function of the valence electrons. In this paper, we present a way to go beyond the FCA by developing pseudopotentials that respond to these changes. In other words, we show how to derive an analytic expression for a pseudopotential that is an explicit function of nuclear coordinates, thus accounting for the polarization effects experienced by atomic cores in different chemical environments. We then use this formalism to develop a coordinatedependent pseudopotential for the bonding electron of the sodium dimer cation molecule and we show how the analytic representation of this potential can be used in oneelectron MQC simulations that provide the accuracy of a fully quantum mechanical HartreeFock (HF) calculation at all internuclear separations. We also show that oneelectron MQC simulations of using our coordinatedependent pseudopotential provide a significant advantage in accuracy compared to frozen core potentials with no additional computational expense. This is because use of a frozen core potential produces a charge density for the bonding electron of that is too localized on the molecule, leading to significant overbinding of the valence electron. This means that FCA calculations are subject to inaccuracies of order ∼10% in the calculated bond length and vibrational frequency of the molecule relative to a full HF calculation; these errors are fully corrected by using our coordinatedependent pseudopotential. Overall, our findings indicate that even for molecules like , which have a simple electronic structure that might be expected to be welltreated within the FCA, the importance of including the effects of the changing core molecular orbitals on the bonding electrons cannot be overlooked.

Modeling molecular response in uniform and nonuniform electric fields
View Description Hide DescriptionThe response of a molecule to an electric field E, often a model of environment, can be expressed in terms of a sum of power series expansions. We investigate the accuracy and limits of applicability of this expression using one, two, and threedimensional models of the hydrogenbonded complex, ClH:NH_{3}. Energetic, structural, and vibrational spectroscopic characteristics are determined at first and secondorder in E and ∇E and compared with ab initio values for a range of uniform and nonuniform electric fields chosen to simulate molecular environments. It is found that even at field strengths large enough to cause dramatic structural change in the complex, energetic, structural, and vibrational spectroscopic characteristics are accurately calculated using only terms linear in E and ∇E. These results suggest that knowledge of the zerofield molecular potential energy, dipole, and quadrupole moment surfaces may be sufficient to accurately model the interaction of a molecule with a wide range of chemical environments.

Generalized Born forces: Surface integral formulation
View Description Hide DescriptionGeneralized Born (GB) models offer a convenient alternative to PoissonBoltzmann based models. In the last decade, the GB radii computed based on the exact results obtained for a charge embedded in a conducting sphere have proven to be accurate also for the complex molecular shapes of proteins. The surface integral formulation of the theory has been much less explored than the volume integral formulation. In this work, we provide the exact equations for the GB solvation forces in the surface integral formulation, which are nontrivial due to the nonnegligible dependence of GB radii on atomic positions and due to the discontinuity in the derivative of the solvent accessible surface point positions with respect to atomic positions. The equations derived here provide a useful reference for developing faster approximations.

NMR chemical shift as analytical derivative of the Helmholtz free energy
View Description Hide DescriptionWe present a theory for the temperaturedependent nuclear magnetic shielding tensor of molecules with arbitrary electronic structure. The theory is a generalization of Ramsey's theory for closedshell molecules. The shielding tensor is defined as a second derivative of the Helmholtz free energy of the electron system in equilibrium with the applied magnetic field and the nuclear magnetic moments. This derivative is analytically evaluated and expressed as a sum over states formula. Special consideration is given to a system with an isolated degenerate ground state for which the size of the degeneracy and the composition of the wave functions are arbitrary. In this case, the paramagnetic part of the shielding tensor is expressed in terms of the g and A tensors of the electron paramagnetic resonance spin Hamiltonian of the degenerate state. As an illustration of the proposed theory, we provide an explicit formula for the paramagnetic shift of the central lanthanide ion in endofullerenes Ln@C_{60}, with Ln = Ce^{3+}, Nd^{3+}, Sm^{3+}, Dy^{3+}, Er^{3+}, and Yb^{3+}, where the ground state can be a strongly spinorbit coupled icosahedral sextet for which the paramagnetic shift cannot be described by previous theories.

Development of polarontransformed explicitly correlated full configuration interaction method for investigation of quantumconfined Stark effect in GaAs quantum dots
View Description Hide DescriptionThe effect of external electric field on electronhole (eh) correlation in gallium arsenide quantum dots is investigated. The electronhole Schrodinger equation in the presence of an external electric field is solved using explicitly correlated full configuration interaction method and accurate exciton binding energy and electronhole recombination probability are obtained. The effect of the electric field was included in the 1particle single component basis functions by performing variational polaron transformation. The quality of the wavefunction at small interparticle distances was improved by using Gaussiantype geminal function that depended explicitly on the electronhole separation distance. The parameters of the explicitly correlated function were determined variationally at each field strength. The scaling of total exciton energy, exciton binding energy, and electronhole recombination probability with respect to the strength of the electric field was investigated. It was found that a 500 kV/cm change in field strength reduces the binding energy and recombination probability by a factor of 2.6 and 166, respectively. The results show that the ehrecombination probability is affected much more strongly by the electric field than the exciton binding energy. Analysis using the polarontransformed basis indicates that the exciton binding should asymptotically vanish in the limit of large field strength.

A relative entropy rate method for path space sensitivity analysis of stationary complex stochastic dynamics
View Description Hide DescriptionWe propose a new sensitivity analysis methodology for complex stochastic dynamics based on the relative entropy rate. The method becomes computationally feasible at the stationary regime of the process and involves the calculation of suitable observables in path space for the relative entropy rate and the corresponding Fisher information matrix. The stationary regime is crucial for stochastic dynamics and here allows us to address the sensitivity analysis of complex systems, including examples of processes with complex landscapes that exhibit metastability, nonreversible systems from a statistical mechanics perspective, and highdimensional, spatially distributed models. All these systems exhibit, typically nonGaussian stationary probability distributions, while in the case of highdimensionality, histograms are impossible to construct directly. Our proposed methods bypass these challenges relying on the direct Monte Carlo simulation of rigorously derived observables for the relative entropy rate and Fisher information in path space rather than on the stationary probability distribution itself. We demonstrate the capabilities of the proposed methodology by focusing here on two classes of problems: (a) Langevin particle systems with either reversible (gradient) or nonreversible (nongradient) forcing, highlighting the ability of the method to carry out sensitivity analysis in nonequilibrium systems; and, (b) spatially extended kinetic Monte Carlo models, showing that the method can handle highdimensional problems.

Accelerated direct semiclassical molecular dynamics using a compact finite difference Hessian scheme
View Description Hide DescriptionThis paper shows how a compact finite difference Hessian approximation scheme can be proficiently implemented into semiclassical initial value representation molecular dynamics. Effects of the approximation on the monodromy matrix calculation are tested by propagating initial sampling distributions to determine power spectra for analytic potential energy surfaces and for “on the fly” carbon dioxide direct dynamics. With the approximation scheme the computational cost is significantly reduced, making ab initio direct semiclassical dynamics computationally more feasible and, at the same time, properly reproducing important quantum effects inherent in the monodromy matrix and the preexponential factor of the semiclassical propagator.

Polarizability effects in molecular dynamics simulations of the graphenewater interface
View Description Hide DescriptionThe importance of including the polarizability of both water and graphene in molecular dynamics simulations of the water/graphene system was quantified. A thin film of either rigid single point charge extended (SPC/E) water or polarizable simple 4site water model with Drude polarizability (SWM4_DP) water on nonpolarizable and polarizable graphene surfaces was simulated. The graphene surface was either maintained neutral or charged, positively and negatively. The results suggest that SPC/E and SWM4_DP water models yield very similar predictions for the water structural properties on neutral nonpolarizable graphene, although they yield slightly different dynamical properties of interfacial water on neutral nonpolarizable graphene. More pronounced were the differences obtained when graphene was modeled with a polarizable force field. In particular, the polarizability of graphene was found to enhance the number of interfacial SWM4_DP water molecules pointing one of their OH bonds towards the neutral surface. Despite this structural difference, the dynamical properties predicted for the interfacial SWM4_DP water were found to be independent on polarizability as long as the polarizability of a carbon atom is smaller than α = 0.878 Å. On charged graphene surfaces, the effect of polarizability of graphene on structural properties and some dynamical properties of SWM4_DP water is negligible because electrostatic forces due to surface charge dominate polarization forces, as expected. For all cases, our results suggest that the hydrogen bond network is insensitive to the polarizability of both water and graphene. Understanding how these effects will determine the accumulation of ions near neutral or charged graphene could have important implications for applications in the fields of energy storage and water desalination.

Quantum dynamical structure factor of liquid neon via a quasiclassical symmetrized method
View Description Hide DescriptionWe apply the phase integration method for quasiclassical quantum time correlation functions[M. Monteferrante, S. Bonella, and G. Ciccotti, Mol. Phys.109, 3015 (Year: 2011)10.1080/00268976.2011.619506] to compute the dynamic structure factor of liquid neon. So far the method had been tested only on model systems. By comparing our results for neon with experiments and previous calculations, we demonstrate that the scheme is accurate and efficient also for a realistic model of a condensed phase system showing quantum behavior.

A variational formulation of electrostatics in a medium with spatially varying dielectric permittivity
View Description Hide DescriptionIn biological and synthetic materials, many important processes involve charges that are present in a medium with spatially varying dielectric permittivity. To accurately understand the role of electrostatic interactions in such systems, it is important to take into account the spatial dependence of the permittivity of the medium. However, due to the ensuing theoretical and computational challenges, this inhomogeneous dielectric response of the medium is often ignored or excessively simplified. We develop a variational formulation of electrostatics to accurately investigate systems that exhibit this inhomogeneous dielectric response. Our formulation is based on a true energy functional of the polarization charge density. The defining characteristic of a true energy functional is that at its minimum it evaluates to the actual value of the energy; this is a feature not found in many commonly used electrostatic functionals. We explore in detail the charged systems that exhibit sharp discontinuous change in dielectric permittivity, and we show that for this case our functional reduces to a functional of only the surface polarization charge density. We apply this reduced functional to study model problems for which analytical solutions are well known. We demonstrate, in addition, that the functional has many properties that make it ideal for use in molecular dynamics simulations.

Combinedhyperbolicinversepowerrepresentation of potential energy surfaces: A preliminary assessment for and
View Description Hide DescriptionThe purpose is to fit an accurate smooth function of the manybody expansion type to a multidimensional large data set using a basisset type method. By adopting a combinedhyperbolicinversepowerrepresentation for the basis, the novel approach is tested in detail for the ground electronic state of trihydrogen and hydroperoxyl systems, assuming that their potential energy surfaces are singlesheeted representable. It is also shown that the method can be easily applicable to potential energy curves by considering as prototypes molecular oxygen and the hydroxyl radical.
