Using CP

This section is intended to explain how to perform basic Car-Parrinello (CP) simulations using the CP codes.

It is important to understand that a CP simulation is a sequence of different runs, some of them used to "prepare" the initial state of the system, and other performed to collect statistics, or to modify the state of the system itself, i.e. modify the temperature or the pressure.

To prepare and run a CP simulation you should:

1. define the system:
   1. atomic positions
   2. system cell
   3. pseudopotentials
   4. number of electrons and bands
   5. cut-offs
   6. FFT grids (CP code only)

2. The first run, when starting from scratch, is always an electronic minimization, with fixed ions and cell, to bring the electronic system on the ground state (GS) relative to the starting atomic configuration. Example of input file (Benzene Molecule):

    &control
       title = ' Benzene Molecule ',
       calculation = 'cp',
       restart_mode = 'from_scratch',
       ndr = 51,
       ndw = 51,
       nstep  = 100,
       iprint = 10, 
       isave  = 100,
       tstress = .TRUE.,
       tprnfor = .TRUE.,
       dt    = 5.0d0,
       etot_conv_thr = 1.d-9,
       ekin_conv_thr = 1.d-4,
       prefix = 'c6h6'
       pseudo_dir='/scratch/acv0/benzene/',
       outdir='/scratch/acv0/benzene/Out/'
    /
    &system
       ibrav = 14, 
       celldm(1) = 16.0, 
       celldm(2) = 1.0, 
       celldm(3) = 0.5, 
       celldm(4) = 0.0, 
       celldm(5) = 0.0, 
       celldm(6) = 0.0, 
       nat  = 12,
       ntyp = 2,
       nbnd = 15,
       nelec = 30,
       ecutwfc = 40.0,
       nr1b= 10, nr2b = 10, nr3b = 10,
       xc_type = 'BLYP'
    /
    &electrons
       emass = 400.d0,
       emass_cutoff = 2.5d0,
       electron_dynamics = 'sd',
    /
    &ions
       ion_dynamics = 'none',
    /
    &cell
       cell_dynamics = 'none',
       press = 0.0d0,
    /
   ATOMIC_SPECIES
    C 12.0d0 c_blyp_gia.pp
    H 1.00d0 h.ps
   ATOMIC_POSITIONS (bohr)
      C     2.6  0.0 0.0
      C     1.3 -1.3 0.0
      C    -1.3 -1.3 0.0
      C    -2.6  0.0 0.0
      C    -1.3  1.3 0.0
      C     1.3  1.3 0.0
      H     4.4  0.0 0.0
      H     2.2 -2.2 0.0
      H    -2.2 -2.2 0.0
      H    -4.4  0.0 0.0
      H    -2.2  2.2 0.0
      H     2.2  2.2 0.0
   

   You can find the description of the input variables in file INPUT_CP in the Doc/ directory. A short description of the logic behind the choice of parameters in contained in INPUT.HOWTO

3. Sometimes a single run is not enough to reach the GS. In this case, you need to re-run the electronic minimization stage. Use the input of the first run, changing restart_mode = 'from_scratch' to restart_mode = 'restart'.

   Important: unless you are already experienced with the system you are studying or with the code internals, usually you need to tune some input parameters, like emass, dt, and cut-offs. For this purpose, a few trial runs could be useful: you can perform short minimizations (say, 10 steps) changing and adjusting these parameters to your need.

   You could specify the degree of convergence with these two thresholds:

   etot_conv_thr: total energy difference between two consecutive steps

   ekin_conv_thr: value of the fictitious kinetic energy of the electrons

   Usually we consider the system on the GS when ekin_conv_thr < 10^-5 . You could check the value of the fictitious kinetic energy on the standard output (column EKINC).

   Different strategies are available to minimize electrons, but the most used ones are:

   • steepest descent:

       electron_dynamics = 'sd'
     

   • damped dynamics:

       electron_dynamics = 'damp', 
       electron_damping = 0.1,
     

     See input description to compute damping factor, usually the value is between 0.1 and 0.5.

4. Once your system is in the GS, depending on how you have prepared the starting atomic configuration, you should do several things:
   • if you have set the atomic positions ``by hand'' and/or from a classical code, check the forces on atoms, and if they are large ( 0.1 - 1.0 atomic units), you should perform an ionic minimization, otherwise the sistem could break-up during the dynamics.
   • if you have taken the positions from a previous run or a previous ab-initio simulation, check the forces, and if they are too small ( 10^-4 atomic units), this means that atoms are already in equilibrium positions and, even if left free, they will not move. Then you need to randomize positions a little bit. see below.

5. Minimize ionic positions.

   As we pointed out in 4) if the interatomic forces are too high, the system could "explode" if we switch on the ionic dynamics. To avoid that we need to relax the system.

   Again there are different strategies to relax the system, but the most used are again steepest descent or damped dynamics for ions and electrons. You could also mix electronic and ionic minimization scheme freely, i.e. ions in steepest and electron in damping or vice versa.

   1. suppose we want to perform a steepest for ions. Then we should specify the following section for ions:

       &ions
          ion_dynamics = 'sd',
       /
      

      Change also the ionic masses to accelerate the minimization:

      ATOMIC_SPECIES
       C 2.0d0 c_blyp_gia.pp
       H 2.00d0 h.ps
      

      while leaving unchanged other input parameters.

      Note that if the forces are really high (> 1.0 atomic units), you should always use stepest descent for the first relaxation steps ( 100 ).

   2. as the system approaches the equilibrium positions, the steepest descent scheme slows down, so is better to switch to damped dynamics:

       &ions
          ion_dynamics = 'damp',
          ion_damping = 0.2,
          ion_velocities = 'zero',
       /
      

      A value of ion_damping between 0.05 and 0.5 is usually used for many systems. It is also better to specify to restart with zero ionic and electronic velocities, since we have changed the masses. Change further the ionic masses to accelerate the minimization:

      ATOMIC_SPECIES
       C 0.1d0 c_blyp_gia.pp
       H 0.1d0 h.ps
      

   3. when the system is really close to the equilibrium, the damped dynamics slow down too, especially because, since we are moving electron and ions together, the ionic forces are not properly correct, then it is often better to perform a ionic step every N electronic steps, or to move ions only when electron are in their GS (within the chosen threshold).

      This can be specified adding, in the ionic section, the ion_nstepe parameter, then the ionic input section become as follows:

       &ions
          ion_dynamics = 'damp',
          ion_damping = 0.2,
          ion_velocities = 'zero',
          ion_nstepe = 10,
       /
      

      Then we specify in the control input section:

          etot_conv_thr = 1.d-6,
          ekin_conv_thr = 1.d-5,
          forc_conv_thr = 1.d-3
      

      As a result, the code checks every 10 electronic steps whether the electronic system satisfies the two thresholds etot_conv_thr, ekin_conv_thr: if it does, the ions are advanced by one step. The process thus continues until the forces become smaller than forc_conv_thr.

      Note that to fully relax the system you need many run, and different strategies, that you shold mix and change in order to speed-up the convergence. The process is not automatic, but is strongly based on experience, and trial and error.

      Remember also that the convergence to the equilibrium positions depends on the energy threshold for the electronic GS, in fact correct forces (required to move ions toward the minimum) are obtained only when electrons are in their GS. Then a small threshold on forces could not be satisfied, if you do not require an even smaller threshold on total energy.

6. randomization of positions.

   If you have relaxed the system or if the starting system is already in the equilibrium positions, then you need to move ions from the equilibrium positions, otherwise they won't move in a dynamics simulation. After the randomization you should bring electrons on the GS again, in order to start a dynamic with the correct forces and with electrons in the GS. Then you should switch off the ionic dynamics and activate the randomization for each species, specifying the amplitude of the randomization itself. This could be done with the following ionic input section:

    &ions
       ion_dynamics = 'none',
       tranp(1) = .TRUE.,
       tranp(2) = .TRUE.,
       amprp(1) = 0.01
       amprp(2) = 0.01
    /
   

   In this way a random displacement (of max 0.01 a.u.) is added to atoms of specie 1 and 2. All other input parameters could remain the same.

   Note that the difference in the total energy (etot) between relaxed and randomized positions can be used to estimate the temperature that will be reached by the system. In fact, starting with zero ionic velocities, all the difference is potential energy, but in a dynamics simulation, the energy will be equipartitioned between kinetic and potential, then to estimate the temperature take the difference in energy (de), convert it in Kelvins, divide for the number of atoms and multiply by 2/3.

   Randomization could be useful also while we are relaxing the system, especially when we suspect that the ions are in a local minimum or in an energy plateau.

7. Start the Car-Parrinello dynamics.

   At this point after having minimized the electrons, and with ions displaced from their equilibrium positions, we are ready to start a CP dynamics. We need to specify 'verlet' both in ionic and electronic dynamics. The threshold in control input section will be ignored, like any parameter related to minimization strategy. The first time we perform a CP run after a minimization, it is always better to put velocities equal to zero, unless we have velocities, from a previous simulation, to specify in the input file. Restore the proper masses for the ions. In this way we will sample the microcanonical ensemble. The input section changes as follow:

    &electrons
       emass = 400.d0,
       emass_cutoff = 2.5d0,
       electron_dynamics = 'verlet',
       electron_velocities = 'zero',
    /
    &ions
       ion_dynamics = 'verlet',
       ion_velocities = 'zero',
    /
   ATOMIC_SPECIES
   C 12.0d0 c_blyp_gia.pp
   H 1.00d0 h.ps
   

   If you want to specify the initial velocities for ions, you have to set ion_velocities = 'from_input', and add the IONIC_VELOCITIES
   card, with the list of velocities in atomic units.

   IMPORTANT: in restarting the dynamics after the first CP run, remember to remove or comment the velocities parameters:

    &electrons
       emass = 400.d0,
       emass_cutoff = 2.5d0,
       electron_dynamics = 'verlet',
       ! electron_velocities = 'zero',
    /
    &ions
       ion_dynamics = 'verlet',
       ! ion_velocities = 'zero',
    /
   

   otherwise you will quench the system interrupting the sampling of the microcanonical ensemble.

8. Changing the temperature of the system.

   It is possible to change the temperature of the system or to sample the canonical ensemble fixing the average temperature, this is done using the Nosè thermostat. To activate this thermostat for ions you have to specify in the ions input section:

    &ions
       ion_dynamics = 'verlet',
       ion_temperature = 'nose',
       fnosep = 60.0,
       tempw  = 300.0,
       ! ion_velocities = 'zero',
    /
   

   where fnosep is 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. tempw is the desired average temperature in Kelvin.

   It is possible to specify also the thermostat for the electrons, this is usually activated in metal or in system where we have a transfer of energy between ionic and electronic degrees of freedom.