Volume 145, Issue 8, 28 August 2016

We use molecular dynamics simulations to compute the spatially resolved static dielectric constant of water in cylindrical and spherical nanopores as occurring, e.g., in protein water pockets or carbon nanotubes. For this, we derive a linearresponse formalism which correctly takes into account the dielectric boundary conditions in the considered geometries. We find that in cylindrical confinement, the axial component behaves similar as the local density akin to what is known near planar interfaces. The radial dielectric constant shows some oscillatory features when approaching the surface if their radius is larger than about 2 nm. Most importantly, however, the radial component exhibits pronounced oscillations at the center of the cavity. These surprising features are traced back quantitatively to the nonlocal dielectric nature of bulk water.
 COMMUNICATIONS


Communication: Evaporation: Influence of heat transport in the liquid on the interface temperature and the particle flux
View Description Hide DescriptionMolecular dynamics simulations are reported for the evaporation of a liquid into vacuum, where a LennardJones type fluid with truncated and shifted potential at 2.5σ is considered. Vacuum is enforced locally by particle deletion and the liquid is thermostated in its bulk so that heat flows to the planar interface driving stationary evaporation. The length of the nonthermostated transition region between the bulk liquid and the interface Ln is under study. First, it is found for the reduced bulk liquid temperature Tl/Tc = 0.74 (Tc is the critical temperature) that by increasing Ln from 5.2σ to 208σ the interface temperature Ti drops by 17% and the evaporation flux decreases by a factor of 4.4. From a series of simulations for increasing values of Ln, an asymptotic value of the interface temperature for Ln → ∞ can be estimated which is 21% lower than the bulk liquid temperature Tl. Second, it is found that the evaporation flux is solely determined by the interface temperature Ti, independent on Tl or Ln. Combining these two findings, the evaporation coefficient α of a liquid thermostated on a macroscopic scale is estimated to be α ≈ 0.14 for Tl/Tc = 0.74.

Communication: Wigner functions in actionangle variables, BohrSommerfeld quantization, the Heisenberg correspondence principle, and a symmetrical quasiclassical approach to the full electronic density matrix
View Description Hide DescriptionIt is pointed out that the classical phase space distribution in actionangle (aa) variables obtained from a Wigner function depends on how the calculation is carried out: if one computes the standard Wigner function in Cartesian variables (p, x), and then replaces p and x by their expressions in terms of aa variables, one obtains a different result than if the Wigner function is computed directly in terms of the aa variables. Furthermore, the latter procedure gives a result more consistent with classical and semiclassical theory—e.g., by incorporating the BohrSommerfeld quantization condition (quantum states defined by integer values of the action variable) as well as the Heisenberg correspondence principle for matrix elements of an operator between such states—and has also been shown to be more accurate when applied to electronically nonadiabatic applications as implemented within the recently developed symmetrical quasiclassical (SQC) MeyerMiller (MM) approach. Moreover, use of the Wigner function (obtained directly) in aa variables shows how our standard SQC/MM approach can be used to obtain offdiagonal elements of the electronic density matrix by processing in a different way the same set of trajectories already used (in the SQC/MM methodology) to obtain the diagonal elements.

Communication: Variation after response in quantum Monte Carlo
View Description Hide DescriptionWe present a new method for modeling electronically excited states that overcomes a key failing of linear response theory by allowing the underlying ground state ansatz to relax in the presence of an excitation. The method is variational, has a cost similar to ground state variational Monte Carlo, and admits both open and periodic boundary conditions. We present preliminary numerical results showing that, when paired with the Jastrow antisymmetric geminal power ansatz, the variationafterresponse formalism delivers accuracies for valence and charge transfer single excitations on par with equation of motion coupled cluster, while surpassing coupled cluster’s accuracy for excitations with significant doubly excited character.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

Vibrational circular dichroism from ab initio molecular dynamics and nuclear velocity perturbation theory in the liquid phase
View Description Hide DescriptionWe report the first fully ab initio calculation of dynamical vibrational circular dichroism spectra in the liquid phase using nuclear velocity perturbation theory (NVPT) derived electronic currents. Our approach is rigorous and general and thus capable of treating weak interactions of chiral molecules as, e.g., chirality transfer from a chiral molecule to an achiral solvent. We use an implementation of the NVPT that is projected along the dynamics to obtain the current and magnetic dipole moments required for accurate intensities. The gauge problem in the liquid phase is resolved in a twofold approach. The electronic expectation values are evaluated in a distributed origin gauge, employing maximally localized Wannier orbitals. In a second step, the gauge invariant spectrum is obtained in terms of a scaled molecular moments, which allows to systematically include solvent effects while keeping a significant signaltonoise ratio. We give a thorough analysis and discussion of this choice of gauge for the liquid phase. At low temperatures, we recover the established double harmonic approximation. The methodology is applied to chiral molecules ((S)d2oxirane and (R)propyleneoxide) in the gas phase and in solution. We find an excellent agreement with the theoretical and experimental references, including the emergence of signals due to chirality transfer from the solute to the (achiral) solvent.

Accurate adiabatic singlettriplet gaps in atoms and molecules employing the thirdorder spinflip algebraic diagrammatic construction scheme for the polarization propagator
View Description Hide DescriptionFor the calculation of adiabatic singlettriplet gaps (STG) in diradicaloid systems the spinflip (SF) variant of the algebraic diagrammatic construction (ADC) scheme for the polarization propagator in third order perturbation theory (SFADC(3)) has been applied. Due to the methodology of the SF approach the singlet and triplet states are treated on an equal footing since they are part of the same determinant subspace. This leads to a systematically more accurate description of, e.g., diradicaloid systems than with the corresponding nonSF singlereference methods. Furthermore, using analytical excited state gradients at ADC(3) level, geometry optimizations of the singlet and triplet states were performed leading to a fully consistent description of the systems, leading to only small errors in the calculated STGs ranging between 0.6 and 2.4 kcal/mol with respect to experimental references.

Electrostatic interactions between charged dielectric particles in an electrolyte solution
View Description Hide DescriptionTheory is developed to address a significant problem of how two charged dielectric particles interact in the presence of a polarizable medium that is a dilute solution of a strong electrolyte. The electrostatic force is defined by characteristic parameters for the interacting particles (charge, radius, and dielectric constant) and for the medium (permittivity and Debye length), and is expressed in the form of a converging infinite series. The limiting case of weak screening and large interparticle separation is considered, which corresponds to small (macro)ions that carry constant charge. The theory yields a solution in the limit of monopole and dipole terms that agrees exactly with existing analytical expressions, which are generally used to describe ionion and ionmolecular interactions in a medium. Results from the theory are compared with DLVO theory and with experimental measurements for the electrostatic force between two PMMA particles contained in a nonpolar solvent (hexadecane) with an added charge control agent.

An adaptive interpolation scheme for molecular potential energy surfaces
View Description Hide DescriptionThe calculation of potential energy surfaces for quantum dynamics can be a time consuming task—especially when a high level of theory for the electronic structure calculation is required. We propose an adaptive interpolation algorithm based on polyharmonic splines combined with a partition of unity approach. The adaptive node refinement allows to greatly reduce the number of sample points by employing a local error estimate. The algorithm and its scaling behavior are evaluated for a model function in 2, 3, and 4 dimensions. The developed algorithm allows for a more rapid and reliable interpolation of a potential energy surface within a given accuracy compared to the nonadaptive version.

Linear and nonlinear dielectric theory for a slab: The connections between the phenomenological coefficients and the susceptibilities
View Description Hide DescriptionThe response of dielectric media to electromagnetic fields can be described by using either the response to a Maxwell field E or to an externally produced field . The former response is measured by phenomenological (dielectric) coefficients and the latter by susceptibilities. With the purpose of clarifying some recent proposals, the connections between the linear (twopoint) and first nonvanishing nonlinear (fourpoint) dielectric coefficients and the susceptibilities for media confined to a slab are examined using a general procedure developed sometime ago. Unlike the relations found for correlations between a local polarization density and the integrated polarization densities (total polarizations), the pointpoint connections give rise to nonvanishing cross correlations between polarization densities which are parallel and perpendicular to the slab surfaces. The cross correlations in the twopoint connections vanish when one polarization density is integrated to form the correlations between a local polarization density and the total polarization thereby losing angular information. The integrated parallel and perpendicular correlations remain different. When the fourpoint connections are similarly integrated most, but not all, cross correlations vanish. The angular correlations induced by the slab surfaces render the use of pointpoint correlations that are valid for isotropic media invalid for use in the integrated slab densities. In addition, the nonlinear fluctuations in the perpendicular components are drastically reduced relative to those in the parallel components or in isotropic media.

A cumulant functional for static and dynamic correlation
View Description Hide DescriptionA functional for the cumulant energy is introduced. The functional is composed of a paircorrection and static and dynamic correlation energy components. The paircorrection and static correlation energies are functionals of the natural orbitals and the occupancy transferred between neardegenerate orbital pairs, rather than the orbital occupancies themselves. The dynamic correlation energy is a functional of the statically correlated ontop twoelectron density. The ontop density functional used in this study is the wellknown ColleSalvetti functional. Using the ccpVTZ basis set, the functional effectively models the bond dissociation of H2, LiH, and N2 with equilibrium bond lengths and dissociation energies comparable to those provided by multireference secondorder perturbation theory. The performance of the cumulant functional is less impressive for HF and F2, mainly due to an underestimation of the dynamic correlation energy by the ColleSalvetti functional.

Condensed phase QM/MM simulations utilizing the exchange core functions to describe exchange repulsions at the QM boundary region
View Description Hide DescriptionIn a recent work, we developed a method [H. Takahashi et al., J. Chem. Phys. 143, 084104 (2015)] referred to as exchangecore function (ECF) approach, to compute exchange repulsion E ex between solute and solvent in the framework of the quantum mechanical (QM)/molecular mechanical (MM) method. The ECF, represented with a Slater function, plays an essential role in determining E ex on the basis of the overlap model. In the work of Takahashi et al. [J. Chem. Phys. 143, 084104 (2015)], it was demonstrated that our approach is successful in computing the hydrogen bond energies of minimal QM/MM systems including a cationic QM solute. We provide in this paper the extension of the ECF approach to the free energy calculation in condensed phase QM/MM systems by combining the ECF and the QM/MMER approach [H. Takahashi et al., J. Chem. Phys. 121, 3989 (2004)]. By virtue of the theory of solutions in energy representation, the free energy contribution δμ ex from the exchange repulsion was naturally formulated. We found that the ECF approach in combination with QM/MMER gives a substantial improvement on the calculation of the hydration free energy of a hydronium ion. This can be attributed to the fact that the ECF reasonably realizes the contraction of the electron density of the cation due to the deficit of an electron.

Locally adaptive method to define coordination shell
View Description Hide DescriptionAn algorithm is presented to define a particle’s coordination shell for any collection of particles. It requires only the particles’ positions and no preexisting knowledge or parameters beyond those already in the force field. A particle’s shell is taken to be all particles that are not blocked by any other particle and not further away than a blocked particle. Because blocking is based on two distances and an angle for triplets of particles, it is called the relative angular distance (RAD) algorithm. RAD is applied to LennardJones particles in molecular dynamics simulations of crystalline, liquid, and gaseous phases at various temperatures and densities. RAD coordination shells agree well with those from a cutoff in the radial distribution function for the crystals and liquids and are slightly higher for the gas.

An exact variational method to calculate rovibrational spectra of polyatomic molecules with large amplitude motion
View Description Hide DescriptionWe report a new fulldimensional variational algorithm to calculate rovibrational spectra of polyatomic molecules using an exact quantum mechanical Hamiltonian. The rovibrational Hamiltonian of system is derived in a set of orthogonal polyspherical coordinates in the bodyfixed frame. It is expressed in an explicitly Hermitian form. The Hamiltonian has a universal formulation regardless of the choice of orthogonal polyspherical coordinates and the number of atoms in molecule, which is suitable for developing a general program to study the spectra of many polyatomic systems. An efficient coupledstate approach is also proposed to solve the eigenvalue problem of the Hamiltonian using a multilayer Lanczos iterative diagonalization approach via a set of direct product basis set in three coordinate groups: radial coordinates, angular variables, and overall rotational angles. A simple set of symmetric top rotational functions is used for the overall rotation whereas a potentialoptimized discrete variable representation method is employed in radial coordinates. A set of contracted vibrationally diabatic basis functions is adopted in internal angular variables. Those diabatic functions are first computed using a neural network iterative diagonalization method based on a reduceddimension Hamiltonian but only once. The final rovibrational energies are computed using a modified Lanczos method for a given total angular momentum J, which is usually fast. Two numerical applications to CH4 and H2CO are given, together with a comparison with previous results.

Hartree potential dependent exchange functional
View Description Hide DescriptionWe introduce a novel nonlocal ingredient for the construction of exchange density functionals: the reduced Hartree parameter, which is invariant under the uniform scaling of the density and represents the exact exchange enhancement factor for one and twoelectron systems. The reduced Hartree parameter is used together with the conventional metageneralized gradient approximation (metaGGA) semilocal ingredients (i.e., the electron density, its gradient, and the kinetic energy density) to construct a new generation exchange functional, termed umetaGGA. This umetaGGA functional is exact for the exchange of any one and twoelectron systems, is sizeconsistent and nonempirical, satisfies the uniform density scaling relation, and recovers the modified gradient expansion derived from the semiclassical atom theory. For atoms, ions, jellium spheres, and molecules, it shows a good accuracy, being often better than metaGGA exchange functionals. Our construction validates the use of the reduced Hartree ingredient in exchangecorrelation functional development, opening the way to an additional rung in the Jacob’s ladder classification of nonempirical density functionals.

From plane waves to local Gaussians for the simulation of correlated periodic systems
View Description Hide DescriptionWe present a simple, robust, and blackbox approach to the implementation and use of local, periodic, atomcentered Gaussian basis functions within a plane wave code, in a computationally efficient manner. The procedure outlined is based on the representation of the Gaussians within a finite bandwidth by their underlying plane wave coefficients. The core region is handled within the projected augment wave framework, by pseudizing the Gaussian functions within a cutoff radius around each nucleus, smoothing the functions so that they are faithfully represented by a plane wave basis with only moderate kinetic energy cutoff. To mitigate the effects of the basis set superposition error and incompleteness at the meanfield level introduced by the Gaussian basis, we also propose a hybrid approach, whereby the complete occupied space is first converged within a large plane wave basis, and the Gaussian basis used to construct a complementary virtual space for the application of correlated methods. We demonstrate that these pseudized Gaussians yield compact and systematically improvable spaces with an accuracy comparable to their nonpseudized Gaussian counterparts. A key advantage of the described method is its ability to efficiently capture and describe electronic correlation effects of weakly bound and lowdimensional systems, where plane waves are not sufficiently compact or able to be truncated without unphysical artifacts. We investigate the accuracy of the pseudized Gaussians for the water dimer interaction, neon solid, and water adsorption on a LiH surface, at the level of secondorder Møller–Plesset perturbation theory.

Learning the mechanisms of chemical disequilibria
View Description Hide DescriptionWhen at equilibrium, largescale systems obey thermodynamics because they have microscopic configurations that are typical. “Typical” states are a fraction of those possible with the majority of the probability. A more precise definition of typical states underlies the transmission, coding, and compression of information. However, this definition does not apply to natural systems that are transiently away from equilibrium. Here, we introduce a variational measure of typicality and apply it to atomistic simulations of a model for hydrogen oxidation. While a gaseous mixture of hydrogen and oxygen combusts, reactant molecules transform through a variety of ephemeral species en route to the product, water. Out of the exponentially growing number of possible sequences of chemical species, we find that greater than 95% of the probability concentrates in less than 1% of the possible sequences. Overall, these results extend the notion of typicality across the nonequilibrium regime and suggest that typical sequences are a route to learning mechanisms from experimental measurements. They also open up the possibility of constructing ensembles for computing the macroscopic observables of systems out of equilibrium.

Geometric integrator for simulations in the canonical ensemble
View Description Hide DescriptionWe introduce a geometric integrator for molecular dynamics simulations of physical systems in the canonical ensemble that preserves the invariant distribution in equations arising from the density dynamics algorithm, with any possible type of thermostat. Our integrator thus constitutes a unified framework that allows the study and comparison of different thermostats and of their influence on the equilibrium and nonequilibrium (thermo)dynamic properties of a system. To show the validity and the generality of the integrator, we implement it with a secondorder, timereversible method and apply it to the simulation of a LennardJones system with three different thermostats, obtaining good conservation of the geometrical properties and recovering the expected thermodynamic results. Moreover, to show the advantage of our geometric integrator over a nongeometric one, we compare the results with those obtained by using the nongeometric Gear integrator, which is frequently used to perform simulations in the canonical ensemble. The nongeometric integrator induces a drift in the invariant quantity, while our integrator has no such drift, thus ensuring that the system is effectively sampling the correct ensemble.

The influence of the “cage” effect on the mechanism of reversible bimolecular multistage chemical reactions proceeding from different sites in solutions
View Description Hide DescriptionManifestations of the “cage” effect at the encounters of reactants have been theoretically treated on the example of multistage reactions (including bimolecular exchange reactions as elementary stages) proceeding from different active sites in liquid solutions. It is shown that for reactions occurring near the contact of reactants, consistent consideration of quasistationary kinetics of such multistage reactions (possible in the framework of the encounter theory only) can be made on the basis of chemical concepts of the “cage complex,” just as in the case of onesite model described in the literature. Exactly as in the onesite model, the presence of the “cage” effect gives rise to new channels of reactant transformation that cannot result from elementary event of chemical conversion for the given reaction mechanism. Besides, the multisite model demonstrates new (as compared to onesite model) features of multistage reaction course.

Eigenvector method for umbrella sampling enables error analysis
View Description Hide DescriptionUmbrella sampling efficiently yields equilibrium averages that depend on exploring rare states of a model by biasing simulations to windows of coordinate values and then combining the resulting data with physical weighting. Here, we introduce a mathematical framework that casts the step of combining the data as an eigenproblem. The advantage to this approach is that it facilitates error analysis. We discuss how the error scales with the number of windows. Then, we derive a central limit theorem for averages that are obtained from umbrella sampling. The central limit theorem suggests an estimator of the error contributions from individual windows, and we develop a simple and computationally inexpensive procedure for implementing it. We demonstrate this estimator for simulations of the alanine dipeptide and show that it emphasizes low free energy pathways between stable states in comparison to existing approaches for assessing error contributions. Our work suggests the possibility of using the estimator and, more generally, the eigenvector method for umbrella sampling to guide adaptation of the simulation parameters to accelerate convergence.

Reverse energy partitioning—An efficient algorithm for computing the density of states, partition functions, and free energy of solids
View Description Hide DescriptionA robust and model free Monte Carlo simulation method is proposed to address the challenge in computing the classical density of states and partition function of solids. Starting from the minimum configurational energy, the algorithm partitions the entire energy range in the increasing energy direction (“upward”) into subdivisions whose integrated density of states is known. When combined with the density of states computed from the “downward” energy partitioning approach [H. Do, J. D. Hirst, and R. J. Wheatley, J. Chem. Phys. 135, 174105 (2011)], the equilibrium thermodynamic properties can be evaluated at any temperature and in any phase. The method is illustrated in the context of the LennardJones system and can readily be extended to other molecular systems and clusters for which the structures are known.

Global structure search for molecules on surfaces: Efficient sampling with curvilinear coordinates
View Description Hide DescriptionEfficient structure search is a major challenge in computational materials science. We present a modification of the basin hopping global geometry optimization approach that uses a curvilinear coordinate system to describe global trial moves. This approach has recently been shown to be efficient in structure determination of clusters [C. Panosetti et al., Nano Lett. 15, 8044–8048 (2015)] and is here extended for its application to covalent, complex molecules and large adsorbates on surfaces. The employed automatically constructed delocalized internal coordinates are similar to molecular vibrations, which enhances the generation of chemically meaningful trial structures. By introducing flexible constraints and local translation and rotation of independent geometrical subunits, we enable the use of this method for molecules adsorbed on surfaces and interfaces. For two test systems, transβionylideneacetic acid adsorbed on a Au(111) surface and methane adsorbed on a Ag(111) surface, we obtain superior performance of the method compared to standard optimization moves based on Cartesian coordinates.