Transcript CPMD : Car-Parrinello molecular dynamics

(a short) CPMD tutorial
F.L. Gervasio
ETH Zurich c/o USI campus,
Via Buffi 13, 6900 Lugano, CH
Acknowledgments:
This tutorial is partially based on that initially put together
with Carme Rovira, Roger Rousseau, Axel Kohlmeyer and
Alessandro Laio presented at Centre for Research in
Theoretical Chemistry of the Parc Cientific de Barcelona in
2004.
Prof. Nicola Marzari (MIT) for useful discussions
Why Ab-initio Molecular Dynamics?
Classical molecular dynamics using predefined potentials is well
established as a powerful tool to investigate many-body condensed matter
systems. Despite overwhelming success a fixed model potential implies serious
drawbacks:
1 Many different atom or molecule types give rise to a myriad of different
inter-atomic interactions that have to be parameterized.
2 The electronic structure/bonding pattern changes qualitatively in the
course of the simulation.
3 Systems at very high temperatures/pressures for which no experimental
data is available.
Why Ab-initio Molecular Dynamics?
Publication and citation analysis.
Squares: number of publications which
appeared up to the year n that contain
the keyword “ab initio molecular
dynamics" or synonyma in title,
abstract or keyword list.
Circles: number of publications which
appeared up to the year n that cite the
1985 paper by Car and Parrinello
[R. Car and M. Parrinello, Phys. Rev.
Lett. 55, 2471 (1985).]
What is CPMD?
The CPMD code is a plane wave/pseudopotential implementation of
DFT for ab-initio molecular dynamics.
First version by Jurg Hutter at IBM Zurich Research Lab.
Many people contributed to the development: Michele Parrinello, Jurg
Hutter, D. Marx, P. Focher, M. Tuckerman, W. Andreoni, A. Curioni, E.
Fois, U. Roetlisberger, P. Giannozzi, T. Deutsch, A. Alavi,
D.Sebastiani, A. Laio, J. VandeVondele, A. Seitsonen, S. Billeter and
others.
What is CPMD?
The current version, 3.11.1, is copyrighted jointly by IBM Corp and by
Max Planck Institute, Stuttgart.
It is distributed free of charge to non-profit organizations
(http://www.cpmd.org/)
CPMD runs on many different computer architectures and it is well
parallelized (MPI and Mixed MPI/SMP).
CPMD characteristics
LDA, LSD and many gradient correction schemes (BLYP,HTCH,PBE,etc)
Norm conserving or ultrasoft pseudopotentials
Free energy density functional implementation
Isolated systems, periodic boundary conditions, k-points
Molecular and crystal symmetry
CPMD capabilities
Wavefunction optimization: direct minimization and diagonalization
Geometry optimization: local optimization and simulated annealing
Molecular dynamics: NVE, NVT, NPT
Path integral MD
Response functions
Excited states
Time-dependent DFT (excitations, MD in excited states)
Coarse-grained non-Markovian metadynamics
Wannier, EPR, Vibrational analysis
QM/MM
See on-line manual at: http://www.cpmd.org/cpmd_on_line_manual.html
Theory: Infos and Literature
Car-Parrinello molecular dynamics (CP-MD) simulations bring together
methods from classical molecular dynamics (MD), solid state physics and
quantum chemistry, so some knowledge in all of these areas is needed.
Reviews:
D.K. Remler and P.A. Madden, Mol. Phys. 70, 921ff. (1990)
M.C. Payne, M.P. Teter, D.C. Allen, T.A. Arias and J.D. Joannopoulos, Rev. Mod.
Phys. 64, 1045-1097 (1992)
D. Marx and J. Hutter, Forschungszentrum Jülich, NIC Series, Vol. 1 (2000), 301449 http://www.fz-juelich.de/nic-series/Volume3/marx.pdf
J. Kohanoff and N. Gidopoulos, Handbook of Molecular Physics and Quantum
Chemistry, ed. Stephen Wilson. Volume 2, Part 5, Chapter 26, pp 532-568 (Wiley,
Chichester, 2003)
J. Kohanoff: “Electronic Structure Calculations for Solids and Molecules” 2006
Cambridge university Press
Webpages:
http://www.pci.unizh.ch/gruppe.hutter/e/information.html
http://www.cpmd.org/cpmd_thecode.html
http://www.theochem.ruhr-uni-bochum.de/~axel.kohlmeyer/cpmd-tutor/
Molecular Dynamics
•Initial geometry, velocities
•Potential U(R)
•Evolve the trajectory on the
Phase space using
Newton’s second law.
The potential depends only on the positions
The total energy is conserved
Initial structure?
The Initial geometry can be obtained from:
•Experiment: PDB database, Cambridge database, etc.
•Drawing it from scratch
It is better that you pre-optimize it by using a MM force field,
If the initial geometry is completely wrong the wavefunction
Optimization will take an extremely long time or even fail
(see WF optimization).
Potentials from electronic structure theory
• Born-Oppenheimer approximation
• Solution to the electronic time-independent Schroedinger eq.
• For a given electronic state we get
Car-Parrinello MD
L
BOMD
I
1
M I RI2   o Hˆ  o
2
DFT

 
 
min
 E n(r ), RI
n(r )
CPMD
LCP

 
 2
 
1
1
  M I RI   i  i  i  E n(r ), RI  constr.
I 2
I 2
Ficticious electronic mass
  
ij
Electronic velocities
kinetic energy of the electrons
i
 j   ij

i, j
orthonormalization
Car-Parrinello MD
Necessary condition for CP-MD
 i (r )
E
 
RI
The electrons have to adiabatically
follow the nuclei
If there is an energy transfer from the nuclei to the elctrons
The electronic temperature increases
The system leaves the BO surface
RICP (t )  RIBO (t )
Car-Parrinello MD
Energy conservation
• BOMD
BO
TOT
E

I
 
1
M I RI2  E  i  , RI 


2
• CPMD
CP
BO
ETOT
 ETOT

1
BO




E

i
i
i
TOT  Tel
2 i
CP
BO
BO
ETOT
 cte.  ETOT
 Tel  ETOT
BO
If Tel  ETOT
Car-Parrinello MD
Why the CP method works?
Nuclear and electronic
subsystems decoupled
No superposition of their
power spectra
E.g.:diatomic molecule
decrease 
increase mnuclear
wnuclear
To control the adiabatic
decoupling
min
increase we
welectronic
wemin 
w
homolumo
E gap

Better adiabaticity
decreasing 
decrease Dtmax (Dtmax 
1
w
max
e


Ecut
Best compromise
)
Car-Parrinello MD
Comparison of CP and BOMD forces (Si crystal )
(Pastore, Smargiassi and Buda, Phys. Rev. A 44, 6334, 1991)
F CP  F BO  0
CP: continous line
BO: points
The difference is small and oscillatory
F CP  F BO
Car-Parrinello MD
In general
Egap ~ 1 eV
 = 400 - 800 a.u.
w
Dt = 5 - 10 a.u. = 0.1 - 0.2 fs
max
nucl .
 2000  3500 cm1
Attention to metals (Egap= 0)
NO adiabatic separation
BOMD or CPMD free energy DFT
Check the following properties
• Tnuclei >> Telectrons
• Telectrons = constant
• Etotal = constant
CPMD vs BOMD
CPMD
Timestep,
Dt ( fs)
BOMD
0.1-0.2

0.5-2

Car-Parrinello MD
Latest results on the effect of 
Car-Parrinello MD
Latest results on the effect of 
(i)[ …] great care is required in the choice of an appropriate fictitious electron mass. Values of  of 800 au
(for H2O) or 1100 au (for D2O) are inappropriate for simulations longer than a few picoseconds. A value of
400 au should be used.
(ii) Simulations in the microcanonical ensemble for small systems near a phase transition can lead to
problematic (nonergodic) behavior that gives reproducibility problems. The same applies to empirical water
models.
(iii) BOMD simulations suffer from the same problems regarding reproducibility as CPMD simulations.
An advantage of the CPMD approach using a Lagrangian including the fictitious electronic degrees of
freedom is that it allows us to conserve the total energy extremely well, better than BOMD.
(iv) Overall, simulations in the canonical ensemble offer distinct advantages: instabilities and sensitivity to
initial conditions of microcanonical simulations are avoided and the use of larger fictitious electron masses
is permissible in CPMD simulations (massive thermostats).
The computational benefits of CPMD sampling over Born-Oppenheimer sampling for structural and
thermodynamic properties are appreciable for liquid water and arise from the substantial reduction in
computer time per step (by about a factor of 20) that is only partially compensated by the requirement for a
smaller time step (by about a factor of 5).
The simulations presented here show that classical trajectories for 64-molecule systems using the BLYPTM and BLYP-GTH descriptions of water (at T = 315 K and F ) 1.0 g/cm3) yield an overstructured liquid
(gOO max too high by about 0.2 units), an underestimated heat capacity, and an underestimated self
diffusion constant. However better than PBE functional (gOO max too high by about 0.8 units).
Using massive thermostats on the ions and electrons gives
reliable results even with relatively large 
The error on the Functional is much larger than any effect of 
Density Functional Theory with PW
BO-Approximation:
Hohenberg and Kohn: proof that the density uniquely determines
the energy of the system. Physical Review 136 B864 (1964).
Kohn and Sham: that the variational search for the density which
provides the lowest energy may be preformed using single particle
wavefunctions. Physical Review 140, A1133 (1965).
Density is sum of square of orbitals.
Electronic energy in terms of electronic density
One can treat this as an eigenvalue problem:
Plane Waves
The KS orbitals in a periodic system may be expanded in a basis set of G vectors
With this basis set the KS equations take on the form:
Where the first term is the kinetic energy operator, Te, and it is diagonal in a basis
of G vectors and we only need to solve an NelecXNplanewaves problem.
Note Te is defined also for non-reciprical lattice vectors k. These vectors
form the Brillouin Zone or k-space and are usually sampled on a grid
of k-points.
Evaluating Kohn Sham Equations:
Planewave expansions
For a periodic 3D potential the kinetic energy operator has eigenfucntions
of planewaves with reciprocal vector G and eigenvalues |G|2.
Since EαG2 then we can define the number of basis functions by the cutoff
energy Ecut. For historical reasons specified in rydberg (Ry)=0.5 a. u.
Estimating the number of plane waves for a cell with a G-space volume Ω
we get:
To construct the electronic density in terms of G vectors we have to use
A doubly large set.
Note: density expansion depends only on the box size and linear in Nelec
-not the number of atoms
When will this be a good basis set:
The pseudopotential.
Plane waves are a good basis set if the functions we want to describe are
smooth-i.e. the eigenfunctions of a smoother potential.
However, atomic wavefunctions are not so nice for core orbitals and
to radial nodes of valence orbitals in the core region.
Solution replace core orbitals by pseudopotential and remove radial nodes.
Steep potential=bad
flat potential=good
The choice of pseudopotentials
CPMD works with several types of pseudopotentials (Norm conserving, Ultra soft)
The most commonly used for non-metallic atoms are
Martin Troullier [Phys. Rev. B, 43, 1993 (1991)].
MT are the preferred norm-conserving potentials.
Work well for light main group elements: C, N, O, S, P etc
They have EXC built in: if you change functional you must change pseudo.
They are subject to “ghost states” from KB seperation.
They are not really flexible enough to obtain transferable potentials for transition
metals.
Typical cut-off: 70Ry
With transition metals you need to use Goedecker-Hutter or Vanderbilt ultra soft.
GH require high cutoff, cover most elements
Vanderbilt: low cut-off~40Ry, more calculations required, in CPMD many features not
implemented.
On www.cpmd.org under /contributed you will find a pseudopotential library.
When you use a new pseudopotential always check it against known properties:
geometry, vibrational properties, lattice constants, bulk modulus (better all-electron).
Evaluating Kohn Sham Equations:
Overview
Input: Structure, box, Ecut
G-vectors and Real space Grid
Initail Guess
Construct: VKS
SCF
Diagonalize KS equations
Calculate EKS
Calculate SCF criterium
OK,
Write Restart
EXIT
Comparing to Gaussians
The Good:
Has a better scaling than Gaussian based methods-without loss of accuracy.
Planewaves are easy to program-fast math libraries. No Pulay Forces.
Can do solids, polymers and molecules all in same framework.
No BSSE (Mad Gaussian Disease).
Other properties, not discussed here, are easy to calculate with PW.
The Bad:
Care must be taken with real space grid and Ecut
Ripple noise makes it hard to converge structures.
To do isolated molecules requires a lot of work: screening of images.
Need lots of planewaves.
Chemists have to learn Solid state physics lingo-culture gap.
The Ugly:
Exact Exchange is non-trivial and very expensive to include in this formalism.
Pseudopotentials!
CPMD basics
 To run this examples you need the compiled executable, the input
and the pseudopotentials
 Compiling CPMD is simple on most architectures
1. Make sure that you have the necessary libraries (MPI, BLAS)
2. Untar your package
3. Run ./mkconfig.sh (PC-IFC-P4) > Makefile
Some hints to compile a scalar version on linux can be found on
http://www.theochem.ruhr−uni−bochum.de/~axel.kohlmeyer/cpmd−linux
.html.
 To run (a scalar version of) cpmd you just call the executable
followed by the input filename
Wavefunction optimization
•Every CPMD run starts by optimizing
the wavefunction or from a RESTART.
•In this example we will optimize the
hydrogen electronic structure
•The input is organized in sections that
begin with &NAME and end with &END
•All commands MUST be in UPPER
case otherwise ignored
•A minimal input should have the
sections &CPMD, &SYSTEM &ATOMS
(see the online manual)
Wavefunction optimization
Wavefunction optimization
Starting from the initial guess based on atomic wavefunctions the
wavefunction for the total system is now calculated with an optimization
procedure. You can follow the progress of the optimization in the output file.
NFI: Step number
GEMAX: largest off−diagonal component
CNORM:average of the off−diagonal comp.
ETOT: total energy
DETOT: change in total energy from previous step
TCPU: (CPU) time for this step.
Wavefunction optimization
The calculation stops after the convergence criterion of 1.0d−7 has been
reached for the GEMAX value.
Although the calculation started
with the experimental H−H bond
length there are still some forces
in the direction of the molecular
axis.
Note, that regardless of the input
units, coordinates in the CPMD
output are always in atomic
units.
Other important OUTPUT files:
RESTART.1, LATEST contain all the information on the final state of the system
GEOMETRY.xyz contains the coordinates of the atoms
Wavefunction optimization
There are several ways to optimize the wavefunction in CPMD
The default uses Direct Inversion of the Iterative Subspace (DIIS)
[P. Pulay, Chem. Phys. Lett. 73, 393 (1980). ] which is usually the faster.
Any optimization that takes more than 100 steps should be considered slow.
If the ODIIS converger gets stuck (more than one reset) stop the run and
restart using
PCG
MINIMIZE
TIMESTEP 20
the conjugate gradient minimizer with line search is much more robust.
Starting a CPMD from a random wavefunction with all atom positions fixed, a
comparatively high electron mass and using ANNEALING ELECTRONS is
another alternative to get to a reasonably converged wavefunction.
Wavefunction optimizations for geometries that are far from equilibrium are
often difficult. You can relax the convergence criteria to or and do some
geometry steps. After that optimization will be easier.
Geometry optimization
A geometry optimization is not much else than repeated single point
calculations, where the positions of the
atoms are updated according to the forces acting on them.
We replaced WAVEFUNCTION with GEOMETRY and added the suboption XYZ to have CPMD write a 'trajectory' of the optimization in a
file name GEO_OPT.xyz .
We also specify the convergence parameter for the geometry.
Please notice:
In the case of large systems this way of optimizing the geometry is very inefficient.
It is much better to perform a Car-Parrinello dynamics with ANNEALING IONS.
Geometry optimization
This run will take a little longer,
than the previous.
Here we have to do multiple
wavefunction optimizations.
In the output you can see, that
after printing the positions and
forces of the atoms there is a
small report block.
The numbers for GNMAX,
GNORM, and CNSTR stand for
the largest absolute component
of the force on any atom,
average force on the atoms,
and the largest component of a
constraint force on the atoms
respectively.
Car-Parrinello MD
Restarting from previous coordinates and optimized wavefunctions, we can now
perform a CPMD
The keyword MOLECULAR DYNAMICS CP defines the job type.
We tell the program to pick up the previously calculated wavefunction and
coordinates from the latest restart file (which is named RESTART.1 by default)
The temperature of the system will be initialized to 50K. The time step is set to 4
atomic units (~0.1 femtoseconds). MAXSTEP limits the MD to 200 steps.
Car-Parrinello MD
In the CPMD code atoms are sometimes referred to as ions. This is due to the
pseudopotential approach, where you integrate the core electrons into the
(pseudo)atom which then could be also described as an ion.
NFI: Step number (number of finite iterations) EKINC: (fictitious) kinetic energy of the
electronic (sub−)system TEMPP: Temperature (= kinetic energy / degrees of freedom) for
atoms (ions) EKS: Kohn−Sham Energy, equivalent to the potential energy in classical MD
ECLASSIC: Equivalent to the total energy in a classical MD (ECLASSIC = EHAM −
EKINC) EHAM: total energy, should be conserved DIS: mean squared displacement of
the atoms from the initial coordinates. TCPU: (CPU) time needed for this step.
Car-Parrinello MD
The plot shows the evolution of the various energies. Little energy from the ionic system
is transferred to the fictitious electron dynamics (the temperature is always less than the
initial). The difference between the orange (EHAM) and the blue (ECLASSIC) graphs is
EKINC, and the difference to the potential energy (EKS) is kinetic energy in the ionic
system.
Other Job Types
There are several further types of calculations possible with CPMD, those above are an
example. Please check out the CPMD manual, the CPMD mailing list archives for more
information on how to perform them.
A great source of useful examples is the CPMD test suite (on www.cpmd.org).
Water Molecule
Now we will study a more ambitious molecule: water.
Since water has a dipole moment, you have to keep in mind, that we are
calculating a system with periodic boundary conditions, so the water molecule
'sees' its images and interacts with them.*
&CPMD
OPTIMIZE GEOMETRY XYZ
HESSIAN UNITY
CONVERGENCE ORBITALS 1.0d-7
CONVERGENCE GEOMETRY 3.0d-4
ODIIS 5
MAXSTEP 100
MAXCPUTIME 1500
STRUCTURE BONDS ANGLES
&END
*There are methods implemented in CPMD to compensate for this effect In &SYSTEM
use SYMMETRY 0 and set POISSON SOLVER HOCKNEY or TUCKERMAN
Water Molecule
This time we will use a gradient
corrected functional (BLYP) instead of
the LDA.
Also note that in the &ATOMS section
the LMAX for the oxygen is set to P
(instead of S for hydrogen) and that the
keyword KLEINMAN−BYLANDER is
required for for the calculation of the
nonlocal parts of the pseudopotential.
&DFT
FUNCTIONAL BLYP
GC-CUTOFF 1.0d-06
&END
&SYSTEM
SYMMETRY 1
CELL 20.0 1.0 1.0 0.0 0.0 0.0
CUTOFF 70.0
&END
&ATOMS
*O_MT_BLYP.psp KLEINMAN-BYLANDER
LMAX=P
1
10.0 10.0 10.0
*H_CVB_BLYP.psp
LMAX=S
2
8.5 9.0 10.0
11.5 9.0 10.0
&END
Water Molecule
We can now do a
properties calculation
using the RESTART
from the previous run
&CPMD
PROPERTIES
RESTART WAVEFUNCTION COORDINATES
WANNIER WFNOUT ALL
WANNIER REFERENCE 10.0 10.0 10.0
&END
&PROP
LOCALIZE
PROJECT WAVEFUNCTION
DIPOLE MOMENT
CHARGES
&END
&DFT
FUNCTIONAL BLYP
GC-CUTOFF 1.0d-06
&END
&SYSTEM …
&ATOMS …
Concluding Remarks
 This has been only a short Introduction to CPMD
and PW-DFT.
 Mixed basis Gaussian/PW codes can combine
advantages of both.
 These methods can be made linear in scalability.
 Currently the community is shifting toward
these approaches (cp2k http://cp2k.berlios.de/index.html).