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.[1] 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.[1] 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.[1]
# 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'scp.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'"
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)
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.