Optimal control of classical molecular dynamics: A perturbation formulation and the existence of multiple solutions

This paper considers the prospect for there being multiple solutions to the control of classically modelled molecular dynamical systems. The research presented here follows up on a parallel study based on quantum mechanics. For polyatomic molecules it is generally expected that a classical mechanical model will be adequate and necessary as a means for designing optical fields for molecular control. The prospect for multiple control field solutions existing in this domain is important to establish in terms of ultimate laboratory realization of molecular control. A general formulation of the multiplicity problem is considered and the existence of a denumerably infinite number of solutions for the control field amplitude is shown to be the case under certain mild limitations on the physical variables.     

Energy conservation in molecular dynamics simulations of classical systems

Classical Newtonian dynamics is analytic and the energy of an isolated system is conserved. The energy of such a system, obtained by the discrete "Verlet" algorithm commonly used in molecular dynamics simulations, fluctuates but is conserved in the mean. This is explained by the existence of a "shadow Hamiltonian" H [S. Toxvaerd, Phys. Rev. E 50, 2271 (1994)], i.e., a Hamiltonian close to the original H with the property that the discrete positions of the Verlet algorithm for H lie on the analytic trajectories of H. The shadow Hamiltonian can be obtained from H by an asymptotic expansion in the time step length. Here we use the first non-trivial term in this expansion to obtain an improved estimate of the discrete values of the energy. The investigation is performed for a representative system with Lennard-Jones pair interactions. The simulations show that inclusion of this term reduces the standard deviation of the energy fluctuations by a factor of 100 for typical values of the time step length. Simulations further show that the energy is conserved for at least one hundred million time steps provided the potential and its first four derivatives are continuous at the cutoff. Finally, we show analytically as well as numerically that energy conservation is not sensitive to round-off errors.     

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     

First principles and classical molecular dynamics simulations of solvated benzene

We have performed extensive ab initio and classical molecular dynamics (MD) simulations of benzene in water in order to examine the unique solvation structures that are formed. Qualitative differences between classical and ab initio MD simulations are found and the importance of various technical simulation parameters is examined. Our comparison indicates that nonpolarizable classical models are not capable of describing the solute-water interface correctly if local interactions become energetically comparable to water hydrogen bonds. In addition, a comparison is made between a rigid water model and fully flexible water within ab initio MD simulations which shows that both models agree qualitatively for this challenging system.     

Transport properties of room-temperature ionic liquids from classical molecular dynamics

Room-temperature ionic liquids (RTILs) have attracted much attention in the scientific community in the past decade due their novel and highly customizable properties. Nonetheless, their high viscosities pose serious limitations to the use of RTILs in practical applications. To elucidate some of the physical aspects behind transport properties of RTILs, extensive classical molecular dynamics calculations are reported. Here, in particular, bulk viscosities and ionic conductivities of butyl-methyl-imidazole based RTILs are presented over a wide range of temperatures. The dependence of the properties of the liquids on simulation parameters, e.g., system-size effects or the choice of the interaction potential, is analyzed in detail.     

Optimal laser control of ultrafast photodissociation of I2- in water: mixed quantum/classical molecular dynamics simulation

A linearized optimal control method in combination with mixed quantum/classical molecular dynamics simulation is used for numerically investigating the possibility of controlling photodissociation wave packets of I(2)(-) in water. Optimal pulses are designed using an ensemble of photodissociation samples, aiming at the creation of localized dissociation wave packets. Numerical results clearly show the effectiveness of the control although the control achievement is reduced with an increase in the internuclear distance associated with a target region. We introduce effective optimal pulses that are designed using a statistically averaged effective dissociation potential, and show that they semiquantitatively reproduce the control achievements calculated by using optimal pulses. The control mechanisms are interpreted from the time- and frequency-resolved spectra of the effective optimal pulses.     

Calculating chemically accurate redox potentials for engineered flavoproteins from classical molecular dynamics free energy simulations

The tricyclic isoalloxazine nucleus of the redox cofactors flavin mononucleotide (FMN) and flavin adenine dinucleotide (FAD) acts as an electron sink in life-sustaining biological electron transfer (eT). The functional diversity of flavin-containing proteins (flavoproteins) transcends that of free flavins. A large body of experimental evidence attributes natural control of flavoprotein-mediated eT to tuning of the thermodynamic driving force by the protein environment. Understanding and engineering such modulation by the protein environment of the flavin redox potential (DeltaE(o)) is valuable in biotechnology and device design. In this study we employed classical molecular dynamics free energy simulations (MDFES), within a thermodynamic integration (TI) formalism, to calculate the change in FMN first reduction potential (DeltaDeltaE(o)(ox/sq)) imparted by 6 flavoprotein active site mutations. The combined performance of the AMBER ff03 (protein) and GAFF (cofactor) force fields was benchmarked against experimental data for mutations close to the isoalloxazine re- and si-faces that perturb the wild-type DeltaE(o)(ox/sq) value in Anabaena flavodoxin. The classical alchemical approach used in this study overestimates the magnitude of DeltaE(o) values, in common with other studies. Nevertheless, chemically accurate DeltaDeltaE(o) values--calculated to within 1 kcal mol(-1) of the experimental value--were obtained for five of the six mutations studied. We have shown that this approach is practical for quantitative in silico screening of the effect of mutations on the first reduction potential where experimental values and structural data are available for the wild-type flavoprotein. This approach promises to be useful as an integral part of future interdisciplinary strategies to engineer desired thermodynamic properties in flavoproteins of biotechnological interest.     

Asymptotic extrapolation scheme for large-scale calculations with hybrid coupled cluster and molecular dynamics simulations

In this paper we discuss a simple extrapolation scheme based on the asymptotic behavior of the electronic energies considered as functions of cutoff factor for orbital energies corresponding to virtual orbitals. The performance of this approach is illustrated in the context of large-scale dynamic simulations for excitation energies of the cytosine molecule in its native DNA environment. We demonstrate that the extrapolation errors are significantly smaller than the excitation-energy fluctuations, due to the fluctuating environment.     

Hybrid quantum/classical molecular dynamics for a proton transfer reaction coupled to a dissipative bath

A hybrid quantum/classical molecular dynamics approach is applied to a proton transfer reaction represented by a symmetric double well system coupled to a dissipative bath. In this approach, the proton is treated quantum mechanically and all bath modes are treated classically. The transition state theory rate constant is obtained from the potential of mean force, which is generated along a collective reaction coordinate with umbrella sampling techniques. The transmission coefficient, which accounts for dynamical recrossings of the dividing surface, is calculated with a reactive flux approach combined with the molecular dynamics with quantum transitions surface hopping method. The hybrid quantum/classical results agree well with numerically exact results in the spatial-diffusion-controlled regime, which is most relevant for proton transfer in proteins. This hybrid quantum/classical approach has already been shown to be computationally practical for studying proton transfer in large biological systems. These results have important implications for future applications to hydrogen transfer reactions in solution and proteins.     

Controlling the shape and flexibility of arylamides: a combined ab initio, ab initio molecular dynamics, and classical molecular dynamics study

Using quantum chemistry plus ab initio molecular dynamics and classical molecular dynamics methods, we address the relationship between molecular conformation and the biomedical function of arylamide polymers. Specifically, we have developed new torsional parameters for a class of these polymers and applied them in a study of the interaction between a representative arylamide and one of its biomedical targets, the anticoagulant drug heparin. Our main finding is that the torsional barrier of a C(aromatic)-C(carbonyl) bond increases significantly upon addition of an o-OCH2CH2NH3+ substituent on the benzene ring. Our molecular dynamics studies that are based on the original general AMBER force field (GAFF) and GAFF modified to include our newly developed torsional parameters show that the binding mechanism between the arylamide and heparin is very sensitive to the choice of torsional potentials. Ab initio molecular dynamics simulation of the arylamide independently confirms the degree of flexibility we obtain by classical molecular dynamics when newly developed torsional potentials are used.     

Bidirectional molecular dynamics: Interpretation in terms of a modern formulation of classical mechanics

The validity and applicability of bidirectional molecular dynamics is shown in terms of modern classical mechanics. A simple interpretation of bidirectional molecular dynamics is given. This interpretation justifies an easy approach to start from an equilibrium configuration, to perform simulations in parallel in “forward” and “backward” time direction, and to combine the two obtained trajectories into a single one for proper evaluation of statistical quantities. Practical results, obtained for liquid ammonia, are presented. © 1996 by John Wiley & Sons, Inc.     

Combined ab initio and classical molecular dynamics simulations of alkyl-lithium aggregates in ethereal solutions

Molecular dynamics simulations of organolithium aggregates in solution are reported for the first time. We use a combined quantum/classical force field (the so-called QM/MM approach) and study ethyl-lithium aggregates in dimethyl ether (DME) solvent. The solutes are described at the Density Functional Theory level while solvent molecules are described using molecular mechanics. NVT Molecular Dynamics simulations at 200 K are carried out in the Born–Oppenheimer approximation. After equilibration, the production phase was run for 80 ps (monomer), 40 ps (dimer) and 26 ps (tetramer). The analysis of the results focuses on Li coordination as a function of aggregate size and we show that the total Li coordination number is always 4. No decoordination has been observed along the simulations. Fluctuations of the structures are predicted to be large in some cases and possible implications on reactivity are discussed.     

Structure of liquid formic acid investigated by first principle and classical molecular dynamics simulations

The structure of liquid formic acid has been investigated by Car-Parrinello and classical molecular dynamics simulations, focusing on the characterization of the H-bond network and on the mutual arrangement of pairs of bonded molecules. In agreement with previous computational studies, two levels of H-bonded structures have been found. Small clusters, characterized by O-H...O bonds, are held together by weak C-H...O bonds to form large branched structures. From the ab initio simulation we infer the importance of cyclic H-bond dimer configurations, typical of the gas phase. Most of these dimer structures are however found to be embedded into H-bonded chains. When only O-H...O bonds are taken into account, linear H-bond chains are detected as basic structures of the liquid. More branched structures occur when C-H...O bonds are also considered. Regarding the arrangement of molecular pairs, we observed that O-H...O bonds favor the occurrence of configurations with parallel molecular planes, whereas no preferential orientation is observed for molecules forming C-H...O bonds.     

Recent applications and developments of charge equilibration force fields for modeling dynamical charges in classical molecular dynamics simulations

With the continuing advances in computational hardware and novel force fields constructed using quantum mechanics, the outlook for non-additive force fields is promising. Our work in the past several years has demonstrated the utility of polarizable force fields, in our hands those based on the charge equilibration formalism, for a broad range of physical and biophysical systems. We have constructed and applied polarizable force fields for small molecules, proteins, lipids, and lipid bilayers and recently have begun work on carbohydrate force fields. The latter area has been relatively untouched by force field developers with particular focus on polarizable, non-additive interaction potential models. In this review of our recent work, we discuss the formalism we have adopted for implementing the charge equilibration method for phase-dependent polarizable force fields, lipid molecules, and small-molecule carbohydrates. We discuss the methodology, related issues, and briefly discuss results from recent applications of such force fields.     

Algorithms for constrained molecular dynamics

In molecular dynamics simulations, the fastest components of the potential field impose severe restrictions on the stability and hence the speed of computational methods. One possibility for treating this problem is to replace the fastest components with algebraic length constraints. In this article the resulting systems of mixed differential and algebraic equations are studied. Commonly used discretization schemes for constrained Hamiltonian systems are discussed. The form of the nonlinear equations is examined in detail and used to give convergence results for the traditional nonlinear solution technique SHAKE iteration and for a modification based on successive overrelaxation (SOR). A simple adaptive algorithm for finding the optimal relaxation parameter is presented. Alternative direct methods using sparse matrix techniques are discussed. Numerical results are given for the new techniques, which have been implemented in the molecular modeling software package CHARMM and show as much as twofold improvement over SHAKE iteration. © 1995 John Wiley & Sons, Inc.     

Quantum statistics and classical mechanics: real time correlation functions from ring polymer molecular dynamics

We propose an approximate method for calculating Kubo-transformed real-time correlation functions involving position-dependent operators, based on path integral (Parrinello-Rahman) molecular dynamics. The method gives the exact quantum mechanical correlation function at time zero, exactly satisfies the quantum mechanical detailed balance condition, and for correlation functions of the form C(Ax)(t) and C(xB)(t) it gives the exact result for a harmonic potential. It also works reasonably well at short times for more general potentials and correlation functions, as we illustrate with some example calculations. The method provides a consistent improvement over purely classical molecular dynamics that is most apparent in the low-temperature regime.     

Investigating the quartz (1010)/water interface using classical and ab initio molecular dynamics

Two different terminations of the (1010) surface of quartz (α and β) interacting with water are simulated by classical (CMD) (using two different force fields) and ab initio molecular dynamics (AIMD) and compared with previously published X-ray reflectivity (XR) experiments. Radial distribution functions between hydroxyl and water show good agreement between AIMD and CMD using the ClayFF force field for both terminations. The Lopes et al. (Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D. J. Phys. Chem. B2006, 110, 2782-2792) force field (LFF), however, underestimates the extent of hydroxyl-water hydrogen bonding. The β termination is found to contain hydroxyl-hydroxyl hydrogen bonds; the quartz surface hydroxyl hydrogens and oxygens that hydrogen bond with each other exhibit greatly reduced hydrogen bonding to water. Conversely, the hydroxyl hydrogen and oxygens that are not hydrogen bonded to other surface hydroxyls but are connected to those that are show a considerable amount of hydrogen bonding to water. The electron density distribution of an annealed surface of quartz (1010) obtained by XR is in qualitative agreement with electron densities calculated by CMD and AIMD. In all simulation methods, the interfacial water peak appears farther from the surface than observed by XR. Agreement among AIMD, LFF, and XR is observed for the relaxation of the near-surface atoms; however, ClayFF shows a larger discrepancy. Overall, results show that for both terminations of (1010), LFF treats the near-surface structure more accurately whereas ClayFF treats the interfacial water structure more accurately. It is shown that the number of hydroxyl and water hydrogen bonds to the bridging Si-O-Si oxygens connecting the surface silica groups to the rest of the crystal is much greater for the α than the β termination. It is suggested that this may play a role in the greater resistance to dissolution of the β termination than that of the α termination.