Volume 139, Issue 21, 07 December 2013

Collective variables (CVs) are lowdimensional representations of the state of a complex system, which help us rationalize molecular conformations and sample free energy landscapes with molecular dynamics simulations. Given their importance, there is need for systematic methods that effectively identify CVs for complex systems. In recent years, nonlinear manifold learning has shown its ability to automatically characterize molecular collective behavior. Unfortunately, these methods fail to provide a differentiable function mapping highdimensional configurations to their lowdimensional representation, as required in enhanced sampling methods. We introduce a methodology that, starting from an ensemble representative of molecular flexibility, builds smooth and nonlinear datadriven collective variables (SandCV) from the output of nonlinear manifold learning algorithms. We demonstrate the method with a standard benchmark molecule, alanine dipeptide, and show how it can be nonintrusively combined with offtheshelf enhanced sampling methods, here the adaptive biasing force method. We illustrate how enhanced sampling simulations with SandCV can explore regions that were poorly sampled in the original molecular ensemble. We further explore the transferability of SandCV from a simpler system, alanine dipeptide in vacuum, to a more complex system, alanine dipeptide in explicit water.
 COMMUNICATIONS


Communication: The correct interpretation of surface hopping trajectories: How to calculate electronic properties
View Description Hide DescriptionIn a recent paper, we presented a road map for how Tully's fewest switches surface hopping (FSSH) algorithm can be derived, under certain circumstances, from the mixed quantumclassical Liouville equation. In this communication, we now demonstrate how this new interpretation of surface hopping can yield significantly enhanced results for electronic properties in nonadiabatic calculations. Specifically, we calculate diabatic populations for the spinboson problem using FSSH trajectories. We show that, for some Hamiltonians, without changing the FSSH algorithm at all but rather simply reinterpreting the ensemble of surface hopping trajectories, we recover excellent results and remove any and all ambiguity about the initial condition problem.

Communication: A reducedspace algorithm for the solution of the complex linear response equations used in coupled cluster damped response theory
View Description Hide DescriptionWe present a reducedspace algorithm for solving the complex (damped) linear response equations required to compute the complex linear response function for the hierarchy of methods: coupled cluster singles, coupled cluster singles and iterative approximate doubles, and coupled cluster singles and doubles. The solver is the keystone element for the development of damped coupled cluster response methods for linear and nonlinear effects in resonant frequency regions.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

Modeling and enhanced sampling of molecular systems with smooth and nonlinear datadriven collective variables
View Description Hide DescriptionCollective variables (CVs) are lowdimensional representations of the state of a complex system, which help us rationalize molecular conformations and sample free energy landscapes with molecular dynamics simulations. Given their importance, there is need for systematic methods that effectively identify CVs for complex systems. In recent years, nonlinear manifold learning has shown its ability to automatically characterize molecular collective behavior. Unfortunately, these methods fail to provide a differentiable function mapping highdimensional configurations to their lowdimensional representation, as required in enhanced sampling methods. We introduce a methodology that, starting from an ensemble representative of molecular flexibility, builds smooth and nonlinear datadriven collective variables (SandCV) from the output of nonlinear manifold learning algorithms. We demonstrate the method with a standard benchmark molecule, alanine dipeptide, and show how it can be nonintrusively combined with offtheshelf enhanced sampling methods, here the adaptive biasing force method. We illustrate how enhanced sampling simulations with SandCV can explore regions that were poorly sampled in the original molecular ensemble. We further explore the transferability of SandCV from a simpler system, alanine dipeptide in vacuum, to a more complex system, alanine dipeptide in explicit water.

Extended Lagrangian BornOppenheimer molecular dynamics in the limit of vanishing selfconsistent field optimization
View Description Hide DescriptionWe present an efficient general approach to first principles molecular dynamics simulations based on extended Lagrangian BornOppenheimer molecular dynamics [A. M. N. Niklasson, Phys. Rev. Lett.100, 123004 (2008)] in the limit of vanishing selfconsistent field optimization. The reduction of the optimization requirement reduces the computational cost to a minimum, but without causing any significant loss of accuracy or longterm energy drift. The optimizationfree first principles molecular dynamics requires only one single diagonalization per time step, but is still able to provide trajectories at the same level of accuracy as “exact,” fully converged, BornOppenheimer molecular dynamics simulations. The optimizationfree limit of extended Lagrangian BornOppenheimer molecular dynamics therefore represents an ideal starting point for robust and efficient first principles quantum mechanical molecular dynamics simulations.

Linearscaling calculation of HartreeFock exchange energy with nonorthogonal generalised Wannier functions
View Description Hide DescriptionWe present a method for the calculation of fourcentre twoelectron repulsion integrals in terms of localised nonorthogonal generalised Wannier functions (NGWFs). Our method has been implemented in the ONETEP program and is used to compute the HartreeFock exchange energy component of HartreeFock and Density Functional Theory (DFT) calculations with hybrid exchangecorrelation functionals. As the NGWFs are optimised in situin terms of a systematically improvable basis set which is equivalent to plane waves, it is possible to achieve large basis set accuracy in routine calculations. The spatial localisation of the NGWFs allows us to exploit the exponential decay of the density matrix in systems with a band gap in order to compute the exchange energy with a computational effort that increases linearly with the number of atoms. We describe the implementation of this approach in the ONETEPprogram for linearscaling first principles quantum mechanical calculations. We present extensive numerical validation of all the steps in our method. Furthermore, we find excellent agreement in energies and structures for a wide variety of molecules when comparing with other codes. We use our method to perform calculations with the B3LYP exchangecorrelation functional for models of myoglobin systems bound with O2 and CO ligands and confirm that the same qualitative behaviour is obtained as when the same myoglobin models are studied with the DFT+U approach which is also available in ONETEP. Finally, we confirm the linearscaling capability of our method by performing calculations on polyethylene and polyacetylene chains of increasing length.

An improved fragmentbased quantum mechanical method for calculation of electrostatic solvation energy of proteins
View Description Hide DescriptionAn efficient approach that combines the electrostatically embedded generalized molecular fractionation with conjugate caps (EEGMFCC) method with conductorlike polarizable continuum model (CPCM), termed EEGMFCCCPCM, is developed for ab initio calculation of the electrostatic solvation energy of proteins. Compared with the previous MFCCCPCM study [Y. Mei, C. G. Ji, and J. Z. H. Zhang, J. Chem. Phys.125, 094906 (2006)], quantum mechanical (QM) calculation is applied to deal with shortrange nonneighboring interactions replacing the classical treatment. Numerical studies are carried out for proteins up to 3837 atoms at the HF/631G* level. As compared to standard full system CPCM calculations, EEGMFCCCPCM shows clear improvement over the MFCCCPCM method for both the total electrostatic solvation energy and its components (the polarized solutesolvent reaction field energy and wavefunction distortion energy of the solute). For large proteins with 1000–4000 atoms, where the standard full system ab initio CPCM calculations are not affordable, the EEGMFCCCPCM gives larger relative wavefunction distortion energies and weaker relative electrostatic solvation energies for proteins, as compared to the corresponding energies calculated by the DivideandConquer PoissonBoltzmann (D&CPB) method. Notwithstanding, a high correlation between EEGMFCCCPCM and D&CPB is observed. This study demonstrates that the linearscaling EEGMFCCCPCM approach is an accurate and also efficient method for the calculation of electrostatic solvation energy of proteins.

HamiltonJacobi method for molecular distribution function in a chemical oscillator
View Description Hide DescriptionUsing the HamiltonJacobi method, we solve chemical FokkerPlanck equations within the Gaussian approximation and obtain a simple and compact formula for a conditional probability distribution. The formula holds in general transient situations, and can be applied not only to a steady state but also to an oscillatory state. By analyzing the long time behavior of the solution in the oscillatory case, we obtain the phase diffusion constant along the periodic orbit and the steady distribution perpendicular to it. A simple method for numerical evaluation of these formulas are devised, and they are compared with Monte Carlo simulations in the case of Brusselator as an example. Some results are shown to be identical to previously obtained expressions.

Resummed thermodynamic perturbation theory for bond cooperativity in associating fluids
View Description Hide DescriptionWe develop a resummed thermodynamic perturbation theory for bond cooperativity in associating fluids by extension of Wertheim's multidensity formalism. We specifically consider the case of an associating hard sphere with two association sites and both pairwise and triplet contributions to the energy, such that the first bond in an associated cluster receives an energy −ɛ^{(1)} and each subsequent bond in the cluster receives an energy −ɛ^{(2)}. To test the theory we perform new Monte Carlo simulations for potentials of this type. Theory and simulation are found to be in excellent agreement. We show that decreasing the energetic benefit of hydrogen bonding can actually result in a decrease in internal energy in the fluid. We also predict that when ɛ^{(1)} = 0 and ɛ^{(2)} is nonzero there is a transition temperature where the system transitions from a fluid of monomers to a mixture of monomers and very long chains.

Can we derive Tully's surfacehopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence
View Description Hide DescriptionIn this article, we demonstrate that Tully's fewestswitches surface hopping (FSSH) algorithm approximately obeys the mixed quantumclassical Liouville equation (QCLE), provided that several conditions are satisfied – some major conditions, and some minor. The major conditions are: (1) nuclei must be moving quickly with large momenta; (2) there cannot be explicit recoherences or interference effects between nuclear wave packets; (3) forcebased decoherence must be added to the FSSH algorithm, and the trajectories can no longer rigorously be independent (though approximations for independent trajectories are possible). We furthermore expect that FSSH (with decoherence) will be most robust when nonadiabatic transitions in an adiabatic basis are dictated primarily by derivative couplings that are presumably localized to crossing regions, rather than by small but pervasive offdiagonal force matrix elements. In the end, our results emphasize the strengths of and possibilities for the FSSH algorithm when decoherence is included, while also demonstrating the limitations of the FSSH algorithm and its inherent inability to follow the QCLE exactly.

Time dependent quantum thermodynamics of a coupled quantum oscillator system in a small thermal environment
View Description Hide DescriptionSimulations are performed of a small quantum system interacting with a quantum environment. The system consists of various initial states of two harmonic oscillators coupled to give normal modes. The environment is “designed” by its level pattern to have a thermodynamic temperature. A random coupling causes the system and environment to become entangled in the course of time evolution. The approach to a Boltzmann distribution is observed, and effective fitted temperatures close to the designed temperature are obtained. All initial pure states of the system are driven to equilibrium at very similar rates, with quick loss of memory of the initial state. The time evolution of the von Neumann entropy is calculated as a measure of equilibration and of quantum coherence. It is pointed out using spatial density distribution plots that quantum interference is eliminated only with maximal entropy, which corresponds thermally to infinite temperature. Implications of our results for the notion of “classicalizing” behavior in the approach to thermal equilibrium are briefly considered.

Efficient basis sets for noncovalent interactions in XDMcorrected densityfunctional theory
View Description Hide DescriptionIn the development and application of dispersioncorrected densityfunctional theory, the effects of basis set incompleteness have been largely mitigated through the use of very large, nearlycomplete basis sets. However, the use of such large basis sets makes application of these methods inefficient for large systems. In this work, we examine a series of basis sets, including Poplestyle, correlationconsistent, and polarizationconsistent bases, for their ability to efficiently and accurately predict noncovalent interactions when used in conjunction with the exchangehole dipole moment (XDM) dispersion model. We find that the polarizationconsistent 2 (pc2) basis sets, and two modifications thereof with some diffuse functions removed, give performance of comparable quality to that obtained with augccpVTZ basis sets, while being roughly 12 to 23 times faster computationally. The behavior is explained, in part, by the role of diffuse functions in recovering small density changes in the intermolecular region. The general performance of the modified basis sets is tested by application of XDM to standard intermolecular benchmark sets at, and away from, equilibrium.

Selfconsistent continuum solvation (SCCS): The case of charged systems
View Description Hide DescriptionThe recently developed selfconsistent continuum solvation model (SCCS) [O. Andreussi, I. Dabo, and N. Marzari, J. Chem. Phys.136, 064102 (2012)] is applied here to charged species in aqueous solutions. Describing ions in solution represents a great challenge because of the large electrostatic interactions between the solute and the solvent. The SCCS model is tested over 106 monocharged species, both cations and anions, and we demonstrate its flexibility, notwithstanding its much reduced set of parameters, to describe charged species in solution. Remarkably low mean absolute errors are obtained with values of 2.27 and 5.54 kcal/mol for cations and anions, respectively. These results are comparable or better than the state of the art to describe solvation of charged species in water. Finally, differences of behavior between cations and anions are discussed.

Proximal distributions from angular correlations: A measure of the onset of coarsegraining
View Description Hide DescriptionIn this work we examine and extend the theory of proximal radial distribution functions for molecules in solution. We point out two formal extensions, the first of which generalizes the proximal distribution function hierarchy approach to the complete, angularly dependent molecular pair distribution function. Second, we generalize from the traditional righthanded solutesolvent proximal distribution functions to the lefthanded distributions. The resulting neighbor hierarchy convergence is shown to provide a measure of the coarsegraining of the internal solute sites with respect to the solvent. Simulation of the test case of a decaalanine peptide shows that this coarsegraining measure converges at a length scale of approximately 5 amino acids for the system considered.

A minimallyresolved immersed boundary model for reactiondiffusion problems
View Description Hide DescriptionWe develop an immersed boundary approach to modeling reactiondiffusion processes in dispersions of reactive spherical particles, from the diffusionlimited to the reactionlimited setting. We represent each reactive particle with a minimallyresolved “blob” using many fewer degrees of freedom per particle than standard discretization approaches. More complicated or more highly resolved particle shapes can be built out of a collection of reactive blobs. We demonstrate numerically that the blob model can provide an accurate representation at low to moderate packing densities of the reactive particles, at a cost not much larger than solving a Poisson equation in the same domain. Unlike multipole expansion methods, our method does not require analytically computed Green's functions, but rather, computes regularized discrete Green's functions on the fly by using a standard gridbased discretization of the Poisson equation. This allows for great flexibility in implementing different boundary conditions, coupling to fluid flow or thermal transport, and the inclusion of other effects such as temporal evolution and even nonlinearities. We develop multigridbased preconditioners for solving the linear systems that arise when using implicit temporal discretizations or studying steady states. In the diffusionlimited case the resulting linear system is a saddlepoint problem, the efficient solution of which remains a challenge for suspensions of many particles. We validate our method by comparing to published results on reactiondiffusion in ordered and disordered suspensions of reactive spheres.

The StokesEinstein relation at moderate Schmidt number
View Description Hide DescriptionThe StokesEinstein relation for the selfdiffusion coefficient of a spherical particle suspended in an incompressible fluid is an asymptotic result in the limit of large Schmidt number, that is, when momentum diffuses much faster than the particle. When the Schmidt number is moderate, which happens in most particle methods for hydrodynamics, deviations from the StokesEinstein prediction are expected. We study these corrections computationally using a recently developed minimally resolved method for coupling particles to an incompressible fluctuating fluid in both two and three dimensions. We find that for moderate Schmidt numbers the diffusion coefficient is reduced relative to the StokesEinstein prediction by an amount inversely proportional to the Schmidt number in both two and three dimensions. We find, however, that the Einstein formula is obeyed at all Schmidt numbers, consistent with linear response theory. The mismatch arises because thermal fluctuations affect the drag coefficient for a particle due to the nonlinear nature of the fluidparticle coupling. The numerical data are in good agreement with an approximate selfconsistent theory, which can be used to estimate finiteSchmidt number corrections in a variety of methods. Our results indicate that the corrections to the StokesEinstein formula come primarily from the fact that the particle itself diffuses together with the momentum. Our study separates effects coming from corrections to noslip hydrodynamics from those of finite separation of time scales, allowing for a better understanding of widely observed deviations from the StokesEinstein prediction in particle methods such as molecular dynamics.

Spinfree DiracCoulomb calculations augmented with a perturbative treatment of spinorbit effects at the HartreeFock level
View Description Hide DescriptionA perturbative approach to compute secondorder spinorbit (SO) corrections to a spinfree DiracCoulomb HartreeFock (SFDCHF) calculation is suggested. The proposed scheme treats the difference between the DC and SFDC Hamiltonian as perturbation and exploits analytic secondderivative techniques. In addition, a costeffective scheme for incorporating relativistic effects in highaccuracy calculations is suggested consisting of a SFDC coupledcluster treatment augmented by perturbative SO corrections obtained at the HF level. Benchmark calculations for the hydrogen halides HX, X = FAt as well as the coinagemetal fluorides CuF, AgF, and AuF demonstrate the accuracy of the proposed perturbative treatment of SO effects on energies and electrical properties in comparison with the more rigorous full DC treatment. Furthermore, we present, as an application of our scheme, results for the electrical properties of AuF and XeAuF.

Diffusion in narrow channels on curved manifolds
View Description Hide DescriptionIn this work, we derive a general effective diffusion coefficient to describe the twodimensional (2D) diffusion in a narrow and smoothly asymmetric channel of varying width, embedded on a curved surface, in the simple diffusion of noninteracting, pointlike particles under no external field. To this end, we extend the generalization of the Kalinay–Percus' projection method [J. Chem. Phys.122, 204701 (2005); Kalinay–Percus', Phys. Rev. E74, 041203 (2006)] for the asymmetric channels introduced in [L. Dagdug and I. Pineda, J. Chem. Phys.137, 024107 (2012)], to project the anisotropic twodimensional diffusion equation on a curved manifold, into an effective onedimensional generalized FickJacobs equation that is modified according to the curvature of the surface. For such purpose we construct the whole expansion, writing the marginal concentration as a perturbation series. The lowest order in the perturbation parameter, which corresponds to the FickJacobs equation, contains an additional term that accounts for the curvature of the surface. We explicitly obtain the firstorder correction for the invariant effective concentration, which is defined as the correct marginal concentration in one variable, and we obtain the first approximation to the effective diffusion coefficient analogous to Bradley's coefficient [Phys. Rev. E80, 061142 (2009)] as a function of the metric elements of the surface. In a straightforward manner, we study the perturbation series up to the nth order, and derive the full effective diffusion coefficient for twodimensional diffusion in a narrow asymmetric channel, with modifications according to the metric terms. This expression is given as , which is the main result of our work. Finally, we present two examples of symmetric surfaces, namely, the sphere and the cylinder, and we study certain specific channel configurations on these surfaces.
 Advanced Experimental Techniques

Observation of strongly forbidden solid effect dynamic nuclear polarization transitions via electronelectron double resonance detected NMR
View Description Hide DescriptionWe present electron paramagnetic resonance experiments for which solid effect dynamic nuclear polarization transitions were observed indirectly via polarization loss on the electron. This use of indirect observation allows characterization of the dynamic nuclear polarization (DNP) process close to the electron. Frequency profiles of the electrondetected solid effect obtained using trityl radical showed intense saturation of the electron at the usual solid effect condition, which involves a single electron and nucleus. However, higher order solid effect transitions involving two, three, or four nuclei were also observed with surprising intensity, although these transitions did not lead to bulk nuclear polarization—suggesting that higher order transitions are important primarily in the transfer of polarization to nuclei nearby the electron. Similar results were obtained for the SABDPA radical where strong electronnuclear couplings produced splittings in the spectrum of the indirectly observed solid effect conditions. Observation of high order solid effect transitions supports recent studies of the solid effect, and suggests that a multispin solid effect mechanism may play a major role in polarization transfer via DNP.
 Atoms, Molecules, and Clusters

Shape resonances in lowenergyelectron collisions with halopyrimidines
View Description Hide DescriptionWe report calculated cross sections for elastic collisions of lowenergy electrons with halopyrimidines, namely, 2chloro, 2bromo, and 5bromopyrimidine. We employed the Schwinger multichannel method with pseudopotentials to compute the cross sections in the staticexchange and staticexchange plus polarization levels of approximation for energies up to 10 eV. We found four shape resonances for each molecule: three of π* nature localized on the ring and one of σ* nature localized along the carbon–halogen bond. We compared the calculated positions of the resonances with the electron transmission spectroscopy data measured by Modelli et al. [J. Phys. Chem. A115, 10775 (2011)]. In general the agreement between theory and experiment is good. In particular, our results show the existence of a π* temporary anion state of A 2 symmetry for all three halopyrimidines, in agreement with the dissociative electron attachment spectra also reported by Modelli et al. [J. Phys. Chem. A115, 10775 (2011)].

Simulation of femtosecond “doubleslit” experiments for a chromophore in a dissipative environment
View Description Hide DescriptionWe performed simulations of the prototypical femtosecond “doubleslit” experiment with strong pulsed laser fields for a chromophore in solution. The chromophore is modeled as a system with two electronic levels and a single FranckCondon active underdamped vibrational mode. All other (intra and intermolecular) vibrational modes are accounted for as a thermal bath. The systembath coupling is treated in a computationally accurate manner using the hierarchy equations of motion approach. The doubleslit signal is evaluated numerically exactly without invoking perturbation theory in the matterfield interaction. We show that the strongpulse doubleslit signal consists of a superposition of Nwavemixing (N = 2, 4, 6…) responses and can be split into population and coherence contributions. The former reveals the dynamics of vibrational wave packets in the ground state and the excited electronic state of the chromophore, while the latter contains information on the dephasing of electronic coherences of the chromophore density matrix. We studied the influence of heat baths with different coupling strengths and memories on the doubleslit signal. Our results show that the doubleslit experiment performed with strong (nonperturbative) pulses yields substantially more information on the photoinduced dynamics of the chromophore than the weakpulse experiment, in particular, if the bathinduced dephasings are fast.