These examples illustrate the applications of cp.x MD code in the Quantum ESPRESSO distribution. Since MD simulations usually take very long, we will look at a small example, a single water molecule, to get familiar with all important options, features and the normal flow of execution for Car-Parrinello simulations with the cp.x code.

Getting ready

For the following tutorial you need to have the cp.x and cppp.x executables available. Download the exercise file examples_cp.tgz in a directory of your choice. Uncompress and unpack the archive and enter the created directory:

  tar -xzvf examples_cp.tgz  cd examples_cp  

Input File Structure

The input file for cp.x has a very similar structure as other programs in the Quantum ESPRESSO distrubution. The most notable differences are

1. the type of calculation is either cp or vc-cp. In both cases a time step dt must be specified (in Hartree atomic units).
2. no symmetry, no k-point list: only Gamma (k=0) is used
3. for ultrasoft pseudopotentials only: you need to specify a subset of the charge density grid (variables "nr1b", "nr2b", "nr3b") conatining the augmentation charges
4. there are two new flags, ndr=XX and ndw=YY. These are appended to the prefix value to form the directory for the restart files to read (=ndr) a restart from or write (=ndw) a restart to. If XX is the same as YY, then the old restart will be overwritten. With continuously incrementing these numbers you can preserve a history of your restart files and have the option to start from previous intermediate results.

Initial Wavefunction Optimization

Before we can start a Car-Parrinello MD, we have to compute the ground state wavefunction. This is achieved via direct minimization of the energy functional, not via slf-consistency + diagonalization. cp.x contains two algorithms to do this: steepest descend and damped dynamics. Both are very robust; steepest descent is notoriously inefficient, damped dynamics is efficient only with a good value of the damping. Since computing the initial wavefunction takes only a small computational effort compared to the ‘production‘ run, efficiency here does not matter that much. During wavefunction optimization, we update the wavefunction, but not the atom positions, thus we set in the &electrons namelist the keyword electron_dynamics=‘sd‘ (for steepest descent), but in the &ions namelist the keyword ion_dynamics=‘none‘. See input file h2o_mol1.in. You can monitor the progress of the calculation in the output file, or look at the energies in the h2o_mol.evp file. The dynamics is determined by a single parameter, dt²/(electron mass).

Mysterious errors in routine "ortho", performing iterative orthonormalization, are usually due to a too large time step. Especially at the beginning of the simulation, it may be useful to perform the first few time steps with Gram-Schmid orthogonalization (keyword ‘gram‘): it is slower but safer than the default ‘ortho‘ (iterative) orthonormalization

In a second step, we want to continue the optimization with the damped dynamics algorithm, so we set the value of restart_mode to ‘restart‘ and electron_dynamics to ‘damp‘. The damping parameter is not set and default to 0.1. See fileh2o_mol2.in. Now we would be ready to start a Car-Parrinello MD.

Electron Density and Atom Position Randomization

If we want to have a look at the electron density corresponding to the starting struction, we have to re-run the wavefunction optimization with the keyword disk_io=‘high‘ (normally the density is not needed and thus not written), see file h2o_mol.cp-dens.in.

Since we were starting from optimized positions of the water molecule, we also want add some potential energy to the system by randomizing the atom positions (we keep the oxygen fixed here) a little bit. So we use in the &ions section the keywords:

  tranp(1) = .FALSE.  tranp(2) = .TRUE.  amprp(2) = 0.1

meaning "randomization of the second atom type only" and "maximum random displacement 10% of the lattice constant".

After this run we have a restart with the re-optimized wavefunction for the randomized atom positions and a file with the charge density in the restart directory. To visualize the charge density, we first need to run the post-processing program, cppp.x to create a .xsf file. The corresponding input file, h2o_mol.cppp-dens.in is pretty self explanatory. Note that both the value of prefix and the value of the ndr keyword has to match the prefix and ndwvalue of the calculation that we want to postprocess. For a list of the available postprocessing options, see the INPUT_CPPP file in the Doc directory of the Quantum ESPRESSO distrubution.

Structural Optimization

Damped dynamics can be used in structural optimization as well. File h2o_mol3.in shows how to find the minimum energy structure using damped dynamics for both ions and electrons. There are two damping parameters, one for electrons, one for ions. As a rule of thumb, the latter should be an order of magnitude smaller than the former. If the convergence becomes sluggish (and it typically does close to the minimum) reduce the damping.

Car-Parrinello Molecular Dynamics.

Now we are ready to start a Car-Parrinello MD: see input file h2o_mol4.in. We don‘t want the outputs from the previous runs of dynamics, so we first clean up the directory where files are written:

   rm -f h2o_mol.???

(don‘t remove the h2o_mol_save directory). Note that we now restart with restart_mode=‘reset_counters‘. This is like a regular restart, but all averages and counters are reset. We also increment ndw, so we do not overwrite the initial wavefunction and have the option to go back in case of problems.

Note that the input data for this run are read from step 2, not 3, i.e. the distorted H₂O molecule, not the equilibrium one. The atomic structure is read from file as well when you restart from file: what is written in input dat is ignored.

For regular CP-dynamics, we set both electron_dynamics and ion_dynamics to ‘verlet‘. We also initialize both the electronic and ionic velocities to 0. We now have to set a fictitious mass for the CP-hamiltonian and have to make sure that this mass and the time step are properly chosen, so we conserve the total energy and stay close to the Born-Oppenheimer surface. This is best done by looking at the evolution of the energies in the file h2o_mol.evp. Column 1 is the time step number, column 2 the fictitious kinetic energy (in atomic units), column 4 the kinetic energy of the system (in Kelvin), column 5 the potential energy of the system, and column 7 contains the total energy of the classical system and column 8 the conserved quantity of the CP-hamiltonian. These are affected by the proper choice of fictitious mass and time step. Play around with these numbers, to see how bad it can go, or how you can improve on the current choice of parameters.

A too large time step will manifest itself as a drift of the constant of motion, column 8. A too heavy electron mass, or a too small energy gap in the system, may cause a loss of adiabaticity, leading to an energy transfer from ionic to electronic degrees of freedom. A more subtle effect of an excessive electron mass is a "drag" effect on ions that spoils the quality of the ionic trajectories.

h2o_mol5.in continues the trajectory by restarting from the previous run. Note that we have to go the restart_mode=‘restart‘ and remove the quench of electronic and ionic velocities. h2o_mol6.in finally restarts the trajectory with a Nose-Hoover thermostat controlling the ionic kinetic energy. This will oscillate around the prescribed temperature ‘tempw‘ (43K in this case). The important parameter here is ‘fnosep‘, the frequency of the thermostat in THz. This should be chosen to be comparable with the center of the vibrational spectrum of the system, in order to excite as many vibrational modes as possible. Don‘t be scared by the large oscillations around the equilibrium temperature: Nose‘ thermostats do not impose a fixed temperature, they add or remove energy from the system so that on the average the temperature is constant.

Quantum ESPRESSO tutorial: Car-Parrinello Molecular Dynamics