July 23 - August 3, 2007
Joint work with Richard A. Bryce.
Alginate copolymers are a key component of the extracellular polymeric
substances (EPS) matrix of bacterial microorganisms such as P. aeruginosa, and
occur in high abundance in nature. Alginate heteropolysaccharides comprise
alternating blocks of alpha-(1->4)-linked L-guluronate and beta-(1->4)-linked
D-mannuronate. To understand the microscopic behaviour and interactions of
these flexible acidic sugars within the EPS matrix, a suitable molecular-level
model is required. We derive a molecular mechanical force field for the two
uronic acids, with validation against available experimental data. We then
construct a detailed study of the structure and dynamics of alginate chains. We
explore alginate models of increasing complexity, from disaccharides to single-
and double-stranded oligomer helices, employing the techniques of molecular
dynamics and replica-exchange molecular dynamics. The condensed phase behaviour
of these systems is discussed, including the role of counterions and the
implications for interaction with other constituents of the EPS matrix.
Same abstract as the 7/24 poster session.
The quantum many-electron problem is easy in principle (solve the N-electron Schrödinger equation) and hard in practice (because the cost of numerical methods typically grows exponentially with the number of variables). However, there are simplifying features. First, the dimensionality can be reduced because electronic Hamiltonians contain only 1-body and 2-body terms. (This leads to reduced density-matrix methods.) Second, the dimensionality can be reduced because electrons are identical particles: if you know everything about one electron, then you know everything about all of the electrons. (This leads to electron-propagator theory and density-functional theory.)
There is a “catch.” Reducing the number of dimensions leads to other problems associated with approximating the energy functional and/or associated with restricting the domain of the variational procedure. Two powerful techniques for resolving these difficulties are the Legendre transform and constrained-search formulations of density functional theory. This talk will discuss these formulations, and show how they can be extended to define "generalized" density-functional theories. I'll conclude with some of my favorite open problems in density-functional theory.
Exchange and correlation effects in electronic systems are rigorously related to
a two-electron function called the exchange-correlation "hole". Modelling of
the hole in real space is a powerful route to development and refinement of
exchange and correlation functionals in DFT. We have developed real-space
models of all correlation types of importance in chemical physics (dynamical,
nondynamical, and dispersion) and these will be reviewed.
In classical molecular dynamics, the accuracy of numerical methods
should be measured with respect to statistical averages, due to the
lack of a meaningful exact trajectory. In this talk I will survey
some results from backward error analysis and show how (under
certain assumptions) these results can be applied to compute
estimates of the error in averages from molecular dynamics
simulations. Results from several test problems will be explored
including examples from constant temperature molecular dynamics
using a Nosé thermostat.
Deformation and fracture are fundamental phenomena with major implications on the stability and reliability of machines, buildings and biological systems. All deformation processes begin with erratic motion of individual atoms around flaws or defects that quickly evolve into formation of macroscopic fractures as chemical bonds rupture rapidly, eventually compromising the integrity of the entire structure. However, most existing theories of fracture treat matter as a continuum, neglecting the existence of atoms or nanoscopic features. Clearly, such a description is questionable. Here we discuss an atomistic approach to describe such processes using ultra large-scale molecular dynamics (MD) simulation implemented supercomputers. MD provides unparalleled insight into the complex atomic-scale deformation processes, linking nano to macro, without relying on empirical input, since all atomic interaction parameters can be derived from fundamental quantum chemical theories. We demonstrate how MD can be used within a multi-scale simulation framework to predict the elastic and fracture properties of hierarchical protein materials, marvelous examples of structural designs that balance a multitude of tasks, representing some of the most sustainable material solutions that integrate structure and function across the scales. Breaking the material into its building blocks enables us to perform systematic studies of how microscopic design features influence the mechanical behavior at larger scales. We review studies of collagen – Nature’s super-glue, spider silk – a natural fiber that can reach the strength of a steel cable, as well as intermediate filaments – an important class of structural proteins responsible for the mechanical integrity of cells, which, if flawed, can cause serious diseases such as the rapid aging disease progeria. The common ground of these examples is the significance of the material properties at large deformation, its alteration under stress, presence of defects or the effect of variation of environmental conditions. Our studies elucidate intriguing material concepts that enable to balance strength, energy dissipation and robustness by selecting nanopatterned, hierarchical features.
We propose an experimental scheme to obtain the phase diagram
of the Hubbard model using cold atoms in optical lattices. The
scheme is based on measuring the total energy for a sequence of
trapping potentials with diﬀerent proﬁle and is independent of di-
mensionality. Its essential ingredient is a two-parameter scaling pro-
cedure that combines a variant of the familiar ﬁnite-size scaling with
a nontrivial additional ’ﬁnite-curvature scaling’ necessary to reach
the homogeneous limit. We illustrate the viability of the scheme
in the one-dimensional case, using simulations based on the Bethe
Ansatz local-density approximation as a substitute for experimental
data, and show that the ﬁlling corresponding to the Mott transition
can be determined with better than 3% accuracy.
Same abstract as the 7/24 poster session.
By means of thermodynamic limit arguments, the Kohn-Sham model can be extended to the case when the nuclear configuration is periodic. The resulting periodic Kohn-Sham model is very useful to understand the electronic structure of perfect crystals, but also to simulate the condensed phase by means of artificial periodic boundary conditions. In the latter case, the periodic cell may contain a very large number of nuclei and electrons.
In the first part of the lecture, I will present the Bloch theory of periodic Schrdinger operators, and recall how the band structure of their spectrum can be related to the insulating or conducting character of the underlying material. I will also present the concept of Wannier functions, which is of growing importance in various applications (computation of dielectric and magnetic properties, local correlation methods, electronic structure of crystals with defects, ...).
In the second part of the lecture, I will focus on the discretization of the periodic Kohn-Sham model in plane-wave (PW) basis sets.
Quantum Chemistry aims at understanding the properties of matter through the modelling of its behaviour at a subatomic scale, where matter is described as an assembly of nuclei and electrons.
At this scale, the equation that rules the interactions between these constitutive elements is the Schrdinger equation. It can be considered (except in few special cases notably those involving relativistic phenomena or nuclear reactions) as a universal model for at least three reasons. First it contains all the physical information of the system under consideration so that any of the properties of this system can be deduced in theory from the Schrdinger equation associated to it. Second, the Schrdinger equation does not involve any empirical parameter, except some fundamental constants of Physics (the Planck constant, the mass and charge of the electron, ...); it can thus be written for any kind of molecular system provided its chemical composition, in terms of natures of nuclei and number of electrons, is known. Third, this model enjoys remarkable predictive capabilities, as confirmed by comparisons with a large amount of experimental data of various types.
Unfortunately, the Schrdinger equation cannot be directly simulated, except for very small chemical systems. It indeed reads as a time-dependent 3(M+N)-dimensional partial differential equation, where M is the number of nuclei and N the number of the electrons in the system under consideration. On the basis of asymptotic and semiclassical limit arguments, it is however often possible to approximate the Schrdinger dynamics by the so-called Born-Oppenheimer dynamics, in which nuclei behave as classical point-like particles. The internuclei (or interatomic) potential can be computed ab initio, by solving the time-independent electronic Schrdinger equation.
The latter equation is a 3N-dimensional partial differential equation (it is in fact a spectral problem), for which several approximation methods are available. The main of them are the wavefunction methods and the Density Functional Theory (DFT). They lead in particular to the Hartree-Fock model and to the Kohn-Sham model, respectively.
In this introductory lecture, I will first present the electronic Schrdinger equation, and show how to derive from this equation the Hartree-Fock and Kohn-Sham models. Although obtained from totally different approaches, these models have similar mathematical structures. They read as constrained optimization problems, whose Euler-Lagrange equations are nonlinear eigenvalue problems. In the second part of the lecture, I will introduce the Linear Combination of Atomic Orbitals (LCAO) discretization method, and briefly discuss two important numerical issues: ensuring self-consistent field (SCF) convergence, and reducing computational scaling.
Quantum locality leads to matrix sparsity, which may be exploited to achieve
an O(N) complexity through the use of an inexact, sparse linear algebra and
iterative methods. The framework for static response is developed around perturbation
of the spectral projector. In the case of dynamic response, an orbital free,
Rayleigh Quotient Iteration (RQI) like method is outlined. Problems with error accumulation,
associated with the inexact sparse linear algebra, are discussed and it is pointed out that
in the static case, the accumulation of errors on iteration quench as idempotence is reached.
In the dynamic case, it is shown that the frequencies are variational with respect to matrix
One of the most challenging issues in materials physics is to
the properties of matter at the nanoscale. In this size
structural and electronic properties exist that resemble
atomic, nor solid state. These altered properties can have
technological implications. Theoretical methods to address
face formidable challenges. Nanoscale systems may contain
electrons and atoms, and often possess little symmetry. I
illustrate some recent advances in this area based on new
methods and apply these techniques to systems ranging from
clusters of a
few dozen atoms to quantum dots containing thousands of
Y. Zhou, Y. Saad, M.L. Tiago, and J.R. Chelikowsky:
Self-Consistent-Field Calculations via
Chebyshev-Filtered Subspace Acceleration,'' Phys. Rev. E 74,
M.L. Tiago, Y. Zhou, M.M.G. Alemany, Y. Saad, J.R.
Evolution of Magnetism in Iron from the Atom to the Bulk,''
Lett. 97, 147201 (2006).
M. Lopez del Puerto, M.L. Tiago, and J.R. Chelikowsky:
effects and optical properties of passivated CdSe clusters,''
Lett. 97, 096401 (2006).
A new simulation method is studied which is capable of calculating the
absolute entropy and free energy of a microstate. The method is based on
calculating the probability of a member in the statistical ensemble by
growing the configuration in a step-by-step manner and expressing it as
a product of conditional probabilities, each of which is calculated by
simulations at every step. This is a potentially exact method with the errors
arising only from inadequate sampling. The method is applied to a small protein loop in
implicit solvent and the results obtained agree quite well with those
obtained from other (approximate) methods. Unlike other
methods, all long range interactions are taken into account in the
Holonomic constraints as coarse graining of (large) molecules. The dynamical evolution of a system subjected to holonomic constraints (SHAKE and beyond): Lagrange & Hamilton eq.s of first type; the standard (unstable) solution; Shake , i.e. an evolution satisfying exactly the constraints at every time step.
The statistical measure and its associated evolution equation (Liouville equation) for mechanical systems with holonomic constraints: nonhamiltonian systems and their statistical behavior; the eq.s of motions of a mechanical systems subjected to constraints as a nonhamiltonian system; generalized Liouville eq. and statistical equilibrium probability.
The linear response theory formulation of a rate constan as a suitable time correlation function; its representation via the transition state theory (TST) term and the transmission coefficient. Calculation of the TST term via the blue-moon approach. i.e. impose a reaction coordinate as a constraint, compute its associated free energy and, by exponentiation, compute the associated probability. Use the unbiased constrained ensemble to generate a proper ensemble of reactive trajectories and compute the transmission coeffcient. The product of the TST term times the transmission coefficient provides the rate constant.
Hubbard U-corrected DFT functionals have been very successful in describing several strongly-correlated
systems for which "standard" approximations to DFT fail. Unfortunately no explicit expression exists for the
effective electronic interaction parameter (the Hubbard U) contained in the corrective ("+U") functional
and semiempirical estimates have been often used to determine its value. In this talk, after a general introduction
to the LDA+U method, I will present our linear response approach to the evaluation of the Hubbard U .
Within this approach the on-site electronic coupling is computed from the response of the considered system
to a shift in the potential acting on its correlated atomic states. Specifically, it is evaluated as the difference between
the inverse of the bare and fully interacting response matrices. The U we obtain thus corresponds to the effective
(atomically averaged) interaction between electrons that are located on the same site.
In this way the strength of the "+U" correction is consistently evaluated from the same DFT scheme we aim to correct;
the LDA+U is transformed in a completely ab-initio method with no need for any empirical evaluation of the effective coupling.
The results are also largely independent on the choice of the localized orbitals: the same occupation matrix that enters
the expression of the "+U" correction is consistently used to compute the effective interaction parameter.
With this approach we successfully studied the structural, electronic, chemical and electrochemical properties of several
transition metals compounds. Examples of applications will include minerals in the Earth's interior , cathode materials
for next-generation lithium-ion batteries  and catalysis reactions on molecules .
 M. Cococcioni and S. de Gironcoli, PRB (2005).
 F. Zhou, M. Cococcioni, A. C. Marianetti, D. Morgan and G. Ceder, PRB (2004).
 H. J. Kulik, M. Cococcioni, D. Scherlis and N. Marzari, PRL (2007).
Electronic correlation play a very important role in the
physical properties of many materials characterized by
very localized electronic states (as transition-metals or
rare-earths compounds). Unfortunately, standard approximations
to Density Functional Theory (like the Local Density Approximation
or the Generalized Gradient Approximation) are not accurate enough
to give a proper description of correlation effects.
While multi-configurational Quantum-Chemistry methods are too
expensive for most systems but very small molecules, several corrections
to DFT have been introduced to give a better description of
correlated states without giving up the relatively low computational cost
of DFT calculations. After a general introduction on electronic correlation,
one of the simplest corrective approaches, named LDA+U, will be described
along with its implementation in a plane-wave pseudopotential
A hands-on tutorial on DFT total energy and force calculations will
cover optimization of input parameters through convergence tests,
electronic relaxation and structural optimization. The simulations
will be performed using the plane-wave pseudopotential implementation of
DFT contained in the Quantum-ESPRESSO package.
Performing constant temperature (and pressure) molecular dynamics simulation requires adequate thermostating (and barostating) algorithms, which couple additional degrees of freedom to the physical phase space variables. Thermostating is particularly important when external work is performed on the system in order to induce a nonequilibrium transition between two states. It appeared recently that nonequilibrium statistical mechanical relations such as the Jarzynski identity and the Crooks theorem naturally emerge as a property of nonautonomous thermostated dynamical systems. This however imposes an additional condition on the form of the thermostating dynamics, which is not required in the equilibrium case and is not fulfilled by all state-of-the-art thermostats.
The potential of mean force (PMF) describes the change in free energy
along a reaction coordinate and determines the strength and likelihood
of association in molecular systems. Current methods, like Thermodynamic
Integration (TI) and WHAM, use piecewise linear and piecewise constant
approximations of the PMF. We propose two methods that allow the use
of arbitrary basis functions, one based on maximum likelihood
estimation and another on the method of weighted residuals. Numerical
comparisons with TI and WHAM are performed.
Joint work with Boris Kozinsky (Department of Physics, MIT),
Nicholas E. Singh-Miller, and Nicola Marzari
(Department of Materials Science and Engineering, MIT).
We address periodic-image errors arising from the use of periodic
boundary conditions to describe systems that do not exhibit full three-
dimensional periodicity. We show that the difference between the
periodic potential, straightforwardly obtained from a Fourier transform,
and the exact potential can be characterized analytically. In light of
this observation, we present an efficient real-space method to correct
periodic-image errors, demonstrating that exponential convergence of the
energy with respect to cell size can be achieved in practical periodic
boundary-condition calculations. Comparing the method with existing
schemes, we find that it is particularly advantageous for studying
charged systems and systems exhibiting partial periodicity.
Free energy computation is one of the main goals of
molecular dynamics simulations. Numerically, these calculations are
made difficult by the presence of free energy barriers separating
meta-stable basins. Several methods have been developed to address
this issue such as the adaptive biasing force (ABF), in which an
external force is progressively adapted. Upon convergence, the
reaction coordinate along which the free energy is computed exhibits a
diffusive behavior in a flat free energy profile, leading to a rapid
convergence of the computed averages and a rapid reduction of
statistical errors. Such methods can serve as the starting point for
computing coarse grained models of large proteins in which the number
of degree of freedom is greatly reduced. We will show the role the
free energy plays in formulating these models. In this case as well,
ABF can be applied to yield faster convergence.
Single molecule magnets are usually based on transition metals with partially filled d shells.
When several metal centers are involved this leads to molecules with many single occupied orbitals
coupled into an intermediate spin state. Conventional methods of quantum chemistry are not able
to deal with this situation, so the Heisenberg model hamiltonian is often used with parameters
estimated from DFT calculations. Practical as well as logical problems with this approach will be
discussed. Some of the difficulties with treating exchange in DFT for open shell systems will be presented.
Objective structures  are a recent classification of atomic and
molecular structures that includes numerous systems of current
interest. Examples are nanotubes and biological virus components.
While these structures are not crystalline, there are many strong
analogies to crystals. These analogies can be exploited to generalize
computational techniques developed in the context of crystals to these
new systems. Further, they may provide strategies for self-assembly
and structure determination.
Some examples of current work exploiting these ideas will be
presented. These include finding the normal modes of these structures
(in analogy with phonons), Objective Molecular Dynamics in analogy
with conventional MD, and DFT and other first-principles techniques in
the Objective setting. Recent work on finding all possible objective
structures will be described.
 R. D. James, Objective Structures, J. Mech. Phys. Solids, 54,
Same abstract as the 7/24 poster session.
Work in collaboration with Eric Cancès (CERMICS-ENPC, France) and Mathieu Lewin (Department of Mathematics, Université de Cergy, France)
We present mathematical results obtained for a new mean-field model
dedicated to the description of interacting electrons in crystals with local
defects. We work with a reduced Hartree-Fock model, obtained from the
usual Hartree-Fock model by neglecting the exchange term. Then we
deduce a new variational model for computing the perturbation in a
basis of precomputed maximally localized Wannier functions of the
reference perfect crystal. Finally we show some numerical results in
which we have applied the variational approximation and used a
one-dimensional model with Yukawa interaction potential.
Same abstract as the 7/24 poster session.
Many systems of significant fundamental and applied interest are irreversible.
These include, but are not limited to, living systems, chemical reactors,
systems with driven flows of matter and energy, and photoactivated systems. For theoretical studies of such non-equilibrium processes, the steady-state
distribution is of central importance because it enables calculation of static
averages of observables for comparison to experimental measurements. For
systems at equilibrium, low probability states can be explored efficiently in
simulations with umbrella sampling methods, in which biasing potentials that
are functions of one or more order parameters are used to enhance sampling of
selected regions of phase space. What complicates extending umbrella sampling
to simulations of non-equilibrium processes is that, by definition, they do not
obey detailed balance (microscopic reversibility). As such, one must account
for the fact that the steady-state probability of observing particular values
of the order parameters can be determined by a balance of flows in phase space
through different possible transitions. In this talk, I will describe the
first algorithm for enforcing equal sampling of different regions of phase
space in an ergodic system arbitrarily far from equilibrium, which enables its
steady-state probability distribution to be determined with high accuracy. The
efficiency of the algorithm will be demonstrated by applying it to a model of a
genetic toggle switch which evolves irreversibly according to a continuous time
Monte Carlo procedure.
A new computational algorithm to extend time
atomically detailed simulations is illustrated. The
(Milestoning) is based on partitioning the dynamics to a
trajectories between “milestones” and constructing a
model for the motion along a reaction coordinate.
Besides "toy" models, two molecular systems are discussed:
kinetics of a conformational transition in a blocked
alanine is computed
and shown to be accurate, more efficient than
Dynamics by a factor of about 9, and non-exponential. (ii)
conformational transition in Scapharca hemoglobin is
results for the rate (about 10+/-9 microseconds) are in
experiment and are obtained about 1,000 times faster than
Molecular Dynamics. No assumption of an activated process
separated by significant barrier is made.
A general scaling argument predicts linear speed-up with
the number of
milestones for diffusive processes, and exponential
transitions over barriers. The algorithm is also trivial to
As a side result Milestoning also produces the free energy
the reaction coordinate, and is able to describe
along one (or a few) degrees of freedom.
Representing the electronic structure in Density Functional Theory (DFT)
by a set of localized wave functions discretized on a real-space mesh
essentially leads to a linear scaling of the computational cost with
the size of the physical system.
This can be achieved by formulating the DFT energy functional
in terms of general non-orthogonal orbitals which are then optimized
under localization constraints (spatial confinement).
Multigrid preconditioning and a block version of Anderson's
extrapolation scheme are used to accelerate convergence towards the
For localization regions --- constraints --- large enough,
one can reduce truncation error to a value smaller than discretization
error and achieve the level of accuracy of a Plane Waves calculation.
Accuracy is improved by allowing for flexible localization regions that
can adapt to the system.
This also reduces problems with local minima and enables energy
conserving Born-Oppenheimer molecular dynamics simulations.
Our implementation of this approach scales on hundreds of processors and
becomes competitive with Plane Waves codes around 500 atoms.
 J.-L. Fattebert and F. Gygi, Phys. Rev. B 73, 115124 (2006)
 J.-L. Fattebert and F. Gygi, Comput. Phys. Comm. 162, 24 (2004)
This work was performed under the auspices of the U.S. Department of
Energy by University of California Lawrence Livermore National
Laboratory under contract No. W-7405-Eng-48.
The accurate representation of the solvent environment is an essential
but also very costly component of biomolecular simulations. As a computationally less expensive alternative, implicit solvent models
based on continuum electrostatic theory have become attractive
alternatives to the standard explicit solvent models. Examples of
biomolecular simulations with implicit aqueous solvent and implicit
biological membranes are presented along with remaining mathematical
challenges to fully take advantage of the much reduced model complexity.
We present a systematic way to construct silicon tips by using
the minima hopping method(MHM) with limited temperature.
In Molecular Dynamics part of the MHM we limit the temperature to
some value used in experiments. We used Tight-Binding
scheme to evaluate energy and forces of silicon and hydrogen atoms
in MHM. Structures with lowest energy obtained with this method
are used in AFM simulation and force between tip and surface is
calculated in terms of distance and some of obtained results
are presented in this poster.
Same abstract as the 7/24 poster session.
A method for applying transient thermal boundary conditions to a
molecular dynamics simulation is presented. The Kapitza effect at the
boundaries is avoided by a diffuse feedback-controlled buffer zone. It
is shown that this technique can be extended to provide a simple
interface for coupling between concurrent atomistic and continuum
simulations for ballistic heat transport.
Wavelets are a systematic localized basis set that is well suited for representing Kohn-Sham orbitals.
I will explain the algorithms that we are using in our new ABINIT wavelet program and show
perfomance results. In the second part I will show how scaling functions can be used to solve
electrostatic problems both for continuous and discrete charge distributions under various boundary conditions.
The multilevel summation method (MSM) computes an approximation to
the pairwise interaction potential and its gradient with an amount of
work that scales linearly as the size of the system. The potential
is smoothly split into a sum of partial potentials of increasing
range and decreasing variability with the longest-range parts
interpolated from grids of increasing coarseness. Multilevel
summation is especially appropriate for dynamical simulations,
because it can produce continuous forces for both nonperiodic and
periodic boundary conditions. A small correction to the grid-based
force approximation can be used to conserve both energy and momentum.
The controlled transport of spin polarised electrons on a 1 nanometre length scale is a realistic prospect and could be the basis for new multifunctional devices with a component density an order of magnitude higher than current VLSI technology. The fundamental materials chemistry challenge is to produce a nano-structured semiconductor that is ferromagnetic at room temperature. Ideally the electronic and magnetic properties need to be robust but tunable through control of composition and structure.
The results of recent high quality theoretical calculations on a number of pure carbon materials will be presented. A novel mechanism for long range magnetic coupling in extended pi-bonded systems will be discussed and documented with explicit calculations on graphene ribbons and defective graphene sheets. A putative ordered defect phase which gives rise to a semiconducting ground state that is ferromagnetic at room temperature will be presented. It will also be shown that the band gap and magnetic coupling may be controlled by varying the defect density.
Joint work with I. Van den Eynde2, J.C. Martins1 and D. Tourwé2
As part of a bigger research effort which consists of the
development of a generic approach for drug development, two
cyclic peptides based on the active sequence of somatostatin
were characterized using nuclear magnetic resonance
spectroscopy (NMR). In order to obtain the 3D conformation of
these molecules, conformational relevant parameters (eg. nOe,
scalar couplings and amide temperature coefficients) were
recorded. A simulated annealing protocol using distance
restraints was used to obtain a set of conformations which were
validated by analysis of the hydrogen bridges present and by
recalculating the measured scalar couplings.
Future efforts will consist of generating mathematical
descriptions (templates) of this kind of small peptide chains
as well as creating templates of non-peptidic small molecules.
This must allow for an intelligent choice for drug templates
upon positive screening results of the small peptide molecules.
1NMR and Structure Analysis Unit, Department of Organic
Chemistry, Ghent University, Krijgslaan 281 S4 B-9000 Gent,
2Unit for Organic Studies, Department of Organic
Chemistry, Vrije Universiteit Brussel, Pleinlaan 2 9G616 B-1050
Same abstract as the 7/24 poster session.
The Gaussian and plane waves (GPW) approach combines
the description of the Kohn-Sham orbitals as a linear combination
of Gaussian functions with a representation of the electron density
in plane waves. The unique properties of Gaussian functions allow
for a fast and accurate calculation of the density in the plane wave
basis. The plane wave representation of the density leads to an
easy solution of Poisson's equation and thereby a representation of
the electrostatic potential. Matrix elements of this potential can
be calculated using the same methods. The auxiliary representation of
the density is further used in the calculation of the exchange-correlation
energy and potential. The resulting approach scales O(N log N) in the
number of electrons and has many additional interesting features,
namely, a small prefactor, early onset of linear scaling, and
a nominal quadratic scaling in the basis set size for fixed system size.
The GPW method is combined with a direct optimization of the subspace
of occupied Kohn-Sham orbitals using an orbital transformation (OT) method.
A variation of this method has recently been implemented that only
requires matrix multiplications. The method combines a small
prefactor with efficient implementation on parallel computers, thereby
shifting the break even point with linear scaling algorithms to
much larger systems. A strategy to combine the OT method with
sparse linear algebra will be outlined.
Joint work with Leonid Gorb,(1
and 2) Mo. Qasim,(2) and
Jerzy Leszczynski. (1)
is one of the most important high energetic nitramines which
are used as explosives and propellants. The decomposition of
CL-20 is very complicated and involves hundreds of elementary
reactions. Accurate knowledge of these reactions and
predictions of their kinetic parameters are important for
modeling these complex processes in combustion and explosion.
However, due to the energetic nature of these materials and the
rapid rates of the intermediate reactions, it is difficult to
monitor these individual reactions experimentally. Recent
advances in first principles modeling have led to enormous
progress toward understanding complex condensed phase chemical
phenomena. Theoretical methods, especially accurate ab initio
molecular dynamics method, provide a viable alternative to
study the dynamics of these reactions. Electronic properties
and the dynamics of the initial thermal decomposition step of
gas-, and e-crystal phases of CL-20 are investigated using the
Car-Parrinello molecular dynamics (CPMD) approach. The
difference in the reaction pathways between gas- and crystal
phases has been accounted for the strong intermolecular
interaction into crystal lattice. The thermal rate constants
have been also estimated.
1) Computational Center of Molecular Structure and Interactions,
Jackson State University, Jackson, MS 39217
2) US Army ERDC, Vicksburg, MS, 39180.
Same abstract as the 7/24 poster session.
A novel Normal-Mode-Partitioned Langevin dynamics integrator is
proposed. The aim is to approximate the kinetics or
thermodynamics of a biomolecule by a reduced model based on a
normal mode decomposition of the dynamical space. The basis set
uses the eigenvectors of a mass re-weighted Hessian matrix
calculated with a biomolecular force field. This particular
choice has the advantage of an ordering according to the
eigenvalues, which has a physical meaning (square-root of the
mode frequency). Low frequency eigenvalues correspond to more
collective motions, whereas the highest frequency eigenvalues
are the limiting factor for the stability of the integrator.
The higher frequency modes are overdamped and relaxed near to
their energy minimum while respecting the subspace of low
frequency dynamical modes. Numerical results confirm that both
sampling and rates are conserved for an implicitly solvated
alanine dipeptide model, with only 30% of the
modes propagated, when compared to the full model. For
implicitly solvated systems the method can be shown to give
improvements in efficiency more than 2 times even for sampling
a small 22 atom (alanine dipeptide) model and in excess of an
order of magnitude for sampling an 882 atom (bovine pancreatic
trypsin inhibitor, or BPTI) model, with good scaling with
system size subject to the number of modes propagated.
The calculation of quantum reaction rates in condensed phase
systems from the perspective of quantum-classical Liouville dynamics
will be described. The rate problem will be cast in a form that
involves quantum equilibrium sampling combined with quantum-classical
dynamics for the evolution. Application to proton transfer in the
condensed phase and polar molecule nanoclusters will be described.
The results of simulations of nonadiabatic reaction rates
using a surface-hopping scheme based on quantum-classical Liouville
dynamics will be presented. The dynamic and structural factors that
influence the rate and reaction mechanism will be discussed.
We propose numerical integration schemes to solve stochastic differential equations driven by two important non-Gaussian noises – dichotomous Markov noise and Poissonian white shot noise. Our formula, which is based on an integral equation, which is equivalent to the stochastic differential equation, utilizes a discrete time approximation with fixed integration time step. We show that our integration formula approaches the Euler formula if these noises approach the Gaussian white noise. We further propose a simplified weak scheme that significantly reduces the computation time, while still satisfying the moment properties up to the required order. Our approach is readily applicable to dynamical systems driven by arbitrary types of noise, provided there exists a way to describe the random increment of the noise accumulated during each integration time step. The accuracy and efficiency of the proposed algorithms are numerically examined.
Same abstract as the 7/24 poster session.
The interfacial free energy is the work required to reversibly form a unit area of interface between two materials and is a principal quantity governing crystal growth kinetics and morphology, nucleation and wetting phenomena. In this talk I will review our recent work on the direct calculation of the interfacial free energy for a number of model systems including crystal/fluid, wall-fluid and wall-crystal systems. For the crystal-melt interface, our method of cleaving walls, based on thermodynamic integration [Phys. Rev. Lett, 85, 4751 (2000)], is of sufficient precision to resolve the often small anisotropy in the interfacial free energy, which is crucial in determining the kinetics and morpholgy of dendritic growth. Our calculations are motivated by the enormous difficulty in obtaining precise interfacial free energies by experiment. We report values of the crystal-melt interfacial free energy for the hard-sphere and Lennard-Jones systems, as well as recent results on the series of inverse-power potentials. In addition, recent work extending the method to wall-fluid and wall-crystal systems is discussed and applied to a recent controversy regarding surface-induced prefreezing in hard-spheres.
To understand biostructures, soft matter, and other abundant sparse
systems, one must account for both strong local atomic bonds and
weak nonlocal van der Waals (vdW) forces between atoms which are sometimes
separated by empty space. A fully nonlocal density functional, vdW-DF [1,3],
now including a self-consistent potential [2,3], will be described.
It has had a number of promising applications , some of which
will be presented, including polymer crystals, metal-organic-framework
structures, and nucleic acids.
 Phys. Rev. Lett. 92, 246401 (2004);  cond-mat/0703442v1;
 Much of the vdW-DF work has been a Chalmers-Rutgers collaboration.
We present a possible approach for the computation of free energies and
ensemble averages in the context of one-dimensional coarse-grained
models in materials science. The approach is based upon a thermodynamic
limit process, and makes use of ergodic type theorems and large
deviation theory. In addition to providing a possible efficient
computational strategy for ensemble averages, the approach in
particular allows for assessing the
validity and the accuracy of some basic assumptions of the
finite-temperature quasicontinuum method.
This is joint work with X. Blanc (Univ. Paris 6), F. Legoll (ENPC Paris), C. Patz (WIAS Berlin).
We introduce a systematic way to construct symplectic schemes for the
numerical integration of a large class of highly oscillatory Hamiltonian
systems. The bottom line of our construction is to consider the Hamilton-Jacobi
form of the Newton equations of motion, and to perform a two-scale
expansion of the solution, for small times and high frequencies.
Several options for the derivation are presented.
The various integrators obtained are tested and compared to several
existing algorithms. The numerical results demonstrate their efficiency.
This is joint work with C. Le Bris (ENPC Paris).
This lecture will introduce the basic types of dynamical models arising in molecular dynamics, and the formulation of structure-preserving numerical methods for integrating the equations of motion. Topics to be covered include: multiple time-scale modelling, the use of constraints, the computation of averages, the role of symplecticness and time-reversal symmetry, the properties of appropriate numerical methods, splitting methods, and SHAKE and RATTLE discretizations. The talk will conclude with a discussion of the use of backward error analysis to correct the numerical bias introduced in dynamical thermostatting methods.
Stochastic dynamics to compute free energy differences are widely used in computational chemistry and biology. Many recent methods rely on complex Markov processes (non-homogeneous or non-linear processes). Examples of such methods are exponential reweighting of non-equilibrium paths (Jarzynski equality) and Adaptive Biasing Force (ABF) techniques. A unifying presentation of adaptive methods is proposed. We also prove the convergence of a certain class of adaptive methods. Finally, we present an efficient implementation of adaptive dynamics using an interacting particle system with birth death processes.
Non-equilibrium molecular dynamics simulations, such as crack propagaion
in solids, requires realistic boundary conditions to guarantee the
accuracy of the computation. In this talk, exact boundary conditions will
be derived using a dimension reduction technique. I will also discuss
practical implementations under this framework.
This presentation describes the "internal” dynamics of
molecular N-body systems in an internal reduced phase space and
describes the coupling of the internal motions with the overall
rotation in the center-of-mass frame. Traditionally, molecular
motions have been separated into translational, rotational, and
internal motions. Symplectic reduction provides a systematic
method for obtaining reduced phase spaces and has been applied
to describe N-body molecular dynamics [1, 2]. This
presentation points out three examples of observations of net
rotation due to internal motions: (a) a net rotation in a
differential geometric study of a triatomic molecule, (b) a net
rotation of 20 degrees in the recoil angle of a departing O
atom in the experimental dissociation of
NO2, and (c) a net
rotation of 42 degrees in 105 reduced
time steps in a
computational study of protein dynamics. For the case of
Jacobi coordinates, the net rotation follows from (a) the
conservation of total rotational angular momentum and,
equivalently, (b) Hamilton’s equations in the center-of-mass
frame . For Eckart generalized coordinates, the net
rotation follows from the conservation of total rotational
angular momentum. A phase space associated with the Eckart
frame is an example of a reduced phase space. Even after
reducing to the “internal” phase space, the reduced internal
dynamics are coupled to the overall rotation in the
center-of-mass frame. In terms of the net rotation, the
strength of the coupling is mass- and coordinate-dependent.
When the total angular momentum vanishes, the coordinates of
overall rotation and internal motions are separated when the
internal angular momentum vanishes, in agreement with a
condition for the separation of energies. Appendix A provides
further applications of the derived net rotation. Appendix B
provides a differential geometric [2, 3] description of the net
 F. J. Lin and J. E. Marsden, Symplectic reduction and
topology for applications in classical molecular dynamics,
Math. Phys., Vol. 33, 1281 – 1294 (1992).
 F. J. Lin, Hamiltonian dynamics of atom-diatomic molecule
complexes and collisions, Discrete and Continuous Dynamical
Systems, Supplement, Vol. 2007, to appear (2007).
 J. E. Marsden, R. Montgomery, and T. Ratiu, Reduction,
symmetry, and phases in mechanics, Mem. Amer. Math. Soc.,
88, No. 436 (American Mathematical Society, Providence,
* This material was presented during the time set aside for Second
Chances on Thursday, July 26, at 4:30 – 4:40 p.m.
The Nose-Hoover thermostat is a deterministic dynamical system
designed for computing phase space integrals for the canonical
Gibbs distribution. Newton's equations are modified by coupling
an additional reservoir variable to the physical variables. The
correct sampling of the phase space according to the Gibbs measure
is dependent on the Nose-Hoover dynamics being ergodic. Hoover
presented numerical experiments that show the Nose-Hoover
dynamics to be non-ergodic when applied to the harmonic
oscillator. We have proven that the Nose-Hoover
thermostat does not give an ergodic dynamics for the
one-dimensional harmonic oscillator when the "mass" of the
reservoir is large. Our proof of non-ergodicity uses KAM theory to
demonstrate the existence of invariant tori for the Nose-Hoover
dynamical system that separate phase space into invariant regions.
We present numerical experiments motivated by our analysis
that seem to show that the dynamics is not ergodic even for a moderate thermostat mass.
We also give numerical experiments of the Nose-Hoover chain (proposed by Martyna, Klein, and Tuckerman) with
two thermostats applied to the one-dimensional harmonic
oscillator. These experiments seem to support the non-ergodicity of the
dynamics if the masses of the reservoirs are large enough and are
consistent with ergodicity for more moderate masses.
Joint work with Frederic Legoll and Richard Moeckel.
Present empirical molecular force fields ignore polarization,
that is the rearrangement of charge distribution due to environment,
and only treat dispersion or vanderWaals forces at the pair level. The
former approximation limits the ability of empirical models to describe
the properties of materials at surface/interfaces. For example, the
shell of a chlorine ion in water can only be predicted by properly
for the polarizability of the water molecules. The latter approximation
is equally problematic at surfaces, yielding large errors in surfaces
which leads to errors in acouning for hydrophobic solvation.
Although ab initio calculations could in principle restore manybody
polarization and manybody dispersion, in fact, the present level of theory
that is sufficiently computational efficient to employ in studies of large
systems, Gradient Corrected Density Functional Theory, does not treat
dispersion even at the pair level. Since the goal of present theoretical
research is to describe the complex interfacial phenomena that drive
processes such as protein folding and transport of materials through
membranes, it is critical to develop novel, linear scale methods, that can
treat both manybody polarization and dispersion. Here, the applied
that underlies the Drude method including a new path integral scheme that
not only models
manybody polarization and dispersion efficient and accurately, but goes
beyond the standard dipole
approximation to include all multipole interactions without requiring an
explicit multipole expansion is presented.
Electronic-structure modeling based on density-functional theory has
become a popular, successful, and reasonably accurate technnique to
characterize and predict material properties. The successful
outcome of a calculations is nevertheless linked to a comprehensive
understanding of the underlying structure of a modern
electronic-structure code, and on the most important parameters that
insure accuracy of the calculation, without sacrificing speed.
We'll overview the basics of the total-energy planewave pseudopotential
method, including issues related to basis set completeness, Brillouin
zone sampling, long-range electrostatic interactions,
exchange-correlation functionals, and minimization techniques.
A hands-on tutorial on first-principles molecular dynamics, using
either Car-Parrinello or Born-Oppenheimer approaches, as
implemented in the CP package of the Quantum-ESPRESSO distribution.
Recent progress in obtaining docked protein complexes will be discussed.
The combination of exhaustive search, clustering and localized global
optimization can reliably find energy minima to highly nonconvex biomolecular
energy functions. Using an energy function that adds desolvation and
screened electrostatics to classical molecular mechanics potentials, the
global minimum is found very near to the observed native state. This is
demonstrated across a large number of benchmark examples.
A number of reasons have resulted in plane-waves becoming one of the
basis sets of choice for simulations based on density-functional theory,
for example: the kinetic energy operator is diagonal in momentum space;
quantities are switched efficiently between real space and momentum space
using fast-Fourier transforms; the atomic forces are calculated by
straightforward application of the Hellmann-Feynman theorem; the
completeness of the basis is controlled systematically with a single
The resulting simulations require a computational effort which scales as
the cube of the system-size, which makes the cost of large-scale
calculations prohibitive. For this reason there has been much interest in
developing methods whose computational cost scales only linearly with
system-size and hence bringing to bear the predictive power of
density-functional calculations on nanoscale systems.
At first sight the extended nature of plane-waves makes them unsuitable
for representing the localised orbitals of linear scaling methods. In
spite of this, we have developed ONETEP (Order-N Electronic Total Energy
Package), a linear-scaling method based on plane-waves which overcomes
the above difficulty and which is able to achieve the same accuracy and
convergence rate as traditional cubic-scaling plane-wave calculations.
A simple algorithm for fast calculation of the Coulombic forces and energies of
point particles with free boundary conditions will be presented. Its
calculation time scales as N log N for N particles. This novel method has lower
crossover point with the full O(N^2) direct summation than the Fast Multipole
Method. The forces obtained by our algorithm are analytical derivatives of the
energy which guarantees energy conservation during a molecular dynamics
simulation. Our method is based on the Poisson solver  for continuous
charge distribution that uses the Deslaurier-Dubuc (interpolating) wavelets.
 L. Genovese, T. Deutsch, A. Neelov, S. Goedecker, and G. Beylkin, Efficient
solution of Poisson's equation with free boundary conditions, J. Chem. Phys.
125, 007415 (2006).
 A. Neelov, S. A. Ghasemi, S. Goedecker, Particle-Particle, Particle-Scaling
function (P3S) algorithm for electrostatic problems in free boundary
conditions, J. Chem. Phys., 2007 (to appear)
Same abstract as the 7/24 poster session.
Joint work with Ben Leimkuhler (School of Mathematics, University of Edinburgh) and Frederic Legoll (LAMI, Ecole Nationale des Ponts et Chaussees, France).
We describe a new dynamical technique for the equilibration of molecular dynamics. Temperature is moderated by a control law and an additional variable, as in Nose dynamics,
but whose influence on the system decreases as the system approaches equilibrium. This device enables approximation of microcanonical averages and autocorrelation functions consistent with given target temperature. Moreover, we demonstrate that the suggested technique is effective for the regulation of heat production in a nonequilibrium setting.
Describing electron correlation effects for large molecules
is a major challenge for quantum chemistry due to the
strong increase of the computational effort with molecular
size. In order to overcome this limitation,
we present a rigorous method based on an AO-formulation of MP2
theory, which allows to avoid the conventional fifth-power scaling
of MO-MP2 theory and to reduce the scaling to linear without
sacrificing accuracy. The key feature of our method are
multipole-based integral estimates (MBIE),
which account for the 1/R coupling in two-electron integrals and
allow to rigorously preselect integral products in AO-MP2 theory.
Here, the magnitude of products decays at least with 1/R**4, so that a
linear-scaling behavior can be achieved by numerical thresholding without
sacrificing any accuracy. The linear-scaling increase of the computational effort is reached
much earlier than for HF or DFT approaches: e.g. the exact behavior of products
indicates a scaling of N**1.0 from one to two DNA base-pairs for a 6-31G* basis.
The number of significant elements in the pseudo-density matrices and of shell pairs
hints to a very similar linear-scaling behavior for larger basis sets studied
up to cc-pVQZ. First results of a preliminary implementation show
that an early crossover to conventional MP2 schemes below two DNA base pairs
is observed, while already for a system with four DNA base pairs wins
are at least a factor of 16.
Joint work with Waldemar Hellmann and Stefan Goedecker (Institute of Physics, University of Basel).
The Bell-Evans-Polanyi principle that is valid for a chemical reaction that
proceeds along the reaction coordinate over the transition state is extended to
molecular dynamics trajectories that in general do not cross the dividing
surface between the initial and the final local minimum at the exact transition
state. Our molecular dynamics Bell-Evans-Polanyi principle states that low
energy molecular dynamics trajectories are more likely to lead into the basin
of attraction of a low energy local minimum than high energy trajectories. In
the context of global optimization schemes based on molecular dynamics our
molecular dynamics Bell-Evans-Polanyi principle implies that using low energy
trajectories one needs to visit a smaller number of distinguishable local
minima before finding the global minimum than when using high energy
Same abstract as the 7/24 poster session.
Density Functional Theory (DFT) is a successful technique used to
determine the electronic structure of matter which is based on a
number of approximations that convert the original n-particle
problem into an effective one-electron system. The end-problem is
essentially a non-linear eigenvalue problem which is solved
iteratively. The challenge comes from the large number of
eigenfunctions to be computed for realistic systems with, say,
hundreds or thousands of electrons. We will discuss techniques for
diagonalization, and for dealing with the nonlinearity. It is
important consider the problem as one of computing an invariant
subspace in the non-linear context of the Kohn-Sham equations. This
viewpoint leads to considerable savings as it de-emphasizes the
accurate computation of individual eigenvectors and focuses instead on
the subspace which they span. We will also discuss other algorithmic
issues encountered in DFT and will offer some thoughts on
Joint work with Volodymyr Babin and Christopher Roland
(Department of Physics, North Carolina State University,
Raleigh, NC 27695-8202).
The Adaptively Biased Molecular Dynamics (ABMD) method
presented here corresponds to the general class of
nonequilibrium dynamics methods that use the history of the
sampling process to bias the dynamics, thereby effectively
flattening the free energy surface of the chosen order
parameter(s). These adaptive methods, are a generalization of
umbrella sampling with a time-dependent potential, and include
the Wang-Landau approach, the adaptive biasing force method,
and non-equilibrium metadynamics (MTD). Recently, we
implemented a variation of MTD in the AMBER package [V. Babin,
C. Roland, T.A. Darden and C. Sagui, J. Chem. Phys.
(2006)]. In spite of the efficient implementation, MTD still
has a square dependence on the sampling time t, i.e., scales as
O(t2) , which can be a significant problem in large systems
with non-negligible entropy. In addition, the large number of
parameters in MTD influences the output in an entangled way.
In the present work, we seek to cure these problems by carrying
out an efficient implementation of ABMD. The method has only
two parameters and scales linearly with time, i.e., it is O(t).
Comparative results for the folding of small peptides are
Jont work with Manoj K. Harbola (Department of Physics, Indian Institute of Technology, Kanpur 208016, India).
The question of whether there exists a mapping from
an excited-state density to a potential is central to performing density-functional calculation for excited states. The motivation of the present work is to establish a unique way of selecting a system (potential) for a given excited-state density. The issue of density – to – potential mapping has been addressed in a series of work by Sahni et al , Harbola  and Gaudoin et al . In the work of  and , it was shown that ground- or excited-state density can be generated by configuration of one’s choice. In the work , they have shown that even with a fixed configuration, one could reproduce same excited-state density from more than one potentials as shown in Fig.1. This implies in addition to the excited-state density one requires some extra information to establish such a unique mapping. This has been done by comparison of ground state densities in Levy-Nagy  theory. Following the earlier attempts we have further explored the density-to-potential mapping based on the work of Levy-Nagy  and Gorling  for excited-states . We have proposed a new criterion , which uniquely establishes such mapping.
Fig.1: Two potentials (middle panel) giving the same excited-state density (upper panel)
along with their corresponding ground-state densities (lower panel) for an excited-state
of the three electron 1D infinitely deep well model system.References:
 V. Sahni, L. Massa, R. Singh and M. Slamet, Phys. Rev. Lett. 87, 113002 (2001).
 M. K. Harbola, Phys. Rev. A 69, 042512 (2004).
 R. Gaudoin and K. Burke, Phys. Rev. Lett. 93, 173001 (2004).
 M. Levy and A. Nagy, Phys. Rev. Lett. 83, 4361 (1999).
 A. Gorling, Phys. Rev. A 59, 3359 (1999).
 P. Samal, M. K. Harbola and A. Holas, Chem. Phys. Lett. 419, 217 (2006)
 P. Samal and M. K. Harbola, J. Phys. B 39, 4065 (2006).
Same abstract as the 7/24 poster session.
Chemists are used to see molecules in three dimensions, and think of the
molecular properties often related to specific regions of space. This is a
good source of inspiration for theoretical methods, but efficient algorithms
for a mathematical treatment of models needing, e.g., integration in
arbitrary, flexible regions in 3D are still needed.
The talk will demonstrate how to extract nanomechanical properties
like stiffnesses and friction parameters of biomolecules from molecular dynamics
simulations. The aim is to construct optimal sets of parameters from molecular
dynamics time series for all important conformations of a biomolecule under
investigation. Parameter optimality will be measured in a maximal likelihood
estimation sense, i.e., such parameters are optimal for which the probability
that the observed time series is an output of equations of motions with these
parameters is the highest possible. It will be demonstrated how to derive the
equation for the optimal parameter set, how these can be implemented
efficiently, and how the resulting algorithm perform for moderately sized
examples (12-alanine with implicit water and some B-DNA 16mer with explicit
Computer simulations of molecular dynamics find many applications in physical chemistry, materials science, biophysics, and structural biology. The basic idea is to simulate the motion of a set of atoms, typically (but not necessarily) in full atomic detail. For this, Newton's equations of motion are used with the forces defined by gradients of a potential energy function that is a rough approximation to that defined by quantum mechanics. A remarkably effective numerical integrator is the simple leapfrog/Stormer/Verlet method. Modeling the boundary of the simulation domain is problematic. Typically, periodicity is assumed, though spherical restraints are sometimes used. The boundary may or may not permit the exchange of energy, momentum, or mass. For the canonical ensemble, only energy is permitted to be exchanged across the boundary. This can be modeled by stochastic boundary conditions. The micro(scopic) state of a system is largely unknown, so initial conditions are chosen at random from a stationary probability distribution of the dynamics. Such a distribution can be obtained for the canonical ensemble by equilibration using Langevin dynamics. Indeed, most quantities of interest are defined only in terms of a stationary distribution and the purpose of the dynamics is to sample configuration space. In some cases, kinetic quantities are of interest and realistic Newtonian dynamics must be used. The practicalities of doing such calculations involves three steps: structure building, simulation, analysis.
Perhaps the greatest computational challenge of molecular dynamics is to generate molecular configurations at random from a prescribed probability distribution, for example, the Boltzmann-Gibbs distribution. Markov chain Monte Carlo methods provide a rigorous solution to this problem, but designing efficient trial moves is challenging. Two systematic approaches are presented: one based on Brownian dynamics, the other on molecular dynamics. The extra computational expense of rejected moves can be avoided by using dynamics instead, the drawback being the introduction of a bias due to the finite stepsize of the integrator. Either deterministic or stochastic dynamics may be used. To achieve reasonable performance, the basic sampling scheme must be combined with more advanced techniques, such as replica exchange and Wang-Landau Monte Carlo.
Density functional theory calculations with a certain class
of approximations to the Kohn-Sham exchange-correlation energy
require an indirect evaluation of the functional derivative
of an implicit functional. Although the formal prescription
for obtaining this derivative is known, there are fundamental
pitfalls in its practical implementation using discrete basis
set representations of the operators. We discuss several
pragmatic solutions to this problem and compare their advantages
in various applications.
(work in collaboration with J.B. Maillet and L. Soulard
(CEA/DAM, Bruyeres le Chatel, France))
I present a model of mesoparticles, in the Dissipative Particle
Dynamics spirit, in which a molecule is replaced by a particle
with an internal thermodynamic degree of freedom (temperature
or energy). The corresponding dynamics, named DPDE and
idependently derived by Avalos and Mackie and by Espanol, is
energy conserving and Galilean invariant. This model is shown
to give quantitatively accurate results for the simulation of
shock waves in a crystalline polymer. I also present an
extension to detonation waves, which are shock waves triggering
exothermic chemical reactions. In this case, an additional
variable per mesoparticle is introduced, namely a progress
variable describing the chemical reaction, and some reversible
second order kinetics depending on the internal energies of the
particles is postulated.
Work in collaboration with Eric Cancès (CERMICS), Ernest R. Davidson (Department of Chemistry, University of Washington), Artur F. Izmaylov, Gustavo Scuseria and Viktor N. Staroverov (Department of Chemistry, Rice University).
This work reviews and presents in a unified fashion several well-known local exchange potentials, such as the Slater potential, Optimized Effective Potentials and their approximations (KLI, CEDA local potentials) and the recently proposed Effective Local Potential. We provide alternative derivations of some of these well-known potentials, mainly based on variational arguments (the local exchange potential being defined as the best approximation of the nonlocal Hartree-Fock operator in some least square sense). The remaining potentials are approximate solutions of the so-called OEP integral equation, and can be recovered through convenient approximations of the resolvent of the Hamiltonian operator.
Hybrid Monte Carlo, introduced in quantum chromodynamics, is a
sampling method that eliminates the bias due to time step
of molecular dynamics (MD). I will discuss the method, give
of its application, and show some recent improvements to its
particularly by discussing the Shadow Hybrid Monte Carlo
Classical modeling depends upon Density Functional Theory (DFT)
for its potential parameters, and quantum modeling usually uses DFT directly. The only significant
improvement in DFT in the last 40 years has been the gradient corrections by Perdew et al., and they
only improved the gross energy errors without improving molecular geometries or forces. Indeed, many
of the talks in the second week address this problem. If we are to get better modeling, we need DFT
The primary problem with DFT functional development is that the total energy errors are quadratic in electron
density errors while the force errors are linear, so that further improvements must get the electron densities
right. The differences between the electron densities predicted by the various theories is at the 1-2% level.
So, in order to get the differences right, one must actually have a numerical method which resolves the
density to 1/10% or the energies to a part per million. Easy-to-use methods which can achieve this accuracy
regularly do not exist, hence my work.
Same abstract as the 7/24 poster session.
The lecture will report on recent progress in combined quantum
mechanical / molecular mechanical (QM/MM) approaches for
chemical reactions in large biomolecules. After a brief outline
the theoretical background and the chosen strategy , we
free-energy QM/MM calculations as well as the use of accurate
correlated ab initio QM methods in QM/MM work. Case studies are
presented for biocatalysis by p-hydroxybenzoate hydroxylase
and cytochrome P450cam [4,5].
 H. M. Senn, W. Thiel, Top. Curr. Chem. 2007, 268,
 H. M. Senn, S. Thiel, W. Thiel, J. Chem. Theory Comput.
 F. Claeyssens, J. N. Harvey, F. R. Manby, R. Mata, A. J.
K. E. Ranaghan, M. Schuetz, S. Thiel, W. Thiel, H.-J.
Angew. Chem. Int. Ed. 2006, 45, 6856-6859.
 J. C. Schoeneboom, F. Neese, W. Thiel, J. Am. Chem. Soc.
 A. Altun, V. Guallar, R. A. Friesner, S. Shaik, W. Thiel,
J. Am. Chem. Soc. 2006, 128, 3924-3925.
In work carried out with Yan Zhao, we have developed a suite of hybrid meta exchange-correlation functionals, including three hybrid meta generalized gradient approximations (hybrid meta GGAs) called M06, M06-2X, and M06-HF and one local meta GGA, called M06-L. The M06 and M06-L functionals are parametrized including both transition metals and nonmetals, whereas the M06-2X and M06-HF functionals are high-nonlocality functionals with double the amount of nonlocal exchange (2X) as compared to M06 and 100% Hartree-Fock exchange, respectively, and they are parametrized only for nonmetals. We have assessed these four functionals by comparing their performance to that of other functionals and other theoretical results for 403 accurate energetic data in 29 diverse databases, including ten databases for thermochemistry, four databases for kinetics, eight databases for noncovalent interactions, three databases for transition metal bonding, one database for metal atom excitation energies, and three databases for molecular excitation energies. We have also tested the performance of these 17 methods for three databases containing 40 bond lengths and for databases containing 38 vibrational frequencies and 15 vibrational zero point energies. We recommend the M06-2X functional for applications involving main-group thermochemistry, kinetics, noncovalent interactions, and electronic excitation energies to valence and Rydberg states. We recommend the M06 functional for applications in organometallic and inorganometallic chemistry and for noncovalent interactions. We recommend the M06-HF functional for all main-group spectroscopy, and we recommend the local M06-L functional for calculations on large systems, where a local functional is very cost efficient. An overview of this work will be presented.
Statistical mechanics provides the connection between the microscopic description
of a system and the macroscopic properties that can be observed in experiments.
In this lecture, the classical statistical mechanics principles that allow us to extract
macroscopic observables from microscopic simulations will be discussed. Statistical
mechanical techniques for analyzing microscopic equations of motion in terms of the
phase-space distributions they generate will be covered. The first part of the lecture
will focus on static equilibrium properties. In the last part of the lecture, time-dependent
statistical mechanics and the calculation of dynamical observables from time correlation
functions will be discussed.
The free energy plays a central role in statistical mechanics and is a key
quantity of experimental interest. The free energy is equal to the reversible
work needed to effect a thermodynamic process and is the
generator of other thermodynamic properties in a given statistical ensemble.
Free energies are directly related to equilibrium constants (for example,
the efficacy of a drug is measured by an equilibrium constant known as the
inhibition constant, which is related to the binding free energy between an inhibitor
and its target enzyme) and chemical potentials, and can be used to estimate
rate constants. Obtaining accurate free energies is challenging because of
the inability of straightforward simulation techniques to sample phase space
sufficiently. In this lecture, we will introduce the concept of free energy and discuss
a number of state-of-the-art techniques for enhancing phase-space sampling,
particularly when the sampling is hindered by so-called rare events.
Same abstract as the 7/24 poster session.
Joint work with Guoai Pan (National Institute for Nanotechnology, Canada),
Thomas Ihle and Daniel M. Kroll (Department of Physics,
North Dakota State University).
A recently introduced particle-based model for fluid dynamics
with continuous velocities is generalized to model
immiscible binary mixtures. Excluded volume interactions
between the two components are modeled by stochastic
multiparticle collisions which depend on the local velocities
and densities. Momentum and energy are conserved
locally, and entropically driven phase separation occurs for
high collision rates. An explicit expression for the equation
of state is derived, and the concentration dependence of the
bulk free energy is shown to be the same as that of the
Widom-Rowlinson model. Analytic results for the phase diagram
are in excellent agreement with simulation data.
Results for the line tension obtained from the analysis of the
capillary wave spectrum of a droplet agree with measurements
based on the Laplace's equation. The dispersion relation for
the capillary waves is derived and compared with the
numerical measurements of the time correlations of the radial
fluctuations in the damped and over-damped limits.
The introduction of "amphiphilic" dimers makes it possible to
model the phase behavior of ternary surfactant mixtures.
Displacement cascades are the primary source of radiation damage of reactor structural materials exposed to fast-neutron irradiation. They are formed by the recoil of primary knock-on atoms with a kinetic energy of more than ~1 keV. The cascade process is characterized by lengths and times of the order of nm and ps, respectively, and it can be modelled by the method of molecular dynamics.
At a time in the early stage of a cascade, only a small proportion of atoms recoil with high velocity while the rest satisfy the equilibrium velocity distribution for the ambient temperature. Because the convergent integration timestep is determined by the velocity of the fastest atom, conventional techniques that use fixed step-size integration over the entire ensemble are impractical at that stage. To obtain reasonable computational efficiency we decomposed the ensemble into evolving subsets of ‘cold’ and ‘hot’ atoms and conducted integration of the equations of motion using variable time-step. Four different techniques were applied to identify point defects and point defect clusters created in displacement cascades and determine the size of clusters and their morphology.
A significant problem in the atomistic simulation of materials is that
molecular dynamics simulations are limited to nanoseconds, while important
reactions and diffusive events often occur on time scales of microseconds
and longer. Although rate constants for slow events can be computed directly
using transition state theory (with dynamical corrections, if desired,
to give exact rates), this requires first knowing the transition state.
Often, however, we cannot even guess what events will occur. For example,
in vapor-deposited metallic surface growth, surprisingly complicated exchange
events are pervasive. I will discuss our accelerated molecular dynamics approach
for treating this problem of complex, infrequent-event processes. The idea is to directly
accelerate the dynamics to achieve longer times without prior knowledge
of the available reaction paths. With these methods, we can often achieve accurate dynamical
evolution on time scales exceeding that available to direct molecular dynamics by many orders
of magnitude. Time permitting, I will also discuss some of our most recent method
Ab initio Molecular Dynamics has been recognized lately as a powerful tool to compute infrared spectra of a variety of systems. The main advantage of this approach with respect to the usual normal mode analysis approach at the optimized geometry is the explicit treatment of temperature and environmental effects, without need for an harmonic approxiamation. However, the interpretation of the simulated spectra in temrs of atomic motions, phonons or vibrational modes is still a challenge.
Here, a general method for obtaining effective normal modes from Molecular Dynamics simulations is presented. The method is based on a localization criterion for the fourier transformed velocity time-correlation functions (FTVCF) of the effective modes. For a given choice of the localization function used, the method becomes equivalent to the principal mode analysis (PMA) based on covariance matrix diagonalization. On the other hand, proper choice of the localization function leads to a novel method with strong analogy with the usual normal mode analysis of equilibrium structures, where the system's Hessian at the minimum energy structure is replaced by the thermal averaged Hessian, although the Hessian itself is never actually calculated. This method does not introduce any extra numerical cost during the simulation and bears the same simplicity as PMA itself. It can thus be readily applied to ab initio Molecular Dynamics simulations. Some examples will be given here.
Density functional theory (DFT) has been firmly established as one of the most widely used first-principles quantum mechanical methods in many fields. Each of the two ways of solving the DFT problem, i.e., the traditional orbital-based Kohn-Sham (KS) and the orbital-free (OF)  schemes, has its own strengths and weaknesses. We have developed a new implementation of DFT, namely orbital-corrected OF-DFT (OO-DFT) , which coalesces the advantages and avoids the drawbacks of OF-DFT and KS-DFT and allows systems within different chemical bonding environment to be studied at a much lower cost than the traditional self-consistent KS-DFT method. For the cubic-diamond Si and the face-centered-cubic Ag systems, OO-DFT accomplishes the accuracy comparable to fully self-consistent KS-DFT with at most two non-self-consistent iterations  via accurately evaluating the total electronic energy before reaching the full self-consistency [2-5]. Furthermore, OO-DFT can achieve linear scaling by employing currently available linear-scaling KS-DFT algorithms and may provide a powerful tool to treat large systems of thousands of atoms within different chemical bonding environment much more efficiently than other currently available linear-scaling DFT methods. Our work also provides a new impetus to further improve OF-DFT method currently available in the literature.
 Y. A. Wang and E. A. Carter, in Theoretical Methods in Condensed Phase Chemistry, edited by S. D. Schwartz (Kluwer, Dordrecht, 2000), p. 117.
 B. Zhou and Y. A. Wang, J. Chem. Phys. 124, 081107 (2006). (Communication)
 “An Accurate Total Energy Density Functional,” B. Zhou and Y. A. Wang, Int. J. Quantum Chem. (in press).
 “The Total Energy Evaluation in the Strutinsky Shell Correction Method,” B. Zhou and Y. A. Wang, J. Chem. Phys. (in press).
 “Accelerating the Convergence of the Total Energy Evaluation in Density Functional Theory Calculations,” B. Zhou and Y. A. Wang, J. Chem. Phys. (submitted).
DFT based approaches permit the determination of structural and thermodynamic properties of materials with sufficiently useful accuracy to allow one to address states and properties of planetary interiors. I will make a brief review of areas in mineral physics problems that have recently experienced much progress and are shedding light on fundamental problems in planetary sciences.
Research supported by NSF/EAR and NSF/ITR programs.
Practical applications of one-electron equations for embedded orbitals (Eqs. 20-21 in Ref. )
hinge on the availability of explicit density functionals to approximate adequately the exchange-correlation energy and the non-additive kinetic energy. The former quantity is defined as in the Kohn-Sham formulation
of density functional theory, whereas the latter one arises from the use of orbitals (/embedded orbitals/) for only
a selected component of the total electron density in the applied formal framework. The quality of the /shifts /of the electronic properties of a chemical species due to its condensed phase
environment calculated by means of Eqs. 20-21 of Ref.  is determined by the kinetic-energy-functional
dependent component of the total effective potential.
In this work, our recent works concerning the development and testing of system-independent
approximations this component of the embedding potential. and selected representative applications
to study details of the electronic structure of embedded systems in condensed phase [2,3] are reviewed.
 T.A. Wesolowski & A. Warshel, /J. Phys. Chem./ *97* (1993), 8050.
 M. Zbiri, C. Daul, and T.A. Wesolowski, /Journal of Chemical Theory and Computation / *2* (2006) 1106.
 J. Neugebauer, C.R. Jacob, T.A. Wesolowski, E.J. Baerends, /J. Phys. Chem. A./ *109* (2005) 7805.
Same abstract as the 7/24 poster session.
This poster reports some progresses we made recently in developing numerical algorithms and program packages for minimizing large scale biomolecular potential energy functions, which is one of the fundamental tasks in biomolecular simulations.
The new algorithms are designed from a sparse incomplete Hessian matrix constructed by a spherical cutoff strategy according to the problem size and computer capability. They are called the incomplete Hessian Newton method (IHN), the truncated-IHN method (T-IHN), and the refined T-IHN method (RT-IHN), respectively. The program package of RT-IHN for solving the biomolecular potential energy minimization problem is developed based on a widely-used biomolecular simulation packages, CHARMM. This poster will give RT-IHN a detailed description, and present its numerical performances for some large scale protein systems.
In theory, these new algorithms are defined and analyzed for a general unconstrained minimization problem with a target function having a dense Hessian matrix. They are proved to be convergent globally. While IHN and T-IHN are shown to have a linear rate of convergence, RT-IHN is proved to have a super-linear rate of convergence. Hence, these new algorithms can be applied to other large scale scientific and engineering applications. For this purpose, we develop a general program package of RT-IHN in C++ based on the optimization program package TAO (http://www-unix.mcs.anl.gov/tao/index.html), and apply it to solve the chemical database optimal projection mapping problem.
This is the joint work with Mazen G. Zarrouk and Jun Wang.
A short overview over solid-solid phase transitions ("martensitic
transformations") will be given, and a one-dimensional lattice model
will be analysed. The equations of motions are given by Newton's
law, where the force is given by the interaction of neighbouring atoms.
The challenge is that the interaction potential is nonconvex, to model
phase transitions. The existence of subsonic travelling waves can be
proved. We will explain how this information can be used to derive
so-called kinetic relations, which relate the velocity of the phase
boundary to a driving force. The interest in kinetic relations can
be explained by their relevance for a continuum modelling of phase
transitions; this will be briefly sketched.
This is joint work with Hartmut Schwetlick (Bath).