Volume 137, Issue 4, 28 July 2012

The displacement of perturbed water upon binding is believed to play a critical role in the thermodynamics of biomolecular recognition, but it is nontrivial to unambiguously define and answer questions about this process. We address this issue by introducing grid inhomogeneous solvation theory (GIST), which discretizes the equations of inhomogeneous solvation theory (IST) onto a threedimensional grid situated in the region of interest around a solute molecule or complex. Snapshots from explicit solvent simulations are used to estimate localized solvation entropies, energies, and free energies associated with the grid boxes, or voxels, and properly summing these thermodynamic quantities over voxels yields information about hydration thermodynamics. GIST thus provides a smoothly varying representation of water properties as a function of position, rather than focusing on hydration sites where solvent is present at high density. It therefore accounts for full or partial displacement of water from sites that are highly occupied by water, as well as for partly occupied and waterdepleted regions around the solute. GIST can also provide a welldefined estimate of the solvation free energy and therefore enables a rigorous endstates analysis of binding. For example, one may not only use a first GIST calculation to project the thermodynamic consequences of displacing water from the surface of a receptor by a ligand, but also account, in a second GIST calculation, for the thermodynamics of subsequent solvent reorganization around the bound complex. In the present study, a first GIST analysis of the molecular host cucurbit[7]uril is found to yield a rich picture of hydration structure and thermodynamics in and around this miniature receptor. One of the most striking results is the observation of a toroidal region of high water density at the center of the host's nonpolar cavity. Despite its high density, the water in this toroidal region is disfavored energetically and entropically, and hence may contribute to the known ability of this small receptor to bind guest molecules with unusually high affinities. Interestingly, the toroidal region of high water density persists even when all partial charges of the receptor are set to zero. Thus, localized regions of high solvent density can be generated in a binding site without strong, attractive solutesolvent interactions.
 COMMUNICATIONS


Communication: Hydration structure and polarization of heavy alkali ions: A first principles molecular dynamics study of Rb^{+} and Cs^{+}
View Description Hide DescriptionHydration structure and polarization of Rb^{+} and Cs^{+} in liquid water at ambient conditions were studied by first principles molecular dynamics. Our systematic analysis of the relevant electronic structures, based on maximally localized Wannier functions, revealed that the dipole moment of H_{2}O molecules in the first solvation shell of the ions slightly increases with increasing the atomic number. We also found that the polarization of heavy alkali ions, particularly Cs^{+}, tends to stabilize a peculiar asymmetric hydration structure with relevant consequences in the extraction of the harmful ^{137}Cs resulting from nuclear wastes.
 Top

 ARTICLES

 Theoretical Methods and Algorithms

Grid inhomogeneous solvation theory: Hydration structure and thermodynamics of the miniature receptor cucurbit[7]uril
View Description Hide DescriptionThe displacement of perturbed water upon binding is believed to play a critical role in the thermodynamics of biomolecular recognition, but it is nontrivial to unambiguously define and answer questions about this process. We address this issue by introducing grid inhomogeneous solvation theory (GIST), which discretizes the equations of inhomogeneous solvation theory (IST) onto a threedimensional grid situated in the region of interest around a solute molecule or complex. Snapshots from explicit solvent simulations are used to estimate localized solvation entropies, energies, and free energies associated with the grid boxes, or voxels, and properly summing these thermodynamic quantities over voxels yields information about hydration thermodynamics. GIST thus provides a smoothly varying representation of water properties as a function of position, rather than focusing on hydration sites where solvent is present at high density. It therefore accounts for full or partial displacement of water from sites that are highly occupied by water, as well as for partly occupied and waterdepleted regions around the solute. GIST can also provide a welldefined estimate of the solvation free energy and therefore enables a rigorous endstates analysis of binding. For example, one may not only use a first GIST calculation to project the thermodynamic consequences of displacing water from the surface of a receptor by a ligand, but also account, in a second GIST calculation, for the thermodynamics of subsequent solvent reorganization around the bound complex. In the present study, a first GIST analysis of the molecular host cucurbit[7]uril is found to yield a rich picture of hydration structure and thermodynamics in and around this miniature receptor. One of the most striking results is the observation of a toroidal region of high water density at the center of the host's nonpolar cavity. Despite its high density, the water in this toroidal region is disfavored energetically and entropically, and hence may contribute to the known ability of this small receptor to bind guest molecules with unusually high affinities. Interestingly, the toroidal region of high water density persists even when all partial charges of the receptor are set to zero. Thus, localized regions of high solvent density can be generated in a binding site without strong, attractive solutesolvent interactions.

Nonequilibrium molecular dynamics simulation of water transport through carbon nanotube membranes at low pressure^{a)}
View Description Hide DescriptionNonequilibrium molecular dynamics (NEMD) simulations are used to investigate pressuredriven water flow passing through carbon nanotube(CNT) membranes at low pressures (5.0 MPa) typical of real nanofiltration (NF) systems. The CNT membrane is modeled as a simplified NF membrane with smooth surfaces, and uniform straight pores of typical NF pore sizes. A NEMD simulation system is constructed to study the effects of the membrane structure (pores size and membrane thickness) on the pure watertransport properties. All simulations are run under operating conditions (temperature and pressure difference) similar to a real NF processes. Simulation results are analyzed to obtain water flux, density, and velocity distributions along both the flow and radial directions. Results show that water flow through a CNT membrane under a pressure difference has the unique transport properties of very fast flow and a nonparabolic radial distribution of velocities which cannot be represented by the HagenPoiseuille or NavierStokes equations. Density distributions along radial and flow directions show that water molecules in the CNT form layers with an oscillatory density profile, and have a lower average density than in the bulk flow. The NEMD simulations provide direct access to dynamic aspects of water flow through a CNT membrane and give a view of the pressuredriven transport phenomena on a molecular scale.

Tensor hypercontraction density fitting. I. Quartic scaling second and thirdorder MøllerPlesset perturbation theory
View Description Hide DescriptionMany approximations have been developed to help deal with the O(N ^{ 4 } ) growth of the electron repulsion integral (ERI) tensor, where N is the number of oneelectron basis functions used to represent the electronic wavefunction. Of these, the density fitting (DF) approximation is currently the most widely used despite the fact that it is often incapable of altering the underlying scaling of computational effort with respect to molecular size. We present a method for exploiting sparsity in threecenter overlap integrals through tensor decomposition to obtain a lowrank approximation to density fitting (tensor hypercontraction density fitting or THCDF). This new approximation reduces the 4thorder ERI tensor to a product of five matrices, simultaneously reducing the storage requirement as well as increasing the flexibility to regroup terms and reduce scaling behavior. As an example, we demonstrate such a scaling reduction for second and thirdorder perturbation theory (MP2 and MP3), showing that both can be carried out in O(N ^{ 4 } ) operations. This should be compared to the usual scaling behavior of O(N ^{ 5 } ) and O(N ^{ 6 } ) for MP2 and MP3, respectively. The THCDF technique can also be applied to other methods in electronic structuretheory, such as coupledcluster and configuration interaction, promising significant gains in computational efficiency and storage reduction.

A multiconfigurational hybrid densityfunctional theory
View Description Hide DescriptionWe propose a multiconfigurational hybrid densityfunctional theory which rigorously combines a multiconfiguration selfconsistentfield calculation with a densityfunctional approximation based on a linear decomposition of the electronelectron interaction. This gives a straightforward extension of the usual hybrid approximations by essentially adding a fraction λ of exact static correlation in addition to the fraction λ of exact exchange. Test calculations on the cycloaddition reactions of ozone with ethylene or acetylene and the dissociation of diatomic molecules with the PerdewBurkeErnzerhof and BeckeLeeYangParr density functionals show that a good value of λ is 0.25, as in the usual hybrid approximations. The results suggest that the proposed multiconfigurational hybrid approximations can improve over usual densityfunctional calculations for situations with strong static correlation effects.

On the precision of quasi steady state assumptions in stochastic dynamics
View Description Hide DescriptionMany biochemical networks have complex multidimensional dynamics and there is a long history of methods that have been used for dimensionality reduction for such reaction networks. Usually a deterministic mass action approach is used; however, in small volumes, there are significant fluctuations from the mean which the mass action approach cannot capture. In such cases stochastic simulation methods should be used. In this paper, we evaluate the applicability of one such dimensionality reduction method, the quasisteady state approximation (QSSA) [L. Menten and M. Michaelis, “Die kinetik der invertinwirkung,” Biochem. Z49, 333369 (1913)] for dimensionality reduction in case of stochastic dynamics. First, the applicability of QSSA approach is evaluated for a canonical system of enzyme reactions. Application of QSSA to such a reaction system in a deterministic setting leads to MichaelisMenten reduced kinetics which can be used to derive the equilibrium concentrations of the reaction species. In the case of stochastic simulations, however, the steady state is characterized by fluctuations around the mean equilibrium concentration. Our analysis shows that a QSSA based approach for dimensionality reduction captures well the mean of the distribution as obtained from a full dimensional simulation but fails to accurately capture the distribution around that mean. Moreover, the QSSA approximation is not unique. We have then extended the analysis to a simple bistable biochemical network model proposed to account for the stability of synaptic efficacies; the substrate of learning and memory [J. E. Lisman, “A mechanism of memory storage insensitive to molecular turnover: A bistable autophosphorylating kinase,” Proc. Natl. Acad. Sci. U.S.A.82, 3055–3057 (1985)]. Our analysis shows that a QSSA based dimensionality reduction method results in errors as big as two orders of magnitude in predicting the residence times in the two stable states.

Replica exchanging selfguided Langevin dynamics for efficient and accurate conformational sampling
View Description Hide DescriptionThis work presents a replica exchanging selfguided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperaturebased replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses selfguided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and/or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.

A computationally efficacious freeenergy functional for studies of inhomogeneous liquid water
View Description Hide DescriptionWe present an accurate equation of state for water based on a simple microscopic Hamiltonian, with only four parameters that are wellconstrained by bulk experimental data. With one additional parameter for the range of interaction, this model yields a computationally efficient freeenergy functional for inhomogeneous water, which captures shortranged correlations, cavitation energies, and, with suitable longrange corrections, the nonlinear dielectric response of water, making it an excellent candidate for the studies of mesoscale water and for use in ab initio solvation methods.

Fourbody interaction energy for compressed solid krypton from quantum theory
View Description Hide DescriptionThe importance of the fourbody contribution in compressed solid krypton was first evaluated using the manybody expansion method and the coupled clustertheory with full single and double excitations plus perturbative treatment of triples. All different fouratom clusters existing in the first and secondnearest neighbor shells of facecentered cubic krypton were considered, and both selfconsistentfield HartreeFock and correlation parts of the fourbody interaction were accurately determined from the ambient conditions up to eightfold volume compression. We find that the fourbody interaction energy is negative at compression ratio lower than 2, where the dispersive forces play a dominant role. With increasing the compression, the fourbody contribution becomes repulsive and significantly cancels the oversoftening effects of the threebody potential. The obtained equation of state (EOS) was compared with the experiments and the densityfunctional theory calculations. It shows that combination of the fourbody effects with two and threebody interactions leads to an excellent agreement with EOS measurements throughout the whole experimental range 0–130 GPa, and extends the prediction to 300 GPa.

Comparison of some dispersioncorrected and traditional functionals as applied to peptides and conformations of cyclohexane derivatives
View Description Hide DescriptionWe compare the energetic and structural properties of fully optimized αhelical and antiparallel βsheet polyalanines and the energetic differences between axial and equatorial conformations of three cyclohexane derivatives (methyl, fluoro, and chloro) as calculated using several functionals designed to treat dispersion (B97D, ωB97xD, M06, M06L, and M062X) with other traditional functionals not specifically parametrized to treat dispersion (B3LYP, X3LYP, and PBE1PBE) and with experimental results. Those functionals developed to treat dispersion significantly overestimate interaction enthalpies of folding for the αhelix and predict unreasonable structures that contain Ramachandran ϕ and ψ and C = O…N Hbonding angles that are out of the bounds of databases compiled the βsheets. These structures are consistent with overestimation of the interaction energies. For the cyclohexanes, these functionals overestimate the stabilities of the axial conformation, especially when used with smaller basis sets. Their performance improves when the basis set is improved from D95** to augccpVTZ (which would not be possible with systems as large as the peptides).

Unrestricted HartreeFock based on the fragment molecular orbital method: Energy and its analytic gradient
View Description Hide DescriptionA consideration of the surrounding environment is necessary for a meaningful analysis of the reaction activity in large molecular systems. We propose an approach to perform unrestricted HartreeFock (UHF) calculations within the framework of the fragment molecular orbital (FMO) method (FMOUHF) to study large systems with unpaired electrons. Prior to an energy analysis one has to optimize geometry, which requires an accurate analytic energy gradient. We derive the FMOUHF energy and its analytic gradient and implement them into GAMESS. The performance of FMOUHF is evaluated for a solvated organic molecule and a solvated metal complex, as well as for the active part of a protein, in terms of energy, gradient, and geometry optimization.

An adaptive pseudospectral method for wave packet dynamics
View Description Hide DescriptionWe solve the timedependent Schrödinger equation for molecular dynamics using a pseudospectral method with global, exponentially decaying, Hagedorn basis functions. The approximation properties of the Hagedorn basis depend strongly on the scaling of the spatial coordinates. Using results from control theory we develop a timedependent scaling which adaptively matches the basis to the wave packet. The method requires no knowledge of the Hessian of the potential. The viability of the method is demonstrated on a model for the photodissociation of IBr, using a Fourier basis in the bound state and Hagedorn bases in the dissociative states. Using the new approach to adapting the basis we are able to solve the problem with less than half the number of basis functions otherwise necessary. We also present calculations on a twodimensional model of CO_{2} where the new method considerably reduces the required number of basis functions compared to the Fourier pseudospectral method.

Precise description of single and double ionization of hydrogen molecule in intense laser pulses
View Description Hide DescriptionA new simulation box setup is introduced for the precise description of the wavepacket evolution of two electronic systems in intense laser pulses. In this box, the regions of the hydrogen molecule H_{2}, and singly and doubly ionized species, and , are well discernible and their timedependent populations are calculated at different laser field intensities. In addition, some new regions are introduced and characterized as quasidouble ionization and their timedependencies on the laser field intensity are calculated and analyzed. The adopted simulation box setup is special in that it assures proper evaluation of the second ionization. In this study, the dynamics of the electrons and nuclei of the hydrogen molecule are separated based on the adiabatic approximation. The timedependent Schrödinger and Newton equations are solved simultaneously for the electrons and the nuclei, respectively. Laser pulses of 390 nm wavelength at four different intensities (i.e., 1 × 10^{14}, 5 × 10^{14}, 1 × 10^{15}, and 5 × 10^{15} W cm^{−2}) are used in these simulations. Details of the central H_{2} region are also presented and discussed. This region is divided into four subregions related to the ionic state H^{+}H^{−} and covalent (natural) state HH. The effect of the motion of nuclei on the enhanced ionization is discussed. Finally, some different timedependent properties are calculated, their dependencies on the intensity of the laser pulse are studied, and their correlations with the populations of different regions are analyzed.

Timedependent quantum transport: An efficient method based on LiouvillevonNeumann equation for singleelectron density matrix
View Description Hide DescriptionBasing on our hierarchical equations of motion for timedependent quantum transport[X. Zheng, G. H. Chen, Y. Mo, S. K. Koo, H. Tian, C. Y. Yam, and Y. J. Yan, J. Chem. Phys.133, 114101 (2010)10.1063/1.3475566], we develop an efficient and accurate numerical algorithm to solve the LiouvillevonNeumann equation. We solve the realtime evolution of the reduced singleelectron density matrix at the tightbinding level. Calculations are carried out to simulate the transient current through a linear chain of atoms, with each represented by a single orbital. The selfenergy matrix is expanded in terms of multiple Lorentzian functions, and the Fermi distribution function is evaluated via the Padè spectrum decomposition. This LorentzianPadè decomposition scheme is employed to simulate the transient current. With sufficient Lorentzian functions used to fit the selfenergy matrices, we show that the lead spectral function and the dynamics response can be treated accurately. Compared to the conventional master equation approaches, our method is much more efficient as the computational time scales cubically with the system size and linearly with the simulation time. As a result, the simulations of the transient currents through systems containing up to one hundred of atoms have been carried out. As density functional theory is also an effective oneparticle theory, the LorentzianPadè decomposition scheme developed here can be generalized for firstprinciples simulation of realistic systems.

Nuclear motion effects on the density matrix of crystals: An ab initio Monte Carlo harmonic approach
View Description Hide DescriptionIn the frame of the BornOppenheimer approximation, nuclear motions in crystals can be simulated rather accurately using a harmonic model. In turn, the electronic firstorder density matrix (DM) can be expressed as the statistically weighted average over all its determinations each resulting from an instantaneous nuclear configuration. This model has been implemented in a computational scheme which adopts an ab initio oneelectron (HartreeFock or KohnSham) Hamiltonian in the CRYSTAL program. After selecting a supercell of reasonable size and solving the corresponding vibrational problem in the harmonic approximation, a Metropolis algorithm is adopted for generating a sample of nuclear configurations which reflects their probability distribution at a given temperature. For each configuration in the sample the “instantaneous” DM is calculated, and its contribution to the observables of interest is extracted. Translational and point symmetry of the crystal as reflected in its average DM are fully exploited. The influence of zeropoint and thermal motion of nuclei on such important firstorder observables as xray structure factors and Compton profiles can thus be estimated.

Linearresponse theory for Mukherjee's multireference coupledcluster method: Static and dynamic polarizabilities
View Description Hide DescriptionThe formalism of response theory is applied to derive expressions for static and dynamic polarizabilities within the statespecific multireference coupledcluster theory suggested by Mukherjee and coworkers (MkMRCC) [J. Chem. Phys.110, 6171 (1998)]. We show that the redundancy problem inherent to MkMRCC theory gives rise to spurious poles in the MkMRCC response functions, which hampers the reliable calculation of dynamic polarizabilities. Furthermore, we demonstrate that in the case of a symmetrybreaking perturbation a working response theory is obtained only if certain internal excitations are included in the responses of the cluster amplitudes. Exemplary calculations within the singles and doubles approximation (MkMRCCSD) are carried out on aryne compounds to illustrate the impact of a multireference ansatz on the polarizability.

Linearresponse theory for Mukherjee's multireference coupledcluster method: Excitation energies
View Description Hide DescriptionThe recently presented linearresponse function for Mukherjee's multireference coupledcluster method (MkMRCC) [T.C. Jagau and J. Gauss, J. Chem. Phys.137, 044115 (2012)]10.1063/1.4734308 is employed to determine vertical excitation energies within the singles and doubles approximation (MkMRCCSDLR) for ozone as well as for obenzyne, mbenzyne, and pbenzyne, which display increasing multireference character in their ground states. In order to assess the impact of a multireference groundstatewavefunction on excitation energies, we compare all our results to those obtained at the singlereference coupledcluster level of theory within the singles and doubles as well as within the singles, doubles, and triples approximation. Special attention is paid to the artificial splitting of certain excited states which arises from the redundancy intrinsic to MkMRCC theory and hinders the straightforward application of the MkMRCCLR method.

Fluctuating hydrodynamics for multiscale modeling and simulation: Energy and heat transfer in molecular fluids
View Description Hide DescriptionThis work illustrates that fluctuatinghydrodynamics (FHD) simulations can be used to capture the thermodynamic and hydrodynamic responses of molecular fluids at the nanoscale, including those associated with energy and heat transfer. Using allatom molecular dynamics (MD) trajectories as the reference data, the atomistic coordinates of each snapshot are mapped onto mass, momentum, and energy density fields on Eulerian grids to generate a corresponding field trajectory. The molecular lengthscale associated with finite molecule size is explicitly imposed during this coarsegraining by requiring that the variances of density fields scale inversely with the grid volume. From the fluctuations of field variables, the response functions and transport coefficients encoded in the allatom MD trajectory are computed. By using the extracted fluid properties in FHD simulations, we show that the fluctuations and relaxation of hydrodynamic fields quantitatively match with those observed in the reference allatom MD trajectory, hence establishing compatibility between the atomistic and field representations. We also show that inclusion of energy transfer in the FHD equations can more accurately capture the thermodynamic and hydrodynamic responses of molecular fluids. The results indicate that the proposed MDtoFHD mapping with explicit consideration of finite molecule size provides a robust framework for coarsegraining the solution phase of complex molecular systems.

Free energy landscapes for the thermodynamic understanding of adsorptioninduced deformations and structural transitions in porous materials
View Description Hide DescriptionSoft porous crystals are flexible metalorganic frameworks that respond to physical stimuli such as temperature, pressure, and gas adsorption by large changes in their structure and unit cell volume. While they have attracted a lot of interest, molecular simulation methods that directly couple adsorption and large structural deformations in an efficient manner are still lacking. We propose here a new Monte Carlo simulation method based on nonBoltzmann sampling in (guest loading, volume) space using the Wang–Landau algorithm, and show that it can be used to fully characterize the adsorptionproperties and the material's response to adsorption at thermodynamic equilibrium. We showcase this new method on a simple model of the MIL53 family of breathing materials, demonstrating its potential and contrasting it with the pitfalls of direct, Boltzmann simulations. We furthermore propose an explanation for the hysteretic nature of adsorption in terms of free energy barriers between the two metastable host phases.

Algebraicdiagrammatic construction polarization propagator approach to indirect nuclear spin–spin coupling constants
View Description Hide DescriptionA new polarization propagator approach to indirect nuclear spin–spin coupling constantans is formulated within the framework of the algebraicdiagrammatic construction (ADC) approximation and implemented at the level of the strict secondorder approximation scheme, ADC(2). The ADC approach possesses transparent computational procedure operating with Hermitian matrix quantities defined with respect to physical excitations. It is sizeconsistent and easily extendable to higher orders via the hierarchy of available ADC approximation schemes. The ADC(2) method is tested in the first applications to HF, N_{2}, CO, H_{2}O, HCN, NH_{3}, CH_{4}, C_{2}H_{2}, PH_{3}, SiH_{4}, CH_{3}F, and C_{2}H_{4}. The calculated indirect nuclear spin–spin coupling constants are in good agreement with the experimental data and results of the secondorder polarization propagator approximation method. The computational effort of the ADC(2) scheme scales as n ^{5} with respect to the number of molecular orbitals n, which makes this method promising for applications to larger molecules.