Basic usage of QM/MM in CPMD

Carme Rovira
Universitat de Barcelona (UB) &
Institució Catalana de Recerca i Estudis Avançats (ICREA)

1
 The CPMD code
([Link])
• Based on the original code by Car and Parrinello [1].

It is a production code with many unique features written in Fortran 77/90 with
about 150,000 lines of code (in 2005). The code is parallelized using MPI and
OpenMP.

• Since January of 2002 the program is freely available for non-commercial use [2].

[1] Car, R.; Parrinello, M.: Unified Approach for Molecular Dynamics and Density-Functional Theory.
Phys. Rev. Lett. 55 (1985) 2471–2474.
[2] CPMD V4.1 Copyright IBM Corp 1990–2004, Copyright MPI für Festkörperforschung Stuttgart
1997–2001, see also [Link]

See also:
J. Hutter and M. Iannuzzi. CPMD: Car-Parrinello molecular dynamics.
Z. Kristallogr. 220 (2005) 549–551
J. Hutter. Car-Parrinello molecular dynamics. WIREs Comput Mol Sci 2012, 2:604-612.
D. Marx, J. Hutter, Ab Initio Molecular Dynamics. Basic Theory and Advanced Methods.
Cambridge University Press. 2012
 Structure of the CPMD input file
SECTIONS
&CPMD

&END

&SYSTEM

&END

&DFT

&END

&ATOMS

&END
 Structure of the CPMD input file
“minimum” &CPMD contents
&CPMD
OPTIMIZE WAVEFUNCTION E = min E [ρ(r), {R I }]
&END ρ(r)

occ. 2
ρ(r) = ∑ ψi (r)
dumps out RESTART.1 and LATEST files i

∑ ci (G) exp[iG.r]
ψi (r )=
1
ΩG
./LATEST

Contains the path for

the restart file
 Structure of the CPMD input file
“minimum” &CPMD contents
&CPMD
OPTIMIZE WAVEFUNCTION E = min E [ρ(r), {R I }]
&END ρ(r)

&CPMD
OPTIMIZE GEOMETRY
RESTART WAVEFUNCTION LATEST
&END RESTART.1 and
LATEST files are
altogether always dumped out

&CPMD
OPTIMIZE WAVEFUNCTION GEOMETRY
&END
Structure of the CPMD input file:
The CPMD section
&CPMD
MOLECULAR DYNAMICS
RESTART WAVEFUNCTION COORDINATES LATEST
&END

&CPMD
MOLECULAR DYNAMICS
RESTART WAVEFUNCTION COORDINATES VELOCITIES LATEST
quench electrons bo In small letters => not used
RESTART ACCUMULATORS LATEST
&END
You can also restart other quantities (not covered in this
tutorial)
RESTART NOSEP NOSEE CELL LATEST
 Structure of the CPMD input file:
The CPMD section
(similar input as the one we’ll work
&CPMD during the afternoon practicals)
MOLECULAR DYNAMICS
RESTART WAVEFUNCTION COORDINATES VELOCITIES LATEST
MIRROR
LSD
minimize the density fuctional before
QUENCH ELECTRONS BO
starting the MD simulation
ODIIS
E = min E [ρ(r), {R I }]
4 ρ(r)
TIMESTEP
5.0 Time step, in a.u.
EMASS (1 au = 0.024 fs, i.e. here we use 0.12 fs)
380.0

….
 Structure of the CPMD input file:
The CPMD section
&CPMD (cont.)
…
TEMPERATURE
300.0
TEMPCONTROL IONS
300.0 40.0 We keep the temperature at 300 ± 40 K
MAXSTEP
10000
TRAJECTORY SAMPLE XYZ
20
STORE
200
…

&END
 Structure of the CPMD input file
SECTIONS
&CPMD

&END

&SYSTEM

&END

&DFT

&END

&ATOMS

&END
 Structure of the CPMD input file:
The SYSTEM section

&SYSTEM
SYMMETRY a2 a3

1 cubic (periodic) box

CELL box enclosing the atoms a1

32.0000 1.0000 1.0000 0.0000 0.0000 0.0000

CUTOFF a1, a2/a1, a3/a1, cos α, cos β, cos γ
70.0 plane wave cutoff
(here cubic box of 32 a.u.**3 )
CHARGE
+0.0 total charge of the system
MULTIPLICITY (M)
3 triplet state (M = 2S+1)
&END
 Structure of the CPMD input file:
The SYSTEM section

&SYSTEM
ANGSTROM
SYMMETRY a2 a3

1 cubic (periodic) box

CELL box enclosing the atoms a1

15.0000 1.0000 1.0000 0.0000 0.0000 0.0000

CUTOFF a1, a2/a1, a3/a1, cos α, cos β, cos γ
70.0 plane wave cutoff
(here cubic box of 15 Å**3 )
CHARGE
0.0 total charge of the system
MULTIPLICITY (M)
3 triplet state (M = 2S+1)
&END
 INSET: BASIS SETS

Atom-centered basis set

Gaussian
function

Plane-wave basis set

Plane
waves
∑ ci (G) exp[iG.r]
ψi (r )=
1
ΩG
Atom-centered basis set

• Give chemical insight

• Depend on atomic position
• Basis set superposition error (BSSE)

Plane-wave basis set

• Less chemical insight
• Independent of atomic position
• No BSSE
∑ ci (G) exp[iG.r] Naturally periodic
ψi (r )=
1 •
ΩG • Many functions needed
• Forces calculation is easier
 &SYSTEM

HOW MANY PLANE WAVES DOES MY SYSTEM NEED?

(i.e. how many basis set functions?)

The number of PWs is defined by Ecutoff and the cell volume

order 3/2 with respect to Ecutoff

1
N PW ≈ 2 Ω E 3/ 2 [a.u.]
cutoff
2π The minimum Ecutoff
Number ot PWs needed depends on
Linear with respect to the cell volume the specific
pseudopotentials
used for each atom
Basis set size:

• Depends on the cell volume

• Does NOT depend on the number of atoms inside the cell

Basic memory depends on Nat, NSP, NPW

≈ (80 + 24 NSP + 4 NAT ) NPW
 &SYSTEM

Working with an isolated cell (not PBC)

… but in QM/MM
method for the solution of the Poisson
&SYSTEM equation for isolated systems ∇2 V (r ) = − 4πρ(r )
H
POISSON TUCKERMAN
SYMMETRY
0 isolated box or cell (also used for QM/MM. I.e. periodic boundary
conditions are not used or the QM region)
CELL
32.0000 1.0000 1.0000 0.0000 0.0000 0.0000
CUTOFF a1, a2/a1, a3/a1, cos α, cos β, cos γ
70.0 plane wave cutoff
(here cubic box of 32 a.u.**3 )
CHARGE
0.0 total charge of the system
MULTIPLICITY
3 a2 a3

&END
a1
Structure of the CPMD input file:
The DFT section

&DFT
FUNCTIONAL HCTH
&END
 Structure of the CPMD input file:
The ATOMS section
&ATOMS
*P_MT_SPR_HCTH.psp
LMAX=D LOC=S Pseudopotential specifications (information on
1 the nonlocality of the pseudopotential).
6.500 6.700 10.000 In general you do not need to change it from
*O_MT_HCTH.psp what is specified in the pseudopotential library
(README file)
LMAX=P LOC=P
2
7.5000 7.000 10.000
6.000 5.000 10.000
*….
CONSTRAINTS
FIX ATOMS
1
23
END CONSTRAINTS
&END
 PSEUDOPOTENTIAL SPECIFICATIONS
./PP_LIBRARY/[Link]

ELEMENT PP TYPE NAME USAGE REFERENCES

H LDA BHS H_BHS_LDA.psp LMAX=S LOC=S ECUT=25 Ry [1]

H LDA MT (Ballone) H_MT_LDA.psp LMAX=S LOC=S ECUT=25 Ry [2]
H BP MT (Boero) H_MT_BP.psp LMAX=S LOC=S ECUT=25 Ry [15]
H BLYP MT (Boero) H_MT_BLYP.psp LMAX=S LOC=S ECUT=25 Ry [3]
 OUTPUT ENERGIES
NFI EKINC TEMPP EKS ECLASSIC EHAM DIS TCPU
1 0.00001 299.6 -545.64887 -473.44754 -473.44754 0.184E-04 18.51
2 0.00010 298.5 -545.37030 -473.44686 -473.44676 0.741E-04 18.25
3 0.00048 296.6 -544.91694 -473.44599 -473.44551 0.173E-03 18.49
4 0.00138 294.1 -544.30697 -473.44522 -473.44384 0.326E-03 18.34
5 0.00301 291.0 -543.56475 -473.44483 -473.44181 0.551E-03 18.43
6 0.00547 287.5 -542.71975 -473.44498 -473.43951 0.866E-03 18.53
7 0.00874 283.7 -541.80538 -473.44575 -473.43701 0.129E-02 18.27
8 0.01273 279.7 -540.85764 -473.44715 -473.43442 0.186E-02 18.43
9 0.01733 275.8 -539.91361 -473.44915 -473.43182 0.258E-02 18.43
10 0.02239 272.0 -539.00979 -473.45151 -473.42912 0.347E-02 20.81

NFI = step number

EKINC = kinetic energy of the fictitious electronic subsystem (in a.u.)
TEMPP = temperature of the system (in K) (i.e. the nuclei)
EKS = Kohn-Sham energy (i.e. the electronic or potential energy)
ECLASSIC = total energy of the nuclei (= EKIN_nuclei + EKS)
EHAM = total energy of the CP system (= ECLASSIC + EKINC)
DIS = mean square displacement of atoms with respect to the first step (in a.u.2)

occ.
EHAM = ECLASSIC +2 ∑ ∫ µ ψ i (r ) 2dr
i
 Structure of the CPMD input file
SECTIONS
&CPMD

&END

&SYSTEM

&END
Other sections:
&DFT & PROPERTIES
& TDDFT
&END & CLASSIC
& VDW
&ATOMS & PIMD
…
&END
 Structure of the CPMD-QMMM input file
In the section &CPMD … &END we simply add the new keyword
QMMM, e.g.
&CPMD
MOLECULAR DYNAMICS
RESTART WAVEFUNCTIONS COORDINATES VELOCITIES LATEST
RESTART ACCUMULATORS LATEST
QMMM
ODIIS
4
…
&END

The sections &DFT … &END, &PROPERTIES … &END, etc..

do not change.
 Structure of the CPMD-QMMM input file
A new section &QMMM … &END must be included
&QMMM
A file that expresses how the force field describes
TOPOLOGY
[Link]
each of the “monomers” of the protein (residues,
COORDINATES ligands …) and their connectivity, the atomic point
[Link] charges etc.
INPUT Name of the file containing the protein coordinates
[Link]
in the appropriate format (gromos format)
AMBER
RCUT_NN
10.
RCUT_MIX force field used (AMBER or GROMOS)
15.
RCUT_ESP
25.
TIMING Gives the time for each part (QM, MM …)
ARRAYSIZES
MAXATT 28
…
&END
 Structure of the CPMD-QMMM input file
A new section &QMMM … &END must be included
&QMMM
TOPOLOGY
[Link]
COORDINATES
[Link]
INPUT
[Link] NN
AMBER
RCUT_NN
10.
RCUT_ESP set the cut-off radii for the 3 regions
25. r1 = 10, r2 = 25
…
&END
QM region
 The 3-regions scheme
NN
Region 1 (“Nearest Neighbor region”):
NN = subset of classical MM atoms inside this region
NN
ρ (x)
∑ qI ∫ d 3 x
I =1 x − RI
r =| x |< r1

Region 2 (“ESP region”):

Classical-RESP charges interaction: QM region

qJRESP ( ρ , R I )
∑ qI ∑
RI − RJ
r1 < r < r2
Note: r1, r2 are taken
I ⸦I ∈
ESP
NN J ∈QM
from any QM atom =>
Region 3 (“MM region”): the regions are not
Multipolar expansion on MM charges: perfectly spherical

℘α ( ρ )(x − R I )
α

∑q ∑
α
∈NN
I
x − RI
3
+ quadrupole r > r2
I ⸦IMM
See Raich et al. Methods Enzymology. 577, 159-183 (2016)
24
for a strategy to chose the radii of the 3 regions
 Structure of the CPMD-QMMM input file

The section &ATOMS … &END looks a bit different

&ATOMS

….

CONSTRAINTS
FIX MM You can fix the whole classical part if you wish
END CONSTRAINTS
&END
 Structure of the CPMD-QMMM input file

The section &ATOMS … &END looks a bit different

&ATOMS
*O_MT_HCTH.psp KLEINMAN-BYLANDER
LMAX=P LOC=P
40
421 435 441 443 444 445 -> give the atom number as listed in
… the [Link] file !
*C_MT_HCTH.psp KLEINMAN-BYLANDER
LMAX=P LOC=P
75
419 422 425 427 429 434 436 438 446 449
…
CONSTRAINTS
FIX MM
END CONSTRAINTS
&END
 OUTPUT ENERGIES IN CPMD
NFI EKINC TEMPP EKS ECLASSIC EHAM DIS TCPU
1 0.00001 299.6 -545.64887 -473.44754 -473.44754 0.184E-04 18.51
2 0.00010 298.5 -545.37030 -473.44686 -473.44676 0.741E-04 18.25
3 0.00048 296.6 -544.91694 -473.44599 -473.44551 0.173E-03 18.49
4 0.00138 294.1 -544.30697 -473.44522 -473.44384 0.326E-03 18.34
5 0.00301 291.0 -543.56475 -473.44483 -473.44181 0.551E-03 18.43
6 0.00547 287.5 -542.71975 -473.44498 -473.43951 0.866E-03 18.53
7 0.00874 283.7 -541.80538 -473.44575 -473.43701 0.129E-02 18.27

NFI = step number

EKINC = kinetic energy of the fictitious electronic subsystem (in a.u.)
TEMPP = temperature of the system (in K) (i.e. the nuclei)
EKS = Kohn-Sham energy (i.e. the electronic or potential energy)
ECLASSIC = total energy of the “classical system” (EKS + EKIN_nuclei)
EHAM = Hamiltonian energy. It is the total energy of the Car-Parrinello system, so
it includes the fictitious kinetic energy of the electrons (EHAM = ECLASSIC + EKINC)
DIS = mean square displacement of atoms with respect to the first step (in a.u.2)
occ.
EHAM = ECLASSIC +2 ∑ ∫ µ ψ i (r ) 2dr
i
 OUTPUT ENERGIES IN CPMD-QMMM
NFI EKINC TEMPP EKS ECLASSIC EHAM EQM DIS TCPU
1 0.00001 299.6 -545.64887 -473.44754 -473.44754 -299.64429 0.184E-04 18.51
2 0.00010 298.5 -545.37030 -473.44686 -473.44676 -299.65252 0.741E-04 18.25
3 0.00048 296.6 -544.91694 -473.44599 -473.44551 -299.66511 0.173E-03 18.49
4 0.00138 294.1 -544.30697 -473.44522 -473.44384 -299.68169 0.326E-03 18.34
5 0.00301 291.0 -543.56475 -473.44483 -473.44181 -299.70173 0.551E-03 18.43
6 0.00547 287.5 -542.71975 -473.44498 -473.43951 -299.72458 0.866E-03 18.53
7 0.00874 283.7 -541.80538 -473.44575 -473.43701 -299.74951 0.129E-02 18.27

NFI = step number

EKINC = kinetic energy of the fictitious electronic subsystem (in a.u.)
TEMPP = temperature of the system (in K) (i.e. the nuclei)
EKS = Kohn-Sham energy (i.e. the electronic or potential energy)
ECLASSIC = total “classical” energy of the system (= EKS + EKIN_nuclei )
EHAM = total energy of the CP system (= ECLASSIC + EKINC)
EQM = energy of the QM region (elecronic + kinetic energy). The QM atoms
are polarized by the MM environment
DIS = mean square displacement of atoms with respect to the first step (in a.u.2)
 OUTPUT ENERGIES IN CPMD-QMMM
NFI EKINC TEMPP EKS ECLASSIC EHAM EQM DIS TCPU
1 0.00001 299.6 -545.64887 -473.44754 -473.44754 -299.64429 0.184E-04 18.51
2 0.00010 298.5 -545.37030 -473.44686 -473.44676 -299.65252 0.741E-04 18.25
3 0.00048 296.6 -544.91694 -473.44599 -473.44551 -299.66511 0.173E-03 18.49
4 0.00138 294.1 -544.30697 -473.44522 -473.44384 -299.68169 0.326E-03 18.34
5 0.00301 291.0 -543.56475 -473.44483 -473.44181 -299.70173 0.551E-03 18.43
6 0.00547 287.5 -542.71975 -473.44498 -473.43951 -299.72458 0.866E-03 18.53
7 0.00874 283.7 -541.80538 -473.44575 -473.43701 -299.74951 0.129E-02 18.27
8 0.01273 279.7 -540.85764 -473.44715 -473.43442 -299.77581 0.186E-02 18.43
9 0.01733 275.8 -539.91361 -473.44915 -473.43182 -299.80291 0.258E-02 18.43
10 0.02239 272.0 -539.00979 -473.45151 -473.42912 -299.83032 0.347E-02 20.81

They refer to the whole system (QM and MM regions)

 OUTPUT ENERGIES IN CPMD-QMMM
NFI EKINC TEMPP EKS ECLASSIC EHAM EQM DIS TCPU
1 0.00001 299.6 -545.64887 -473.44754 -473.44754 -299.64429 0.184E-04 18.51
2 0.00010 298.5 -545.37030 -473.44686 -473.44676 -299.65252 0.741E-04 18.25
3 0.00048 296.6 -544.91694 -473.44599 -473.44551 -299.66511 0.173E-03 18.49
4 0.00138 294.1 -544.30697 -473.44522 -473.44384 -299.68169 0.326E-03 18.34
5 0.00301 291.0 -543.56475 -473.44483 -473.44181 -299.70173 0.551E-03 18.43
6 0.00547 287.5 -542.71975 -473.44498 -473.43951 -299.72458 0.866E-03 18.53
7 0.00874 283.7 -541.80538 -473.44575 -473.43701 -299.74951 0.129E-02 18.27
8 0.01273 279.7 -540.85764 -473.44715 -473.43442 -299.77581 0.186E-02 18.43
9 0.01733 275.8 -539.91361 -473.44915 -473.43182 -299.80291 0.258E-02 18.43
10 0.02239 272.0 -539.00979 -473.45151 -473.42912 -299.83032 0.347E-02 20.81

These energies refer to the QM region

 QM/MM simulation protocol
([Link] PDB structure

• add H atoms
• choose protonation states system
• add counterions preparation
• solvate
• run MM program (structure relaxation)

Relax structure (classical MD)

• run MM program (MD at constant T/P)

rmsd

Equilibrate structure (classical MD)

• prepare the QM/MM system, choose QM region

• run QM/MM program (structure relaxation)

Relax + equilibrate structure (QM/MM MD)

• run QM/MM dynamics QM/MM trajectory

 Specific for AMBER followed by CPMD QM/MM

Coordinates
1. System preparation Topology

2. Classical MD simulation ./amber_to_gromos

(convertor)

3. Format conversion
Gromos input file
New topology
4. QM/MM simulation New coordinates

cpmd input CPMD

cpmd output files