The sign problem and population dynamics in the full configuration interaction quantum Monte Carlo method

The recently proposed full configuration interaction quantum Monte Carlo method allows access to essentially exact ground-state energies of systems of interacting fermions substantially larger than previously tractable without knowledge of the nodal structure of the ground-state wave function. We investigate the nature of the sign problem in this method and how its severity depends on the system studied. We explain how cancellation of the positive and negative particles sampling the wave function ensures convergence to a stochastic representation of the many-fermion ground state and accounts for the characteristic population dynamics observed in simulations.     

Myoglobin-CO substate structures and dynamics: multidimensional vibrational echoes and molecular dynamics simulations

Spectrally resolved infrared stimulated vibrational echo data were obtained for sperm whale carbonmonoxymyoglobin (MbCO) at 300 K. The measured dephasing dynamics of the CO ligand are in agreement with dephasing dynamics calculated with molecular dynamics (MD) simulations for MbCO with the residue histidine-64 (His64) having its imidazole epsilon nitrogen protonated (N(epsilon)-H). The two conformational substate structures B(epsilon) and R(epsilon) observed in the MD simulations are assigned to the spectroscopic A(1) and A(3) conformational substates of MbCO, respectively, based on the agreement between the experimentally measured and calculated dephasing dynamics for these substates. In the A(1) substate, the N(epsilon)-H proton and N(delta) of His64 are approximately equidistant from the CO ligand, while in the A(3) substate, the N(epsilon)-H of His64 is oriented toward the CO, and the N(delta) is on the surface of the protein. The MD simulations show that dynamics of His64 represent the major source of vibrational dephasing of the CO ligand in the A(3) state on both femtosecond and picosecond time scales. Dephasing in the A(1) state is controlled by His64 on femtosecond time scales, and by the rest of the protein and the water solvent on longer time scales.     

Matching-pursuit/split-operator Fourier-transform simulations of nonadiabatic quantum dynamics

A rigorous and practical approach for simulations of nonadiabatic quantum dynamics is introduced. The algorithm involves a natural extension of the matching-pursuitsplit-operator Fourier-transform (MPSOFT) method [Y. Wu and V. S. Batista, J. Chem. Phys. 121, 1676 (2004)] recently developed for simulations of adiabatic quantum dynamics in multidimensional systems. The MPSOFT propagation scheme, extended to nonadiabatic dynamics, recursively applies the time-evolution operator as defined by the standard perturbation expansion to first-, or second-order, accuracy. The expansion is implemented in dynamically adaptive coherent-state representations, generated by an approach that combines the matching-pursuit algorithm with a gradient-based optimization method. The accuracy and efficiency of the resulting propagation method are demonstrated as applied to the canonical model systems introduced by Tully for testing simulations of dual curve-crossing nonadiabatic dynamics.     

Nonadiabatic molecular dynamics simulations of correlated electrons in solution. 1. Full configuration interaction (CI) excited-state relaxation dynamics of hydrated dielectrons

The hydrated dielectron is composed of two excess electrons dissolved in liquid water that occupy a single cavity; in both its singlet and triplet spin states there is a significant exchange interaction so the two electrons cannot be considered to be independent. In this paper and the following paper,we present the results of mixed quantum/classical molecular dynamics simulations of the nonadiabatic relaxation dynamics of photoexcited hydrated dielectrons, where we use full configuration interaction (CI) to solve for the two-electron wave function at every simulation time step. To the best of our knowledge, this represents the first systematic treatment of excited-state solvation dynamics where the multiple-electron problem is solved exactly. The simulations show that the effects of exchange and correlation contribute significantly to the relaxation dynamics. For example, spin-singlet dielectrons relax to the ground state on a time scale similar to that of single electrons excited at the same energy, but spin-triplet dielectrons relax much faster. The difference in relaxation dynamics is caused by exchange and correlation: The Pauli exclusion principle imposes very different electronic structure when the electrons' spins are singlet paired than when they are triplet paired, altering the available nonadiabatic relaxation pathways. In addition, we monitor how electronic correlation changes dynamically during nonadiabatic relaxation and show that solvent dynamics cause electron correlation to evolve quite differently for singlet and triplet dielectrons. Despite such differences, our calculations show that both spin states are stable to excited-state dissociation, but that the excited-state stability has different origins for the two spin states. For singlet dielectrons, the stability depends on whether the solvent structure can rearrange to create a second cavity before the ground state is reached. For triplet dielectrons, in contrast, electronic correlation ensures that the two electrons do not dissociate, even if the dielectron is artificially kept from reaching the ground state. In addition, both singlet and triplet dielectrons change shape dramatically during relaxation, so that linear response fails to describe the solvation dynamics for either spin state. In the following paper (Larsen, R. E.; Schwartz, B. J. J. Phys. Chem. B 2006, 110, 9692), we use these simulations to calculate the pump-probe spectroscopic signal expected for photoexcited hydrated dielectrons and to predict an experiment to observe hydrated dielectrons directly.     

Exciton migration dynamics in a dendritic molecule: quantum master equation approach using ab initio molecular orbital configuration interaction method

We investigate the exciton migration dynamics in a dendritic molecular model composed of pi-conjugation linear-leg units (acetylenes and diacetylene) and a benzene ring (branching point) using the quantum master equation approach with the ab initio molecular orbital (MO) configuration interaction (CI) method. The efficient migration of exciton from short-length linear legs (acetylenes) to long-length linear leg (diacetylene) via a benzene ring is observed. As predicted in previous studies, the exciton (electron and hole) distributions are relatively well localized in each generation segmented by the meta-branching point (meta-substituted benzene ring) though the electron and hole distributions are delocalized and are somewhat spatially different from each other within each generation. It is found that the excitons localized in the generation composed of short linear legs occupy in higher-lying exciton states, while those in the generation composed of long linear legs do in lower-lying ones. These features suggest the decoupling of pi-conjugation at the meta-branching point. On the other hand, the relaxation effect between exciton states is found to be caused by the exciton-phonon coupling, in which the existence of common configurations (electron-hole pairs) in CI wave functions between adjacent exciton states (having primary distributions on short and long linear-leg regions, respectively) is important for the relaxation between their exciton states. This feature indicates the importance of partial penetration of pi-conjugation through the meta-substituted benzene ring in excited states for such exciton migration.     

Hydrogen bond dynamics in heavy water studied with quantum dynamical simulations

The structure and dynamics of the hydrogen-bond network in heavy water (D(2)O) is studied as a function of the temperature using quantum dynamical simulations. Our approach combines an ab initio-based representation of the water interactions with an explicit quantum treatment of the molecular motion. A direct connection between the calculated linear and nonlinear vibrational spectra and the underlying molecular dynamics is made, which provides new insights into the rearrangement of the hydrogen-bond network in heavy water. A comparison with previous calculations on liquid H(2)O suggests that tunneling does not effectively contribute to the dynamics of the water hydrogen-bond network above the melting point. However, the effects of nuclear quantization are not negligible at all temperatures and become increasingly important near the melting point, which is in agreement with recent experimental analysis of the structural properties of liquid water as well as of the proton momentum distribution in supercooled water.     

Real space Hartree-Fock configuration interaction method for complex lateral quantum dot molecules

We present unrestricted Hartree-Fock method coupled with configuration interaction (CI) method (URHF-CI) suitable for the calculation of ground and excited states of large number of electrons localized by complex gate potentials in quasi-two-dimensional quantum dot molecules. The method employs real space finite difference method, incorporating strong magnetic field, for calculating single particle states. The Hartree-Fock method is employed for the calculation of direct and exchange interaction contributions to the ground state energy. The effects of correlations are included in energies and directly in the many-particle wave functions via CI method using a limited set of excitations above the Fermi level. The URHF-CI method and its performance are illustrated on the example of ten electrons confined in a two-dimensional quantum dot molecule.     

Thermal Gaussian molecular dynamics for quantum dynamics simulations of many-body systems: application to liquid para-hydrogen

A new method, here called thermal Gaussian molecular dynamics (TGMD), for simulating the dynamics of quantum many-body systems has recently been introduced [I. Georgescu and V. A. Mandelshtam, Phys. Rev. B 82, 094305 (2010)]. As in the centroid molecular dynamics (CMD), in TGMD the N-body quantum system is mapped to an N-body classical system. The associated both effective Hamiltonian and effective force are computed within the variational Gaussian wave-packet approximation. The TGMD is exact for the high-temperature limit, accurate for short times, and preserves the quantum canonical distribution. For a harmonic potential and any form of operator A?, it provides exact time correlation functions C(AB)(t) at least for the case of B, a linear combination of the position, x, and momentum, p, operators. While conceptually similar to CMD and other quantum molecular dynamics approaches, the great advantage of TGMD is its computational efficiency. We introduce the many-body implementation and demonstrate it on the benchmark problem of calculating the velocity time auto-correlation function for liquid para-hydrogen, using a system of up to N = 2592 particles.     

Trajectory-based solution of the nonadiabatic quantum dynamics equations: an on-the-fly approach for molecular dynamics simulations

The non-relativistic quantum dynamics of nuclei and electrons is solved within the framework of quantum hydrodynamics using the adiabatic representation of the electronic states. An on-the-fly trajectory-based nonadiabatic molecular dynamics algorithm is derived, which is also able to capture nuclear quantum effects that are missing in the traditional trajectory surface hopping approach based on the independent trajectory approximation. The use of correlated trajectories produces quantum dynamics, which is in principle exact and computationally very efficient. The method is first tested on a series of model potentials and then applied to study the molecular collision of H with H(2) using on-the-fly TDDFT potential energy surfaces and nonadiabatic coupling vectors.     

Nuclear quantum effects at aqueous metal interfaces captured by molecular dynamics simulations

In this review, we summarize the recent development in modeling nuclear quantum effects at aqueous metal interfaces. First, we review the nuclear quantum effects on the water-metal interface at ultrahigh vacuum. Then, we illustrate the nuclear quantum effects at the potential of zero charge conditions. At last, we give some outlook for the perspective work in modeling the nuclear quantum effects at electrochemical interfaces and some practical simulation strategies.     

Molecular dynamics simulations and quantum chemical calculations on β-cyclodextrin spironolactone complex

Molecular dynamics simulations on β-cyclodextrin in vacuo, with water and complexed with spironolactone (SP) were performed at a temperature of 300 K over a period of 1 ns. Two different orientations of SP in the cavity were considered. Along with conformational parameters, the formation of hydrogen bonds has been monitored during the whole simulation time. Cyclodextrins have the capability to form hydrogen bonds with the surrounding water molecules or intramolecular ones. The incorporation of ligands into the hydrophobic interior of β-cyclodextrin changes the preference of hydrogen bonds significantly and results in a contribution to the decrease of flexibility. Quantum chemical calculations on SP β-CD inclusion complex were performed to determine the interaction energy and to prove the applicability of various methods. Although all applied methods describe reasonable geometries for the association complex, higher level methods (e.g., B3LYP/6-31G(d,p)) seem to be necessary to determine reliable interaction energies.     

Nonspecific interaction forces at water-membrane interface by forced molecular dynamics simulations

Nonspecific interactions are the main driving forces for the behavior of molecules with great affinity for biologic membranes. To investigate not only the molecular details of these interactions but to estimate their magnitude as well, the theoretical method of Forced Molecular Dynamics Simulations, based on the Atomic Force Spectroscopy experimental technique, was applied. In this approach, an additional one-dimensional elastic force, representing the cantilever probe, was incorporated to the force field of a Molecular Dynamics computational program. This force represents a spring fixed on one end to a selected atom of the molecule; the other end of the spring is displaced at constant velocity to pull the molecule out of the membrane. The force experimented by the molecule due to the spring, is proportional to the spring elongation relative to its equilibrium position. This value is registered during the entire simulation, and its maximum value will determine the molecule-membrane interaction force. Nonexplicit medium simulations were carried out. Polar and apolar media were considered according to their polarizability degree and a specific dielectric constant value was assigned. In this approach, the membrane was considered as the apolar region limited by two flat surfaces with a polar aqueous medium. The potential energy discontinuity at the interfaces was smoothed by considering the polarization-induced effects using the image method. The results of this methodology are presented using a small system, a single Alanine amino acid model, which enables extended simulations in a microsecond time scale. The confinement of this amino acid at the interface reduces its degrees of freedom and forces it to adopt one of the six defined conformations. A correlation between these stable structures at the water-membrane interface and the interaction force value was determined.     

Finite element interpolation for combined classical/quantum mechanical molecular dynamics simulations

A method is presented to interpolate the potential energy function for a part of a system consisting of a few degrees of freedom, such as a molecule in solution. The method is based on a modified finite element (FE) interpolation scheme. The aim is to save computer time when expensive methods such as quantum-chemical calculations are used to determine the potential energy function. The expensive calculations are only carried out if the molecule explores new unknown regions of the conformation space. If the molecule resides in regions previously explored, a cheap interpolation is performed instead of an expensive calculation, using known neighboring points. We report the interpolation techniques for the energies and the forces of the molecule, the handling of the FE mesh, and an application to a simple test example in molecular dynamics (MD) simulations. Good performance of the method was obtained (especially for MD simulations with a preceding Monte Carlo mesh generation) without losing accuracy. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1484–1495, 1997     

Brownian dynamics simulations of polyelectrolyte adsorption in shear flow with hydrodynamic interaction

The adsorption of single polyelectrolyte molecules in shear flow is studied using Brownian dynamics simulations with hydrodynamic interaction (HI). Simulations are performed with bead-rod and bead-spring chains, and electrostatic interactions are incorporated through a screened Coulombic potential with excluded volume accounted for by the repulsive part of a Lennard-Jones potential. A correction to the Rotne-Prager-Yamakawa tensor is derived that accounts for the presence of a planar wall. The simulations show that migration away from an uncharged wall, which is due to bead-wall HI, is enhanced by increases in the strength of flow and intrachain electrostatic repulsion, consistent with kinetic theory predictions. When the wall and polyelectrolyte are oppositely charged, chain behavior depends on the strength of electrostatic screening. For strong screening, chains get depleted from a region close to the wall and the thickness of this depletion layer scales as N(1/3)Wi(2/3) at high Wi, where N is the chain length and Wi is the Weissenberg number. At intermediate screening, bead-wall electrostatic attraction competes with bead-wall HI, and it is found that there is a critical Weissenberg number for desorption which scales as N(-1/2)kappa(-3)(l(B)|sigmaq|)(3/2), where kappa is the inverse screening length, l(B) is the Bjerrum length, sigma is the surface charge density, and q is the bead charge. When the screening is weak, adsorbed chains are observed to align in the vorticity direction at low shear rates due to the effects of repulsive intramolecular interactions. At higher shear rates, the chains align in the flow direction. The simulation method and results of this work are expected to be useful for a number of applications in biophysics and materials science in which polyelectrolyte adsorption plays a key role.     

Conformational properties of penicillins: quantum chemical calculations and molecular dynamics simulations of benzylpenicillin

Herein, we present theoretical results on the conformational properties of benzylpenicillin, which are characterized by means of quantum chemical calculations (MP2/6-31G* and B3LYP/6-31G*) and classical molecular dynamics simulations (5 ns) both in the gas phase and in aqueous solution. In the gas phase, the benzylpenicillin conformer in which the thiazolidine ring has the carboxylate group oriented axially is the most favored one. Both intramolecular CH. O and dispersion interactions contribute to stabilize the axial conformer with respect to the equatorial one. In aqueous solution, a molecular dynamics simulation predicts a relative population of the axial:equatorial conformers of 0.70:0.30 in consonance with NMR experimental data. Overall, the quantum chemical calculations as well as the simulations give insight into substituent effects, the conformational dynamics of benzylpenicillin, the frequency of ring-puckering motions, and the correlation of side chain and ring-puckering motions.     

Phenol-benzene complexation dynamics: quantum chemistry calculation, molecular dynamics simulations, and two dimensional IR spectroscopy

Molecular dynamics (MD) simulations and quantum mechanical electronic structure calculations are used to investigate the nature and dynamics of the phenol-benzene complex in the mixed solvent, benzene/CCl4. Under thermal equilibrium conditions, the complexes are continuously dissociating and forming. The MD simulations are used to calculate the experimental observables related to the phenol hydroxyl stretching mode, i.e., the two dimensional infrared vibrational echo spectrum as a function of time, which directly displays the formation and dissociation of the complex through the growth of off-diagonal peaks, and the linear absorption spectrum, which displays two hydroxyl stretch peaks, one for the complex and one for the free phenol. The results of the simulations are compared to previously reported experimental data and are found to be in quite reasonable agreement. The electronic structure calculations show that the complex is T shaped. The classical potential used for the phenol-benzene interaction in the MD simulations is in good accord with the highest level of the electronic structure calculations. A variety of other features is extracted from the simulations including the relationship between the structure and the projection of the electric field on the hydroxyl group. The fluctuating electric field is used to determine the hydroxyl stretch frequency-frequency correlation function (FFCF). The simulations are also used to examine the number distribution of benzene and CCl4 molecules in the first solvent shell around the phenol. It is found that the distribution is not that of the solvent mole fraction of benzene. There are substantial probabilities of finding a phenol in either a pure benzene environment or a pure CCl4 environment. A conjecture is made that relates the FFCF to the local number of benzene molecules in phenol's first solvent shell.     

One- and multidimensional conformational free energy simulations

A new thermodynamic integration approach to conformational free energy simulations is presented. The method is applicable both to one-dimensional cases (reaction coordinates) and multidimensional situations (free energy surfaces). Analysis of the properties of the thermodynamic integration algorithm is used to formulate methods of calculating multidimensional free energy gradients. The method is applied to calculate the free energy profile for rotation around the central C—C bond of n-butane in the gas and liquid phase and to generate maps of the 18-dimensional free energy gradient with respect to all nine ϕ and nine ψ dihedrals of the decaalanine and deca-α-methylalanine peptides in vacuum. For n-butane essentially no change in the gauche–trans equilibrium between the gas and liquid is predicted within the CHARMM explicit hydrogen model, with the thermodynamic integration, thermodynamic perturbation, and direct simulation methods yielding free energy profiles that are identical within errors. For the decapeptides the right-handed helical region of conformational space is investigated. For decaalanine a minimum on the free energy surface is found in the vicinity of (ϕ, ψ) = (-64.5°, -42.5°) in the α-helix region; no minimum exists for 3₁₀-helix-type conformers. For deca-α-methylalanine free energy minima corresponding to both the α-helix at ( - 55.5°, - 51.5°) and the 3₁₀-helix at ( - 54°, - 29°) are found; the α-helix state is favored by about 4 kcal/mol and the barrier for the concerted 3₁₀-helix → α-helix transition is about 3 kcal/mol. The α-methylation also considerably increases the rigidity of the α-helix with respect to deformations. The computational efficiency, ease of generalization to calculations of multidimensional gradients, and analytical capability due to component analysis of free energy differences make the method a novel, powerful tool to improve the basic understanding of conformational equilibria of flexible molecules in condensed phases. A related scheme for energy minimization in the presence of holonomic constraints is also presented, allowing generation of adiabatic energy surfaces in constrained systems. © 1996 by John Wiley & Sons, Inc.