December 8 - 12, 2008
Continuum electrostatics methods have become increasingly popular due to their ability to provide approximate descriptions of solvation energies and forces without the expensive sampling required by all-atom solvent models. In particular, the Poisson–Boltzmann equation (PBE) provides electrostatic potentials, solvation energies, and forces by modeling the solvent as a featureless dielectric material and the mobile ions as a continuous distribution of charge. Polar solvation forces and energies obtained from the PBE are often supplemented with simple solvent-accessible surface area (SASA) models of nonpolar solvation. However, while polar and nonpolar continuum models have been assessed on their ability to reproduce global properties, such as solvation free energies, their ability to provide accurate representations of local solvation properties such as forces has not previously been adequately studied. We have developed efficient software for describing polar biomolecular solvation by solving the Poisson-Boltzmann equation using multigrid and adaptive finite element methods. Additionally, we have implemented new models to describe nonpolar solvation phenomena. These models have been used to study solvation forces for protein, RNA, and alkane systems. In particular, we have performed comparisons of continuum and all-atom representations of solvation forces for these very different molecular systems in order to assess the performance of continuum models in the presence of widely varying charge densities. The results of these comparisons show that current implementations of the PBE are capable of generating polar solvation forces that correlate well with explicit solvent forces for protein systems but provide significantly less accurate representations of polar solvation forces for RNA systems. Conversely, SASA-based nonpolar forces are found to have no significant correlation with nonpolar explicit solvent forces for either protein or RNA molecules. Good correlation between explicit and continuum nonpolar forces is only obtained when area, volume, and attractive dispersion forces are included in the continuum model. We discuss the implications of these studies in the context of molecular simulation as well as the impact of this work on basic models for understanding experimental observations of biomolecular binding and folding.
Solvation makes major contributions to the energetics of protein folding. In an unfolded protein the free energy of solvating nonpolar side chains is unfavorable while solvating polar peptide groups is favorable. The classical model for the energetics of burying nonpolar side chains through folding is Kauzmann's 1959 proposal to use transfer data for model hydrocarbons from water to an organic solvent. The simple picture is that 50 square angstroms of water-accessible surface area (ASA) per average side chain is buried via folding and -25 cal mol-1 is gained per square angstrom or -1.25 kcal mol-1 per residue. Today this model is regarded as seriously over-simplified and side chain burial is modeled by a 2-step process. (1) First remove the nonpolar side chain from water, using liquid to gas phase transfer data; (2) pack the folded protein side chains using the Lennard-Jones potential and protein structural coordinates to find the packing energy. A major unsolved problem is the highly approximate proportionality between solvation free energy and ASA. Model compound data give a huge value for the solvation enthalpy of the peptide group (-14.2 kcal/mol, Makhatadze & Privalov, 1993), more than 10 times larger than the free energy change per residue for burying nonpolar side chains. Its size shows that more work is badly needed on determining accurate energetics for peptide solvation and forming peptide H-bonds. Recent work shows that the principle of group additivity is not valid for the polar peptide group.
Two classes of protein groups in an unfolded protein interact strongly with water and make major contributions to the energetics of folding. Nonpolar side chains interact unfavorably with water and become buried in the protein interior during folding. Polar peptide groups interact favorably and very strongly with water, and as the protein folds the peptide groups become desolvated and form hydrogen bonds (H-bonds) with each other.
Two models are used currently to represent desolvation of nonpolar side chains during folding. Kauzmann’s model (1959) treats the protein interior as an organic liquid and uses liquid-liquid transfer data to estimate the enerrgetics. The packing-desolvation model (Privalov and Gill, 1988) treats the protein interior as a semi-crystalline solid and uses gas to liquid transfer data to estimate desolvation while using the Lennard-Jones potential and the protein structural coordinates to estimate the packing energetics.
Solvation energetics of the peptide group are based on calorimetric data for the solvation enthalpy of dipeptide analogs and mono-amides together with electrostatic calculations of solvation energetics for longer peptides. These results indicate that, contrary to the longstanding assumption that group additivity is valid, it is in fact entirely invalid for the polar peptide groups. These new results may explain the 2-fold discrepancy between the 1995 simulations of folding energetics by Lazaridis, Archontis and Karplus versus the calorimetric results of Makhatadze and Privalov.
Implicit-solvent models commonly decompose the molecular solvation free energy into a sum of electrostatic and non-polar terms, with separate theoretical and numerical procedures employed for each component. In this talk I will describe a new method for estimating the electrostatic free energy of solvation, called BIBEE (boundary-integral-based electrostatics estimation). The BIBEE method rapidly approximates the solution of the apparent-surface-charge (ASC) or polarizable-continuum-model (PCM) integral equation formulation of the mixed-dielectric continuum electrostatics problem that is often solved using finite-element or finite-difference methods. The method draws inspiration from the surface-generalized Born model, but avoids the unphysical interpolations inherent to generalized-Born method. As a result, BIBEE methods more accurately reproduce solvent-screened interactions between charges. Furthermore, it is straightforward to show that variants of the BIBEE method can generate rigorous upper and lower bounds for the actual electrostatic free energy of solvation.
Calcium channels conduct Na ions in the absence of Ca, but they selectively
conduct Ca ions when Ca ions are present at physiological concentrations. In
an experiment when Ca is added to NaCl gradually, even a micromolar amount of
Ca ions effectively blocks Na current. In our model of the selectivity filter
of Ca channels, the terminal groups of the side chains of amino acids—four
glutamates--in the selectivity filter are represented as mobile ions that are
restricted so they move inside the filter These structural ions form a
liquid-like self-adjusting environment for the passing ions so that the
system assumes minimum free energy. They also fill part of the pore so the
counterions have to compete for space in the crowded selectivity filter
(charge/space competition (CSC) mechanism). In this picture electrostatic
attraction and repulsive entropic excluded volume effects compete with each
other to determine which ions can enter the selectivity filter. We argue that
this competition is crucial in explaining the selectivity mechanism of Ca
channels. We show grand canonical Monte Carlo simulation results for
competition between ions of different valence and diameter. We couple our
Monte Carlo simulations to the integrated Nernst-Planck equation to compute
current from equilibrium profiles. Our results are in good agreement with
We draw a comparison between the quantum density functional
theory for electronic
structure calculations and the ''classical'' density
functional theory of molecular
liquids and we show how, borrowing ideas and techniques from
classical DFT can be used as a useful chemist's tool to
provide, at a microscopic
level, the solvation properties of complex molecules in polar
includes the determination of absolute solvation
free-energies, as well as
three-dimensional microscopic solvation structures. The
proposed strategy is
as follows: we first compute the homogeneous-fluid, position
direct correlation function, the c-function. To this end, we
carry out extensive MD
simulations of the pure solvent, and compute this way the
position and angle-dependent
pair distribution function (the h-function). We invert
subsequently the so-called
Ornstein-Zernike equation to go from the h-function to the
c-function. This direct
correlation can then be used as the definition of the
free-energy in the expression of the free-energy functional. In
the presence of a given
molecular solute, which provides the external potential, this
functional can be minimized
with respect to the position and angle-dependent density,
using a 3D cartesian grid
for positions and a Gauss-Legendre angular grid for
orientations, to obtain, at the minimum,
the absolute solvation free-energy of the solute and the
equilibrium solvent density
profile around it. The DFT results can be compared to direct MD
simulations of the solute/solvent
system or experimental data The procedure is shown to be
efficient and accurate for
polar solvents such as acetonitrile.
 Rosa Ramirez, Ralph Gebauer, Michel Mareschal, and Daniel
Borgis, Phys. Rev. E, 66, 031206 (2002).
 Rosa Ramirez and Daniel Borgis, J. Phys. Chem B 109, 6754
 Rosa Ramirez, Michel Mareschal, and Daniel Borgis, Chem.
Phys. 319, 261 (2005).
 Lionel Gendre, PhD Thesis, Université d'Evry, July 2008.
Manuscript in preparation.
The increasing efforts devoted to study linear and non linear optical properties (NLO) of solvated molecules have followed the success of modern quantum mechanical (QM) methods in the forecast of the NLO properties of isolated systems. The approaches adopted uses the same QM methodology developed for the isolated systems with the additional introduction of the features due to the solute-solvent interactions. However, even when all the solvent effects are included in the solvation model the computed NLO quantities are still microscopic in nature and cannot be directly compared with their macroscopic manifestation, i.e. the macroscopic scusceptibility.
In recent years we have introduced an effective framework to treat local field effects in linear and non-linear NLO properties in solution within the quantum mechanical Polarizable Continuum Model [1,2]. In this framework effective polarizabilities are defined to include the effect due to the difference between the local field acting on the solute molecules and the macroscopic field in the medium (Maxwell field).
In this poster we will focus on the vibrational components of the effective polarizabilities within the quantum mechanical Polarizable Continuum Model framework. In particular, we will give a detailed overview of the method for the calculation of the vibrational effective polarizabilities in the double harmonic approximation  . As will be show, the double harmonic procedure can be reformulated within the PCM so as to obtain the effective vibrational polarizabilities in terms of the derivatives with respect to normal coordinates of effective electronic molecular properties . We will also discuss the connection between the PCM formulation of the vibrational effective polarizabilities and the parallel approach given within the framework of semiclassical Onsager-Wortmann-Bishop method .
 R. Cammi, B. Mennucci, and J. Tomasi, J. Phys. Chem. A, 102(1998)870.
 R. Cammi, B. Mennucci, and J. Tomasi, J. Phys. Chem. A, 104(2000)4690.
 J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 105(2005)2999.
 R. Cammi, B. Mennucci, and J. Tomasi, in M.G. Papadopulos (ed.), Nonlinear Optical Response of Molecules Solids and Liquids: Methods and Applications, Research Signpost, Kerala, India, 2003, p.113.
 R. Wortmann and D. Bishop, J. Chem. Phys., 108(1998)1001.
Quantum chemistry computational tools coupled with a continuum description of the solvent offer an effective approach to study molecular properties and processes in condensed phase due to modest increase of the computational effort with respect a QM calculation of isolated molecules .
However, this coupling is by no meas free of problems. The QM continuum models are characterized by the use of a non-linear Hamiltonian, being the solute Hamiltonian term representing by the solute-solvent interaction dependent on the wave-function of the solute it selves. As a consequence all the conventional Quantum Chemical models to be able to describe in a coherent way the effects of the solute-solvent interaction have to be reconsidered from their basic principles and suitably modified to take into account of this non linearity. The neglecting of this effects may lead to a systematic errors in the numerical outcome of the QM solvation models.
In this talk we will give an overview of the systematic work done within the Polarizable Continuum Model group to introduce the non-linearity effects into the Quantum Chemical methods. The emphasis will be given on the basic principles and problems involved. In particular we will focus on the definition of a new basic energetic functional  , different, from that used for isolated molecules, which has to be taken as a starting point to develop QM methods for continuum solvation models.
A wide sample of QM-solvation methods will be also discussed, ranging from the ab-initio time independent methods (Hartree-Fock, MC-SCF, Moeller Plesset) to the corresponding time-dependent response function formalisms, and from the the Density Functional Theory to the Time dependent DFT approach to the excited states of solvated chromophores.
Final remarks will regards some still open issues in the coupling of QM methods with continuum solvation models.
 Tomasi J., B. Mennucci, Cammi, R. "Quantum Mechanical Continuum Solvations Models", Chem. Rev., 2005, 105, 2999.
 Cammi, R., Tomasi, J. J. Comp. Chem. 16, 1449, 1994
Understanding the electrical response of electrochemical convertors, such as fuel cells or batteries, involves elucidating the effect of the macroscopic voltage on the microscopic charge distribution at the electrode-electrolyte interface.
I will present a quantum/classical model which couples a quantum molecular description of the electrode-electrolyte interface with a polarizable-continuum representation of the long-range effects of the ionic solvent. I will mainly focus on the mathematical and numerical aspects. In the last part of my talk, I will present numerical calculations of the vibrational Stark effect for chemisorbed CO, which demonstrate the efficiency of this approach. This is a joint work with I. Dabo, Y. Li and N. Marzari.
Several phenomena take place at the interface between liquid water and the gold (111) surface. For example, it is a widespread system for controlled electrochemical investigations, for the study of molecular conduction and of the interaction between proteins and the solid surfaces in water. Despite this great and interdisciplinary importance, the wettability properties of the gold surface are microscopically not well understood. For example, the contact angle between water and clean gold is zero, thus the gold surface is considered hydrophilic . On the other hand, computational DFT studies on single molecules (or model monolayers) on the gold surface pointed out that the gold-water interaction is very small . The term hydrophobic is often used in this context to classify Au(111).
To shed light on the microscopic structure of the water-Au(111) interface, we performed classical molecular dynamics simulations of water on a gold (111) surface. These simulations are based on a force field for gold-water and gold-protein interactions that we have recently developed . It is based on available experimental results and quantum mechanical (DFT and MP2) calculations. Gold polarization effects are also taken into account.
In this talk I will discuss the developed force field and the results of our simulations for a neutral gold surface, with particular emphasis on the microscopic nature of the water-gold interface
 K.W. Bewing, and W. A. Zisman, J. Phys. Chem 69, 4238 (1965).
 A. Michaelides et al., Phys. Rev. Lett. 90, 216102 (2003).
 F. Iori, and S. Corni J. Comp. Chem. 29, 1656 (2008); F. Iori et al. J. Comp. Chem. in press.
Understanding solvent effects, which cause line broadening and screen electronic interactions, is fundamentally important because of its central role in the control of EET dynamics. Recent work has furthermore shown that simple models for solvation may not be sufficient to explain effects that go beyond Förster theory, such as coherent contribution to energy transfer. Here results obtained from novel quantum-mechanical methodologies will be reported as a starting point for an atomistic description of the interplay between solvation and electronic energy transfer. First, we will present results obtained with a quantum-mechanical method that accounts for the actual shape of the molecules inside the dielectric environment through the Polarizable Continuum Model (PCM), thus overcoming the point-dipole assumption of Förster theory. This model predicts an exponential distance-dependent attenuation of the solvent screening in chromophore pairs taken from light harvesting antenna proteins, thus indicating a substantial underestimation of EET rates by Förster theory at separations less than about 20 Å. As a further step towards a realistic model of solvation in EET, we will present a combined quantum mechanics/molecular mechanics (QM/MM) method that adopts an atomistic description of the solvent, thus going beyond the continuum dielectric approximation.
 (a) Jang, S.; Newton, M. D.; Silbey, R. J. Phys. Rev. Lett. 2004, 92, 218301. (b) Engel, G. S.; Calhoun, T. R.; Read, E. L.; Ahn, T.-K.; Mancal, T.; Cheng, Y.-C.; Blankenship, R. E.; Fleming, G. R. Nature 2007, 446, 782.
 (a) Scholes, G. D.; Curutchet, C.; Mennucci, B.; Cammi, R.; Tomasi, J. J. Phys. Chem. B 2007, 111, 6978. (b) Curutchet, C.; Scholes, G. D.; Mennucci, B.; Cammi, R. J. Phys. Chem. B 2007, 111, 13253.
Solvation free energy determines the thermostability of molecules in solution.
Accurate prediction of solvation free energy from atomistic simulations
have great impact on the understanding of many biologically important
processes. What makes the simulation of the solvation free energy difficult
is that it consists of competing effects from opposing interactions
at a molecular level. To accelerate convergence, separations of the competing
forces are necessary. Traditionally formulation of solvation free energy
separates van der Waals and electrostatics. We find that further separation
of the attractive and repulsive forces within the van der Waals interactions
is also necessary. Test calculations with constant pressure molecular
dynamics (MD) in water show that the seemingly small nonpolar
solvation free energy
is made up of relatively large unfavorable contribution from repulsive forces
and a canceling favorable contribution from the attractive forces. Our
results show good agreement with experimental solvation free energies and
those of other computational studies. Alternative to the constant pressure
simulation with a fluctuating volume, the same results can be generated with
fixed volume but fluctuating number of water molecules --- grand canonical
ensemble. With the grand canonical Monte Carlo (GCMC) moves on solvent water,
we are able to obtain accurate solvation free energy as well.
The GCMC/MD method is shown to yielded accurate results for ligand binding in
cytochrome p450, in which case water displacement from a confined
binding pocket can not be account for in typical MD simulation time.
Water at the surface of proteins plays a crucial biological
role. The study of hydration is therefore fundamental for
achieving a complete description of key factors determining the
protein biology. In this framework, Molecular Dynamics
simulations have provided precious tools for elucidating the
structure and dynamics of waters in the protein hydration
shell. We employed Molecular Dynamics to address on the
hydration properties of the Prion Protein, whose misfolding and
aggregation is associated to Transmissible Spongiform
Encephalopathies. The investigations focused on different
states along a likely aggregation pathway specifically
analyzing native state, misfolded state and amyloid-like state.
The presentation will discuss on the influence that water
exerts on protein structural stability [1, 2], intermolecular
interactions , misfolding  and self-assembly [5, 6].
1. PNAS (2005) 102:7535-7540.
2. FEBS Letters (2006) 580:2488-2494.
3. Biophysical Journal (2006) 90:3052-61.
4. Biophysical Journal (2007) 93:1284-1292.
5. Proteins (2008) 70:863-872.
6. BBRC (2008) 366:800-806.
Hydration shells of normal proteins display both structured and bulk-
like waters. Isomers with considerable bulk-like hydration tend to
aggregate. In my talk, I will present both theoretical and experimental
data showing that different morphological states of aggregated isomers
differ by hydration distribution profiles and water magnetic resonance
(MR) signals. The results help explain the MR contrast patterns of
amyloids, a subject of long controversy, and suggest a new approach for
identifying unusual protein aggregation related to disease. As an
example, I will present MRI data and cell toxicity measurements revealing
the relationship between the structural conformations of several
amyloidogenic peptides and the mechanisms of cellular dysfunction.
Selectivity is one of the most important properties of living systems. One of the founders of molecular biology (Nobel Laureate Aaron Klug) recently said (with some hyperbole I suspect) "There is only one word that matters in biology and that is specificity."
My collaborators and I study selectivity in ion channels. Ion channels are proteins with a hole down their middle that are the (nano nearly pico)valves of life. Ion channels control an enormous range of biological function in health and disease. A large amount of data is available about selectivity in many channels. Selectivity in ion channels occurs without structural change of the channel protein (on the biological time scale of 10-5 sec or longer) and does not involve changes in covalent bonds (i.e., changes in shape of electron orbitals). Selectivity in channels involves only electrodiffusion—usually of charged hard spheres. Thus, physical analysis of selectivity in ion channels is easier than analysis of specificity in enzymes or many other proteins while being at least as important biologically.
A simple pillbox model with two adjustable parameters accounts for the selectivity of both DEEA Ca channels (Aspartate Glutamate Glutamate Alanine) and DEKA Na channels (Aspartate Glutamate Lysine Alanine) in many ionic solutions of different composition and concentration. The predicted properties of the Na and Ca channels are very different even though 'Pauling' crystal radii are used for ions and all parameters are the same for both channels in all solutions. Only the side chains are different in the model of the Ca and Na channels. No information from crystal structures is used in the model. Side chains of the channel protein are grossly approximated as spheres.
How can such a simple model give such powerful results when chemical intuition says that selectivity depends on the precise relation of ions and side chains? We use Monte Carlo simulations of this model that determine the most stable-lowest free energy-structure of the ions and side chains. Structure is the computed consequence of the forces in this model and so is different in different conditions. The relationship of ions and side chains vary with ionic solution and are very different in simulations of the Na and Ca channels. Selectivity is a consequence of the 'induced fit' of side chains to ions and depends on the flexibility (entropy) of the side chains as well as their location. The induced fit depends on the concentrations of ions in the surrounding solutions in a complex way. Thus, calculations in a single set of conditions are of limited use. In particular, calculations of 'free energy of binding' in infinitely dilute or ideal solutions are not likely to give useful estimates of binding in physiological solutions. Physiological solutions are typically ~ 200 mM, far from dilute.
The self-organized induced-fit model captures the relation of side chains and ions well enough to account for selectivity of both Na channels and Ca channels in the wide range of conditions measured in experiments, even though the components of the model are grossly oversimplified. Perhaps the simplified model works because the structures in both the model and the real channel are the most stable, self-organized, and at their free energy minimum, different in different conditions.
It seems that an important biological function can be understood by an oversimplified model if the model calculates the 'most stable' structure as it changes from solution to solution, and mutation to mutation.
The study of chromophores located at surfaces and interfaces is
important both in material science and in biological chemistry: coated
materials, sensors, nanoparticles, cell membranes, micelles are
systems where the interface is not a mere boundary between different
bulk regions but the central part of the investigated
system. Chromophores are often employed to probe and investigate the
surface by making use of surface-specific techniques as second
harmonic generation (SHG) and sum-frequency generation (SFG). It is
therefore important to develop the theoretical tools necessary to
model chromophores in such a peculiar environment in order to
interpret and understand the spectroscopic findings.
Continuum solvation models have in recent years been extended to
surfaces and interfaces. A brief account of the challenges faced and
the theoretical developments needed to achieve the goal will here be
given. Recent results making use of the developed theoretical models
will also be presented.
Traditionally, molecular dynamics and Monte Carlo simulations of
condensed phase systems including biopolymers are carried out using molecular mechanics or force fields that describe intermolecular interactions. In fact, the formalism of these force fields for biomolecular simulatiosn was established in the 1960s. Although the computational accuracy has increased
enormously through parameterization, little has changed in
the functional terms used in the force fields. In this talk, I will describe an explicit polarization
(X-Pol) method based on quantum mechanics for developing force fields in condensed phase simulations. The theory, algorithm and application of this approach will be illustrated and its feasibility is demonstrated.
Water molecules are ubiquitous in living organisms and have therefore been
viewed more as an environment for biomolecules rather than as a collection
of interacting molecules. Water molecules make up an integral part of
protein structures, while assisting in catalysis, providing stability and
controlling the plasticity of binding sites. In order to realistically mimic
the environment of biomolecules, molecular dynamics simulations are
routinely done in explicit water. Unfortunately, most of the computational
resources go into computing the interactions between these water molecules.
Therefore, many implicit solvation models pursued over the years have only
viewed the presence of water as a continuum dielectric. However, critical
questions do arise about the inherent faster dynamics that are usually
obtained with implicit solvation models. We show that explicit water does
not only slow down protein dynamics by increasing the frictional drag, but
also by increasing the local energetic roughness of the energy landscape by
as much as 1.0 kcal/mol, an effect which is lacking in typical implicit
solvation models. The possible implications of this effect in catalysis will
also be discussed.
Few years ago, we have succeeded to “probe” water molecules bound in a cavity of a protein by means of the statistical mechanics of molecular liquids, or the RISM/3D-RISM theory. This is the first finding in the history of the statistical mechanics to show that the theory is applicable to such fluids in an extremely inhomogeneous field in atomic scale. On the other hand, the finding implies that we got a powerful machinery in our hand to clarify the “molecular recognition” which is the most fundamental process in living cells. It will become a “bridgehead” in our theoretical challenge on life phenomena.
In the succeeding applications of the method, we have focused on a variety of molecular processes in bio-systems, in which the molecular recognition plays an essential role: the selective ion-binding by protein, an enzymatic reaction, water channels (aquaporin), the preferential binding of inert gas by protein, and the pressure denaturation of protein.
In the talk, I will present our recent studies on the chemical processes concerned with the molecular recognition, focusing especially on the possibility of ion permeation through aquaporins. The talk will include also prospects of the theory to be extended to the temporal fluctuation of protein structure coupled with the solvent dynamics.
Poisson-Boltzmann (PB) electrostatics is a well established model in biophysics, however its application to the study of large scale biosystem dynamics such as the protein-protein encounter is still limited by the efficiency and memory constraints of existing numerical techniques. In this poster, we present an efficient and accurate scheme which incorporates recently developed novel numerical techniques to further enhance the present computational ability. In particular, a boundary integral equation approach is applied to discretize the linearized Poisson-Boltzmann (LPB) equation; the resulting integral formulas are well conditioned and are extended to systems with arbitrary number of biomolecules; the solution process is accelerated by the Krylov subspace methods and an adaptive new version of fast multipole method (FMM); and in addition to the full electrostatic interaction energy, the forces and torques are computed in the post-processing procedure using interpolation. The Adaptive Fast Multipole Poisson-Boltzmann (AFMPB) solver will be released under open source license agreement.
We will discuss a class of fast algorithms for linear and
nonlinear integral equations. These are two-level algorithms
based on the classic Atkinson-Brakhage method from the 1970s.
We will present more efficient approach which uses a matrix-free
Newton-Krylov iteration on the coarse mesh and does the
fine-to-coarse intergrid transfer with an average. We will then
apply the approach to the Ornstein-Zernike (OZ) equations for atomic
fluids and some extensions of the OZ equations for molecular fluids.
Joint work with A.C. Pan and B. Roux.
Understanding the mechanism of conformational changes in macromolecules requires the knowledge of the intermediate states. A version of the string method, which uses multiple short dynamics trajectories to propagate the pathway, was recently developed by Pan et al. We use data from swarms of trajectories calculated at discrete points in phase space to interpolate the average displacement and variance of the system at arbitrary points. This is tested on model potentials using statistics from actual swarms of trajectories. We use the interpolated parameters to compute the Markovian propagators from one point on the transition path to the next. We use them to obtain a time-dependent action of a path, which can be optimized to produce the highest probability pathway. We describe the optimization protocol and demonstrate that in artificial flat potentials the existing string method cannot correct problems such as loops in the initial path, while the new method produces the correct pathway. The method will be tested by applying it to protein conformational transitions, such as the KcsA potassium channel, and comparing its performance to existing transition pathway methods.
Joint work with Alexander E. Kobryn (National Institute for Nanotechnology, National Research Council of Canada).
Motivated by the fundamental questions raised by the most recent experimental
achievements in nanofluidics, we propose the first-ever derivation and calculation
of the hydrodynamic slip length from the first principles of statistical mechanics,
based on a combination of linear response theory and equilibrium molecular theory
of solvation. The slip length obtained is independent of the type of flow and
is related to the fluid organization near the solid surface, as governed by the
solid-liquid interaction. In the wide range of shear rates and surface-liquid
interactions, the slip length is expressed in terms of the Green-Kubo-Nakano
relations as a function of the anisotropic inhomogeneous time correlation
function of density fluctuations of the liquid in contact with the surface.
The time dependence of the correlation function is factored out by treating it
in the hydrodynamic limit. And the spatially inhomogeneous two-body correlation
function is represented in the Kirkwood-like approximation as a product of the
three-dimensional density distributions of interaction sites of the liquid
near the surface and the site-site pair correlations of the bulk liquid.
The presented treatment generalizes the phenomenological definition of the
friction coefficient (as well as the slip length) to a tensor quantity,
which reflects an anisotropic nature of an ordered crystalline surface.
This enables theoretical prediction of friction forces acting aslant to
the liquid flow direction for such surfaces. We derive generic analytical
expressions for the liquid-surface friction coefficient (and slip length)
for an arbitrary surface-liquid interaction potential. We further illustrate
it by numerical calculations for the case of a laminar flow of acetonitrile,
benzene, ethanol, ethanolamine, dimethylsulfoxide, glycerol, methanol,
tert-butyl alcohol, and water, at ambient conditions in contact with the
(100) FCC surface of gold, copper and nickel modeled by using all-atom or
united-atom models for liquids and the Steele potential for crystalline surfaces.
The obtained values for slip length range from few to hundreds of nanometers
and are consistent with experimental measurements. We complement calculations
by obtaining pressure and temperature dependence of the slip length for water
in wide range of these thermodynamic parameters.
A salient feature of nanoscale objects and phenomena is that they are very different
from the atomic constituents as well as from those described by the conventional,
macroscopic laws of continuous media. The behavior of nanostructures can be changed
in a wide range, which constitutes the promise of nanoscience on control over
properties of materials. Examples are self-assembly of synthetic organic supramolecular
nanoarchitectures in solution, nanocatalysis, and electrochemistry in disordered
nanoporous electrodes. Modeling of nanosystems should operate at length scale from
one to hundreds nanometers and time scale up to microseconds and more, and yet
derive their properties from the chemical functionalities of the constituents.
Explicit molecular modeling of such nanosystems involves long-time description
of millions of molecules which is by far not feasible with ab initio methods,
very problematic for molecular simulations, and thus requires multiple-scale approaches.
Statistical-mechanical, molecular theory of solvation, a.k.a. three-dimensional
reference interaction site model (3D-RISM), predicts from the first principles
the molecular structure and thermodynamics of nanosystems, with proper account of
their chemical functionalities . The theory operates with 3D distributions of species
averaged over the statistical ensemble rather than with trajectories of individual
molecules, and yields the thermodynamics and volumetrics of solvation analytically.
It represents both long-range electrostatic and short-range steric features of the
solvation structure of chemical specificities, such as hydrogen bonding, hydrophobicity,
and electrochemical effects. The 3D-RISM theory can be regarded as a molecular extension
to the Poisson-Boltzmann solvation electrostatics. This presentation describes
the 3D-RISM methodology for the solvation structure and thermodynamics, and its
self-consistent combination with ab initio quantum chemical methods in multiscale
theory of electronic structure in solution. The methodology is illustrated with
applications to solvation of different molecules and chemical reactions in solution ,
self-assembly and conformational stability of synthetic organic supramolecular
nanoarchitectures , and supercapacitor and electrosorption cell devices with
nanoporous carbon electrodes .
 A. Kovalenko, Three-dimensional RISM theory for molecular liquids and solid-liquid interfaces,
in: Molecular Theory of Solvation, F. Hirata (Ed.) Series: Understanding Chemical Reactivity,
P. G. Mezey (Ed.), vol.24, (Kluwer Academic Publishers, Dordrecht, 2003, 360 p.) pp.169-275.
 S. Gusarov, T. Ziegler, A. Kovalenko, Self-Consistent Combination of the Three-Dimensional
RISM Theory of Molecular Solvation with Analytical Gradients and the Amsterdam Density Functional
Package, J. Phys. Chem. A, 110, 6083 (2006); D. Casanova, S. Gusarov, A. Kovalenko, T. Ziegler,
Evaluation of the SCF Combination of KS-DFT and 3D-RISM-KH; Solvation Effect on Conformational
Equilibria, Tautomerization Energies, and Activation Barriers, J. Chem. Theory Comput., 3, 458 (2007).
 R.S. Johnson, T. Yamazaki, A. Kovalenko, H. Fenniri, Molecular Basis for Water-Promoted
Supramolecular Chirality Inversion in Helical Rosette Nanotubes, J. Am. Chem. Soc., 129, 5735 (2007);
J.G. Moralez, J. Raez, T. Yamazaki, R.K. Motkuri, A. Kovalenko, H.Fenniri, Helical Rosette Nanotubes
with Tunable Stability and Hierarchy, J. Am. Chem. Soc. Communications, 127, 8307 (2005).
 A. Kovalenko, Molecular description of electrosorption in a nanoporous carbon electrode,
J. Comput. Theor. Nanosci., 1, 398 (2004); A. Tanimura, A. Kovalenko, F. Hirata, A molecular
theory of a double layer formed by aqueous electrolyte solution sorbed in a carbonized
polyvinylidene chloride nanoporous electrode, Chem. Phys. Lett., 378, 638 (2003).
Recent studies have questioned the consistency and applicability of the existing implicit-solvent models in which solvent accessible surfaces (SAS) or solvent excluded surfaces (SES) are pre-defined, and used for the calculation of solvation free energies. As an emerging concept and theory, the variational implicit solvation determines equilibrium solute-solvent interfaces and solvation free energies by minimizing a mean-field free-energy functional. This free energy couples the surface energy, dispersive interaction, and electrostatic contribution. It can also incorporate molecular mechanical interactions. A level-set method is developed to numerically relax an initial interface of high free energy to an equilibrium in the direction of steepest descent. Numerical results demonstrate the initial success of the coupling of the level-set method with the variational theory of solvation, particularly in capturing the hydrophobic interaction. Such interactions are crucial to biomolecular structuring and dynamics, but are not well accounted for by SAS/SES implicit-solvent models. This is joint work with Li-Tien Cheng, Zhongming Wang, Yang Xie, Jianwei Che, Joachim Dzubiella, and J. Andrew McCammon.
We present a novel Monte Carlo method for solving partial differential
equation (PDE) boundary-value problems (BVPs) involving the
Poisson-Boltzmann equation (PBE). Such BVPs arise in many situations
where the calculation of electrostatic properties of solvated large
molecules is desired. These techniques are very accurate, and are
fundamentally different than their deterministic counterparts. They are
based on using the Feynman-Kac formula to solve the Poisson and the
linearized PBE, as well as the introduction of a stochastic method to
enforce the boundary conditions on the surface of the molecule. In
addition, we present an algorithm for computing the electrostatic free
energy for all salt concentrations of interest, simultaneously. This,
coupled with the naturally parallel nature of Monte Carlo methods, makes
these techniques very appealing for studying solvation-based effects.
Joint work with Anna I. Krylov.
Experimental photoelectron angular distributions (PADs) contain
about the excited states of less stable compounds and their
dissociation mechanism, but
the interpretation of the results is difficult for molecules
with complex electronic
structure, e.g. NO dimer, CS2 photodissociation. Ionization
cross-sections and PADs can
be calculated from the corresponding electronic transition
dipole matrix elements.
The wavefunction of the leaving electron is given by Dyson
orbitals, which are
the overlap between an N electron molecular wavefunction and
the N-1/N+1 electron
wavefunction of the corresponding cation/anion. Our
implementation of the Dyson
orbitals calculation within the high level ab initio
EOM-IP/EA-CCSD method, allows the
calculation of these terms for the ionization of electronically
excited states and open-shell
species, when correlation effects become important. As a first
approximation, the states
of the ionized electron are described by plane waves expressed
in the bases of spherical
waves, |E, l, m>.
Currently, the calculation of photodetachment cross-sections
and anisotropies is
benchmarked on atoms, and small molecules. Calculated PADs for
photodissociation allow qualitative comparison with
experimental distributions. In the
case of larger, anisotropic molecules, the development of
better approximations for the
description of the ionized electron is necessary.
A macroscopic pattern consisting of liquid crystal layers in hydorgel
is self-organized. The pattern forms when potassium ion diffuses into
a kappa-carrageenan solution from one end of the glass capillary that
contains the kappa-carrageenan solution. Both the distance between two
adjacent liquidcrystalline layers (spacing, xn+1 −xn,) and thickness
of liquid crystalline layers (width, wn) depend linearly on the
distance from the diffusing end of the potassium ion xn. The time
prior to the formation of the nth liquid crystalline layer tn depends
linearly on x2 n. In addition, the spacing coefficient p
is inversely proportional to the concentration of potassium ions.
These results are in good agreement with the Liesegang henomenon. In
this system the kappa-carrageenan solution behaves as a supporting
medium for the spatial pattern, as well as the pattern forming
substance. The lower values of p, rather than the common Liesegang
pattern, in this system could be attributed to the large molecular
weight of the kappa-carrageenan–potassium complex. The pattern
consisted of discrete liquid crystal phases must be formed due to the
much larger diffusion constant of potassium ion (diffusant) than that
of kappa-carrageenan (product).
A modified Poisson-Boltzmann equation is developed from statistical mechanical considerations to describe the influence of the transmembrane potential on macromolecular systems . Using a Green's function formalism, the electrostatic free energy of a protein associated with the membrane is expressed as the sum of three terms: a contribution from the energy required to charge the system's capacitance, a contribution corresponding to the interaction of the protein charges with the membrane potential, and a contribution corresponding to a voltage-independent reaction field free energy. The membrane potential, which is due to the polarization interface, is calculated in the absence of the protein charges, whereas the reaction field is calculated in the absence of transmembrane potential. Then, a theoretical framework is elaborated to account for the effect of a transmembrane potential in computer simulations with explicit solvent. It is shown that a simulation with a constant external electric field applied in the direction normal to the membrane is equivalent to the influence of surrounding infinite baths maintained to a voltage difference via ion-exchanging electrodes connected to an electromotive force . It is also shown that the linearly-weighted displacement charge within the simulation system tracks the net flow of charge through the external circuit comprising the electromotive force and the electrodes. Using a statistical mechanical reduction of the degrees of freedom of the external system, three distinct theoretical routes are formulated and examined for the purpose of characterizing the free energy of a protein embedded in a membrane that is submitted to a voltage difference. The two methods are compared and applied to the voltage-gated potassium channel.
1. Roux, B. 1997. Influence of the membrane potential on the free energy of an intrinsic protein. Biophysical J. 73:2980-2989.
2. Roux, B. 2008. The membrane potential and its representation by a constant electric field in computer simulations. Biophys J. 95:4205-4216.
Descriptions of implicit solvent models like the Poisson-Boltzmann
equation express the electrostatic potential as the solution of an
elliptic PDE with delta function source terms. These equations are
particularly good fits with iterative solvers that use a fast N-body
solver as a preconditioner. Also, modern architectures favor such
algorithms due to their high ratio of floating-point operations
to memory references. Two types of N-body solvers can be distinguished:
hierarchical-clustering algorithms, such as the celebrated
fast multipole method, and kernel-splitting algorithms,
such as the popular particle--mesh Ewald method.
By formulating the problem as that of computing a matrix--vector
product, the basic structure of these algorithms is elucidated.
Additionally, evidence is presented indicating
that kernel-splitting algorithms are much to be preferred
for molecular simulations and that the virtually unknown multilevel
summation method of Brandt and Lubrecht is the best among these.
This method uses hierarchical interpolation of interaction potentials
on nested grids to calculate energies and forces in linear time
for both periodic and nonperiodic boundary conditions.
This is joint work with David Hardy. Its application to the
Poisson(-Boltzmann) equation is being pursued in a collaboration
with Stephen Bond.
Hybrid solvation models can be used in molecular dynamics simulations to extract the salient features of both explicit and implicit solvent methods, maintaining the atomic detail of solvent near regions of interest while removing such degrees of freedom elsewhere. We present such a hybrid solvation model that, in the spirit of previous methods, decomposes the simulation domain into explicit and bulk regions separated by an energetic barrier. The use of a Grand Canonical Monte Carlo scheme at the boundary is used to represent particle exchange with the bulk reservoir and extends the applicability of this model by allowing the use of a dynamic explicit domain. Current testing of this model is aimed at obtaining thermodynamic quantities within the explicit region with a focus on removing boundary artifacts for bulk liquid simulations.