# The Car-Parrinello method

### Introduction

In this short tutorial I will show how I did some of my Car-Parrinello simulations, assuming that you know the theory. This is a result of many months of work (trial and error) done by Riccardo Bertossa and Federico Grasselli during the years 2017 and 2018 at SISSA, Trieste (Italy).

The main concepts of the method can be found in the original article (why don't we use a Lagrangian?), here (why the method works?) and here (you must be careful about the fictitious mass of the electron!). Oversimplifying the method: suppose you have calculated a DFT ground state, given the position of the atoms. Now think of this ground state electronic density as a fluid, as the water in a bowl of your kitchen, and identify the bowl with our atomic potential. We can associate the DFT ground state, where the energy is at minimum, to the water-bowl system at rest: also in this case the energy is at the minimum. If I move a little bit the bowl, the water will start moving, following some classical Lagrangian law of motion. Now suppose that I have many bowls that interact with each other, and the interaction depends on the shape of the water surface in the bowls. They will start to move, and then also the water will follow! If the mass of the water (that is the fictitious mass of the electrons) is in the correct range and the movements of the bowls are not too fast, the shape of the surface will not change too much with respect to the system at rest, the calculated forces (that remember, depends on the shape of the water surface) will not be too wrong and the water will not go out of the bowls! In this case I can run a classical dynamic simulation, without solving a self-consistent equation at each time step. The electronic wave functions and the atomic positions will evolve together. And we have for free a conserved quantity: the energy of this fictitious system.

So we have many ingredients. The standard DFT ones (pseudo-potentials, plane waves, cutoffs, ...) that defines our DFT energy functional. This is a functional of the electronic density and the atomic positions. Then we have the Car-Parrinello Lagrangian, with the associated fictitious mass of the electronic degrees of freedom. Let's see how we can do a molecular dynamic simulation with this!

### Setup of the software

Before starting, we need to setup the software on a UNIX system (it is not mandatory, but I don't know how to do this on Windows). In this tutorial I will use Quantum Espresso and LAMMPS. You need to have installed on your (super)computer the MPI library, a compatible FORTRAN compiler and a C++ one.

Download Q.E. here. Unpack it. Go inside the distribution directory and open a terminal. Type:
`./configuremake cp`,
Then hope that everything goes well. Check that the MPI library is detected in the output of `./configure`. If not, install it on the system. Note that `make` can be executed in parallel for example with `make -j12` for a faster compilation (here -j12 means: compile up to 12 files in parallel). In the quantum espresso distribution this does not always work as expected, while the serial compilation is always successful. You will end up with the executable `CPV/src/cp.x`. Try to run it and see what happens (Ctrl + C to terminate the code).

Download LAMMPS here. Go inside the distribution directory and open a terminal. Type:
`mkdir build; cd buildcmake -D PKG_KSPACE=yes ../cmakemake -j12`
Check again in the output of the `cmake` command (on some systems it is `cmake3`) that MPI was detected. `-D PKG_KSPACE=yes` tells cmake to compile the routines that are needed to compute the Ewald sums. The program is located in the current directory. Try to run `./lmp` and see what happens! (Do you remember the command to terminate the application?)

### Initial conditions?

The Car-Parrinello molecular dynamic is tricky to start. Think about it: you have to solve the newton equation of motion for a zillion of variables. To start the numerical integration you have to give the velocities of the atoms, the positions of the atoms, the initial wave functions (that are given by a minimization procedure), and the velocities of the wave functions. Every point in the space has to have a initial value for the wave functions and a starting velocity.

Then, suppose you want to simulate a liquid. How to choose the initial positions? If you do a easy random choice there is a little chance that everything will work smoothly, as in classical simulations, but here is worse. Remember the story of the water and the bowls? What will happen to the water if a bowl is moved too fast? And then the forces (that depends on the shape of the water)? For sure you will get very expensive random number as soon as the second time step is computed.

You have many choices to solve these problems: start with a more stable Born-Oppenheimer dynamics (solve the self-consistent DFT equation at every time step) and slowly go to the configuration that you need, or, if you have a cheap force field for your system, use it to evolve an other system into a thermalized state and hope that this state is good, after few Born-Oppenheimer steps, for your Car-Parrinello (I will call it CP from now on) simulation. This will solve (maybe) the problem for the atoms, but what about the velocity of the wave functions? You can say: easy, let's put them to zero at the beginning! Unfortunately this is not a good idea. Remember the water and the bowl? If the bowl is moving, and the water is not, then when the dynamic starts the water will receive a little kick. It is better if the water moves from the beginning with the correct velocity associated to the motion of the bowl. In practice, a single electronic ground state is sufficient to start a CP run, but the simulation will work better if you calculate the correct starting velocity of the wave functions. It is as easy as taking the difference of two ground states of two consecutive time steps (with something more that I will not talk about, see Cold restart of Car–Parrinello molecular dynamics).

OK, it is a mess. Let's see a example of a cheap procedure of getting the initial condition for a molten salt: KCl.

### Initial atomic position and velocities

So let's start with a simple system: 4 atoms of Cl and 4 atoms of K, using the Fumi-Tosi classical potential (very simple and cheap). The LAMMPS input script will look something like this:

`````` # Molten KCl simulation; R. Bertossa e F. Grasselli, 13 ottobre 2017
echo both
# User defined variables
# Fundamental constants
variable     kB         equal  1.3806504e-23    # [J/K]
variable     kB_eV      equal  8.617343e-5
variable     NA         equal  6.02214e23       # [1/mol]
variable     massunit   equal  1.660538921e-27  # [kg]
variable     V          equal  vol
variable     eV         equal  1.60218e-19      # eV in Joules
variable     dt         equal  0.001            # [ps] -- MD time step
#  KCl Potential; from Ref. dx.doi.org/10.1021/jp5050332 | J. Phys. Chem. B 2014, 118, 10196−10206
#  1 = +
#  2 = -
variable        b       equal  0.338e-19/\${eV}
print           \${b}
variable        A11     equal  1.25*\$b
variable        A12     equal  1.0*\$b
variable        A22     equal  0.75*\$b
variable        rho     equal  0.3367   # in Angstrom
variable        sigma11 equal  2.926    # in Angstrom
variable        sigma12 equal  3.048    # in Angstrom
variable        sigma22 equal  3.170    # in Angstrom
#
variable        C11     equal  24.3e-79/\${eV}/1e-60   # eV Angstrom^6
print           \${C11}
variable        C12     equal    48.0e-79/\${eV}/1e-60   # eV Angstrom^6
variable        C22     equal   124.5e-79/\${eV}/1e-60   # eV Angstrom^6
# note: In Ref. are defined with opposite sign with respect to LAMMPS
variable        D11     equal   -24.0e-99/\${eV}/1e-80   # eV Angstrom^8
print           \${D11}
variable        D12     equal   -73.0e-99/\${eV}/1e-80   # eV Angstrom^8
variable        D22     equal  -250.0e-99/\${eV}/1e-80   # eV Angstrom^8
variable        alat    equal     6.292 # Angstrom; initial guess; from Wikipedia
# Simulation cell setup
dimension       3
units           metal
boundary        p p p
atom_style      charge
lattice         fcc \${alat} origin 0 0 0        # Lattice for K system
region          box block 0 1 0 1 0 1 units lattice
create_box      2 box # "2": atom-type number, "box": name of the box
create_atoms    1 box
lattice         fcc \${alat} origin 0.5 0 0      # Lattice for Cl system
create_atoms    2 box
# Masses and charges
mass 1 39.098   #  K
mass 2 35.453   #  Cl
set type 1 charge +1
set type 2 charge -1
# Ewald sums
kspace_style ewald 1.0e-5  # as in Ref.
# Fumi-Tosi potential (Born-Huggins-Mayer)
pair_style born/coul/long 13.0
pair_coeff 1 1 \${A11} \${rho} \${sigma11} \${C11} \${D11}
pair_coeff 1 2 \${A12} \${rho} \${sigma12} \${C12} \${D12}
pair_coeff 2 2 \${A22} \${rho} \${sigma22} \${C22} \${D22}

#now we minimize the energy to find the lattice parameter with this potential
#I want that the simulation box is able to expand during the minimization
fix     1 all   box/relax iso 0.0 vmax 0.01
#"1": name of the fix, "all": every atom is affected
#"vmax 0.01": maximum delta V / V in the minimization, where V is the current volume of the simulation box
#"iso 0.0": target pressure

min_style       cg #Conjugate gradient for the minimization.
#	do it! ( https://lammps.sandia.gov/doc/minimize.html )
minimize        0.0 1.0e-15 1000000 1000000
unfix 1
#initialize velocities for the MD run
velocity        all create 1000 87287 mom yes rot yes
timestep        \${dt}
reset_timestep  0
# Langevin thermostat
#

thermo 50
thermo_style	custom step time etotal pe temp press lx
dump 		FILM all atom 50 KCLequil.lammpstrj

variable	lange_int	equal	\${dt}*1.0e2
variable	temp_min	equal	1000
variable	temp_max	equal	2800
variable	temp_equ	equal	1200
variable	npt_int		equal	\${dt}*1.0e2

fix		NPT all npt temp \${temp_min} \${temp_max} \${npt_int} iso 0 0 \${npt_int}
run		10000
unfix		NPT
fix		NVE all nve
fix 		LANGEVIN all langevin \${temp_max} \${temp_max} \${lange_int} 699483
run		30000
unfix		LANGEVIN
fix 		LANGEVIN all langevin \${temp_max} \${temp_equ} \${lange_int} 699483
run		10000
unfix		LANGEVIN
unfix		NVE

#zero c.m. velocity and angular momentum
velocity all zero linear
velocity all zero angular

#some Nose-Hoover

fix		NPT all npt temp \${temp_equ} \${temp_equ} \${npt_int} iso 0 0 \${npt_int}
run		10000
unfix		NPT

variable        npt_int         equal   \${dt}*1.0e3

fix		NPT all npt temp \${temp_equ} \${temp_equ} \${npt_int} iso 0 0 \${npt_int}
run		10000
unfix		NPT

variable        npt_int         equal   \${dt}*1.0e4

fix		NPT all npt temp \${temp_equ} \${temp_equ} \${npt_int} iso 0 0 \${npt_int}
run		10000
unfix		NPT

fix		NVE all nve
run 		10000
unfix		NVE

write_dump all custom 'KCl-LAMMPStoQE-VELOCITIES.dat' element vx vy vz modify element K Cl
write_dump all custom 'KCl-LAMMPStoQE-POSITIONS.dat' element xs ys zs modify element K Cl
``````

It is commented, but if you don't understand the script, I suggest you to read the LAMMPS documentation. In general a LAMMPS simulation is ran with `mpirun`, but in this case the simulation is so small that it is more efficient to run it with a single process. So it is sufficient to run:
`lmp -in input_file_that_you_see_above.in`
and wait for the (very fast) result.

### Initial wave functions and wave functions velocities

Now that we have a starting configuration for the atomic degrees of freedom, let's prepare the input for Quantum Espresso's `cp.x` program, so we can calculate the starting generalized coordinates for the electronic degrees of freedom! The script is the following (it generates 2 inputs: one for starting the simulation and one for running it):
`````` #!/bin/sh

echo "This example generates the script for finding the"
echo "electronic GS in CP from a given lammps configuration"
echo "of bulky KCl."

#Be sure to modify the following line
PSEUDO_DIR="/your/path/to/the/pseudopotentials/directory"

TMP_DIR=./tmp-cgstart
ECUT=50
MYPREFIX="KCl.\$ECUT"
MYINPUTFILE="\$MYPREFIX-cgstart.in"
MYINPUTFILE2="\$MYPREFIX-cp.in"
#LAMMPS output
LAMMPS_IN_POS="KCl-LAMMPStoQE-POSITIONS.dat"
LAMMPS_IN_VEL="KCl-LAMMPStoQE-VELOCITIES.dat"
CONV_VEL=2.1876912633e4 # Angstrom/ps

NAT=\$(grep -A 1 'NUMBER OF ATOMS' \$LAMMPS_IN_POS | tail -n 1 )

#lattice parameter in bohr radius.
#If not sure check the output of LAMMPS and put it there!
MYALAT=\$(echo \$(echo \$(tail -n 100 log.lammps | head -n 1 | awk '{print \$7}')/0.5292 | bc -l | cut -c 1-6 ) )

echo "!!!!!!!!!!!!!!!!!!!!"
echo "CHECK THE FOLLOWING!"
echo \$NAT atoms
echo alat \$MYALAT Bohr radius
echo ecutwfc \$ECUT Ry cutoff energy
echo "!!!!!!!!!!!!!!!!!!!!"

cat > \$MYINPUTFILE << EOF
&CONTROL
title = ' Potassium Cloride ',
calculation = 'cp',
restart_mode = 'from_scratch',
ndr = 70,
ndw = 71,
verbosity = 'medium',
prefix = '\$MYPREFIX',
outdir = '\$TMP_DIR',
pseudo_dir = '\$PSEUDO_DIR',
tprnfor = .true.,
tstress = .true.,
nstep = 30,
max_seconds = 42500,
dt = 0.5d0,
etot_conv_thr = 0.5d-5,
ekin_conv_thr = 0.5d-4
/
&SYSTEM
ibrav = 1,
celldm(1) = \$MYALAT,
nat = \$NAT,
ntyp = 2,
ecutwfc= \$ECUT,
/
&ELECTRONS
emass = 400.d0,
emass_cutoff = 2.5d0,
orthogonalization = 'Gram-Schmidt',
electron_dynamics = 'cg',
/
&IONS
ion_dynamics = 'verlet',
ion_velocities = 'from_input',
/

&CELL
cell_dynamics = 'none',
press = 0.0d0
/

ATOMIC_SPECIES
K   39.098 K_ONCV_PBE-1.0.upf
Cl  35.453 Cl_ONCV_PBE-1.0.upf

AUTOPILOT
ON_STEP = 2 :  dt = 1.0d0
ON_STEP = 10 : dt = 2.0d0
ON_STEP = 20 : dt = 5.0d0
ENDRULES

ATOMIC_POSITIONS {crystal}
\$(tail -n \$NAT \$LAMMPS_IN_POS)

ATOMIC_VELOCITIES { a.u }
\$(tail -n \$NAT \$LAMMPS_IN_VEL | awk "{ print \\$1, \\$2/\$CONV_VEL/\$MYALAT, \\$3/\$CONV_VEL/\$MYALAT, \\$4/\$CONV_VEL/\$MYALAT }" )

EOF

cat > \$MYINPUTFILE2 << EOF
&CONTROL
title = ' Potassium Cloride ',
calculation = 'cp',
restart_mode = 'restart',
ndr = 71,
ndw = 72,
verbosity = 'medium',
prefix = '\$MYPREFIX',
outdir = '\$TMP_DIR',
pseudo_dir = '\$PSEUDO_DIR',
tprnfor = .true.,
tstress = .true.,
nstep = 50,
max_seconds = 42500,
dt = 5.0d0,
etot_conv_thr = 0.5d-5,
ekin_conv_thr = 0.5d-4
/
&SYSTEM
ibrav = 1,
celldm(1) = \$MYALAT,
nat = \$NAT,
ntyp = 2,
ecutwfc= \$ECUT,
/
&ELECTRONS
emass = 400.d0,
emass_cutoff = 2.5d0,
orthogonalization = 'ortho',
electron_dynamics = 'verlet',
/
&IONS
ion_dynamics = 'verlet',
ion_velocities = 'default',
/

&CELL
cell_dynamics = 'none',
press = 0.0d0
/

ATOMIC_SPECIES
K   39.098 K_ONCV_PBE-1.0.upf
Cl  35.453 Cl_ONCV_PBE-1.0.upf

ATOMIC_POSITIONS {crystal}
\$(tail -n \$NAT \$LAMMPS_IN_POS)

ATOMIC_VELOCITIES { a.u }
\$(tail -n \$NAT \$LAMMPS_IN_VEL | awk "{ print \\$1, \\$2/\$CONV_VEL/\$MYALAT, \\$3/\$CONV_VEL/\$MYALAT, \\$4/\$CONV_VEL/\$MYALAT }" )

EOF

echo your first input file for cp.x is "'\$MYINPUTFILE'"
echo your second input file for cp.x is "'\$MYINPUTFILE2'"

``````

Be careful to the parameters of your system! In particular the lattice parameter and the number of atoms. Make sure that `ecutwfc` is converged. You can do it with normal `pw.x` calculations as you know (if you don't know this, look at other tutorials). Download the ONCVPSP pseudo-potentials here, extract them from the archive and set the directory in the script. Execute the script. Now let's run `cp.x`! Yeah! Note: I assume that the executable is in your path
`mpirun -np 12 cp.x -in KCl.50-cgstart.in > KCl.50-cgstart.out `
where `-np 12` means 'hey MPI library, start 12 processes'. Typically to have a number of processes equal to the number of cores of your computer is a good idea. `KCl.50-cgstart.in` is the fresh input file generated by the bash script that you can find above, and ` > KCl.50-cgstart.out` means 'hey bash, put all the output of the program inside the file KCl.50-cgstart.out, overwriting it'. This run will do 30 steps (did you see the nstep variable in the input?) using the Verlet algorithm to update the atomic positions and velocities, while the electronic degrees of freedom will be calculated with a conjugate gradient minimization at every time step (expensive). You can say... Where the hell is the wave function velocity calculated? I can't see anything in the input script! But... Use the Code Luke! In the complicated CG routine there is a flag that allows the code to calculate the wave function velocity, by doing two CG minimization and a projection without telling you anything, so we are safe! At the end of the calculation in the restart file the wave function and its velocity will be present, together with the atomic positions and velocities. Hopefully the system will be in a state that allows us to finally run the Car-Parrinello equations of motion, so...

## Run it now!

`mpirun -np 12 cp.x -in KCl.50-cp.in > KCl.50-cp.out `
Done! What am I doing?

### Oh no, why is it exploding?

Reading the article of Pastore, Smargiassi and Buda can be a good starting point. In few words the method works only when the electronic degrees of freedom are decoupled from the atomic one, and there is no flow of kinetic energy from the atoms to the electrons.

By the way, did you notice all those columns in the output of the code? The most important are:

• ekinc : the kinetic energy of the electronic degrees of freedom (how much does the electronic fluid move inside the potential)
• tempp : temperature of the nuclei
• etot : DFT energy (the potential energy of this messy equation of motion)
• econt : constant of motion of the Car-Parrinello Lagrangian (no physics in this number!)
• econs : econt - ekinc (the physical energy of the system)
If you want to do a good simulation you want a low and constant ekinc. This means that the electronic fluid is doing more or less always the same movements, and it will never go out of the bowl. When the coupling between the electrons and the ions becomes significant, you will see an exponential divergence of ekinc, and the simulation will not be physically significant at all. To have a small coupling between the electronic and ionic system you need very different masses in the two systems.

There is an other possible reason that transforms your simulation in expensive low quality random numbers: the time step. The time step must be big enough so you don't wait forever the end of the computation and small enough to account for the fastest movements of the electronic degrees of freedom. You can check that the time step is OK by looking at the behavior of the conserved quantity econt, that is a constant in the limit of the time step going to zero. To use a bigger time step you need a bigger mass for the fastest generalized coordinate.

### Finding the best parameters!

We have to find a compromise between many competing effects in the simulation. Luckily, we have two parameters to play with: the time step and the fictitious mass of the electronic degrees of freedom. Wait. In principle every point in the space can have a different mass, because it is artificial, this mass does not exist in the real world. In practice this is done in the reciprocal space: the mass of the electron wave function depends on his k component (related to the multigrid method). This is controlled by the parameter `emass_cutoff`.

But does the artificial mass of the electrons produce other effects? Yes: this mass contributes to the inertia of the atoms (Tangney). We have an other actor in this game: the systematic error due to the increased inertia of the nucleai.

So what can you do? Start with a very small time step (1 a.u.) and a small wave function mass (say 25 electron masses). Compare the calculated forces with a plain DFT self consistent calculation for some of the snapshots. How much, on average, it is wrong? (remember that the CP force should oscillate around the 'true' one) Increase the mass. Continue until you like the systematic error. Then start to play with the time step (you can use AUTOPILOT and pilot.mb). If you don't like it, you can try some different emass_cutoff, and/or repeat everything. It is up to you!

Good luck!

### Notes on the Quantum Espresso's cp.x code

The documentation of the `cp.x` code was not always correct, and there were around some bugs that have been corrected. One of the tricky things that you can find if you use an old version of the documentation is that the units of measure of the velocity were wrongly stated. They are not in atomic units, but in the same units of the positions divided by time in atomic units (documentation corrected in version 6.3). Note that in version 6.3 was corrected the change of time step with the AUTOPILOT module (the old code did not rescale the velocities if you change dt with AUTOPILOT). If you use older versions, you have to manually do a restart and use the option tolp (read the documentation). In this version is added also the rescaling of the electronic degrees of freedom when changing time step (before it was missing), making the transition smoother. If you use a restart file you have to specify a new option to use this feature, while with autopilot is done automatically. In this version of the code the AUTOPILOT module has the ability to switch from Verlet to cg the dynamic of the electrons on the fly, with some limitation on the parallelization (the conjugate gradient is parallelized only on plane waves).

This AUTOPILOT... I think it is very useful, at least at the beginning of the simulation, but its documentation file is hidden in the `CPV/Doc` directory of the distribution (here). It is also very easy to add a new autopilot feature, because it is well documented and commented in the code (use the Code Luke!).

Changing a variable on the fly (when the simulation is running!) is as easy as placing a file called `pilot.md` with a content that is easy to understand. For example:
`PILOT NOW : ISAVE = 50 NOW : DT = 5.0 NOW+50 :IONS_TEMPERATURE='damped' NOW + 150: TEMPW=350.0 ENDRULES `
or, if you want to do a step of CG with the 'cold restart', you can put at the end of the input the following:
` AUTOPILOT on_step = 10000 : orthogonalization = 'Gram-Schmidt' on_step = 10000 : electron_dynamics = 'cg' on_step = 10001 : electron_dynamics = 'verlet' on_step = 10001 : orthogonalization = 'ortho' ENDRULES `
or in pilot.mb as in the previous example. I encourage you to read the documentation.