Car-Parrinello Molecular Dynamics

    Introduction

    The standard procedure to perform an ab-initio, Car-Parrinello molecular dynamics simulation consists of a number of steps that have to be executed in the correct order. The first step consists in preparing a statistically meaningful initial configuration.
    In the following example we will prepare the initial configuration of a thermally disordered Si crystal by randomly displacing the atoms from their ideal crystalline positions. We will start by considering the ideal Si crystal in a cubic supercell containing 8 atoms, and finding the corresponding electronic ground state. We will then impart small random displacements to the atoms, and will let the MD run start. This simulation (8 Si atoms in cubic supercell) was the first application of the Car-Parrinello method, and is described in the original paper of R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).

    Electronic minimization of an ideal Si crystal

      Download the sample input file "si.in.00" here. Most of the input parameters are already familiar to us. We will focus on the parameters that are new, i.e. specific to molecular dynamics.

      Molecular dynamics generally consists of several simulations carried out in a sequence, so at the end of each simulaton it is typically useful to store all information necessary to continue a calculation (wfc's, positions, etc) in a separate directory, which is by convention labeled with a number between 50 and 99. The input variable "ndw" labels the directory where the program writes the final configuration, 
              ndw = 50
      The parameter "dt" determines the time step for the integration of the MD equations. It is expressed in atomic units (1 a.u. = 2.4 10^-17 s). In our case, we set 
              dt = 10.0d0
      "nstep" fizes the total number of time steps that will be executed 
              nstep  = 300
      The parameter "iprint" determines after how many MD steps the codes prints out a summary of several interesting variables, like atomic positions, forces, velocities, etc. 
              iprint = 10
      Notice that a number of quantities are printed in the standard output, but the most relevant ones are also printed out in separate files for different physical quantities. Positions are printed out in file "prefix".pos, forces in "prefix".for, various energies in "prefix".evp, and so on. Notice that the code appends (and does not overwrite) the output of consecutive runs to the same .pos, .for, etc, files.
      The directory "ndw" (see above) can be saved not only at the end of the run, but also every "isave" steps. This is to make sure that if the simulation terminates abruptly (e.g. power failure) you can restart it from the last saved step. 
              isave  = 100
      Before we can start a MD simulation we must ensure that the electrons are in their ground state. In this example, we will use a steepest descent ("sd") algorithm to find the electronic ground state 
              electron_dynamics = 'sd'
      The electron dynamics is controlled by the "fictitious" mass "emass", 
              emass = 400.d0
      The ions will not move while we determine the ground state, 
              ion_dynamics = 'none'
      We place the atoms in the ideal positions of a face centered cubic conventional cell, 
              ATOMIC_POSITIONS (crystal)
          Si     0.0000    0.0000    0.0000  1 1 1
          Si     0.0000    0.5000    0.5000  1 1 1
          Si     0.5000    0.5000    0.0000  1 1 1
          Si     0.5000    0.0000    0.5000  1 1 1
          Si     0.2500    0.2500    0.2500  1 1 1
          Si     0.2500    0.7500    0.7500  1 1 1
          Si     0.7500    0.7500    0.2500  1 1 1
          Si     0.7500    0.2500    0.7500  1 1 1
      We are now ready to execute the electronic minimization: 
              ../bin/cp.x < si.in.00 > si.out.00
      Notice that the resulting forces on the ions are zero, due to symmetry. This is not a convenient starting point for a MD simulation because the ions will not move away from this configuration.


    Random displacement of the atoms from their ideal positions

      After having determined the electronic ground state of the ideal system, we now diplace the atoms randomly from their ideal positions. Let's construct the input file for the next calculation ("01"). Copy the previous file into a new file 
              cp  si.in.00  si.in.01
      and modify the new file "si.in.01" as follows. We need to read the old configuration from the stored unit "50", so we need to change the restart mode from "from_scratch" to 
              restart_mode = "restart"
      and we need to add the line 
              ndr = 50
      If we want to keep the information (positions, wfc's, etc) of the "ideal" crystal, we better write on a different unit, so we need to change the value of "ndw" as 
              ndw = 51
      Notice that unit "50" will not be touched if ndr.ne.ndw.
      The randomization is obtained by adding to the "&ions" namelist the following parameters: 
              tranp = .true.        
        amprp = 0.1d0
      All ionic positions will be displaced with respect to the positions read from unit "ndr", by a random quantity of the order of "amprp" (in atomic units).
      The code will now displace the ions randomly, and perform an electronic minimization on the distorted configuration: 
              ../bin/cp.x < si.in.01 > si.out.01


    Let's dance!

      We are now ready to perform a MD simulation. Let's construct the new input file ("02"). Copy the previous file into a new file 
              cp  si.in.01  si.in.02
      and modify the new file "si.in.01" as follows.
      For "safety" reasons, it is convenient to update the read-write units as 
              ndr = 51 
        ndw = 52
      It is also essential to remove the randomization parameters (we don't want to randomize the ions again, they are already displaced). The "tranp" and "amprp" parameters must be removed.
      Both electrons will now move according to a second-order (Verlet) dynamics, and we want them to start with zero velocities, so we need to set 
              electron_dynamics = 'verlet'
        electron_velocities = 'zero'
      in the "&electrons" namelist, and 
              ion_dynamics = 'verlet'
        ion_velocities = 'zero'
      in the "&ions" namelist. Finally, we may want to change the total number of steps 
              nstep  = 1000
      With a time step of 10 a.u., this implies that we will be simulating the dynamics of the system for a total time of 0.24 ps. For convenience, we may also want to reset the time step counter 
              restart_mode = 'reset_counters',
      We also clean up the directory by removing files that we now want to overwrite: 
              rm si.pos si.cel si.str si.evp si.for si.vel si.spr si.the si.nos si.eig  si.con si.pol
      The code will now perform a Car-Parrinello MD simulation configuration: 
              ../bin/cp.x < si.in.02 > si.out.02
      We will consider here two typical analyses of a MD trajectory. First, we will learn how to visualize the trajectories using "xcrysden". We will then also learn how to compute pair-correlation functions.


    Visualizing the trajectory with XCrysDen

      Visualizing the MD trajectory using XCrysDen requires a post-processing tool that you need to download. Download the post-processing source code cp2xsf.f90 and and the sample input file cp2xsf.in. Compile xsf.f 
              gfortran cp2xsf.f90 -o cp2xsf.x
      Take a look at "cp2xsf.in". It contains information about the prefix of the .pos and .cel files of which you want to calculate the trajectory, and the total number of frames that you want to include in the movie. If you want to visualize the run described above, then nframes = nstep/iprint = 1000/10 = 100 (NB: positions and cell are printed in files .pos and .cel every iprint).
      Now you can execute cp2xsf.x 
              ./cp2xsf.x  < cp2xsf.in
      This will produce an output file named "si.axsf", where the suffix "axsf" stands for "animated" xsf. Now launch xcrysden... File... Open Structure... si.axsf, and enjoy the movie!


    Generating a pair-correlation function

      Generating a pair-correlation function requires a post-processing tool that you need to download. Download the post-processing source code gofr.f and and the sample input file gofr.in. Compile gofr.f 
              gfortran gofr.f -o gofr.x
      Take a look at "gofr.in". It contains information about the names of the .pos and .cel files of the trajectory for which you want to calculate the pair correlation function, and the total number of frames that you want to include in the statistics.
      Now you can execute gofr.x 
              ./gofr.x  < gofr.in
      This will produce an output file named "gofr.out", which contains three columns. Column 1 is the radius, column 2 contains the pair-correlation function, and column 3 contains the integral of the pair correlation function, normalized in such a way that it corresponds to the total number of atoms found within that radius.


    Setting a target temperature with a thermostat

      In order to bring your MD trajectory to a given target temperature, you may wish to use a thermostat, for example a Nose-Hoover thermostat. The thermostat can be activated by adding the following input variables in the "&ions" namelist, 
              ion_temperature='nose'
        tempw=4000
        fnosep=3.0
      where "tempw" is the target temperature (in Kelvin) and "fnosep" is the characteristic frequency of the thermostat (in THz!). Such frequency should be tuned within the range of characteristic frequencies of the ionic dynamics of your system. Vibrational modes with frequencies closer to the thermostat frequency will thermalize faster, but the final thermodynamical state of your system should not depend on the thermostat frequency, as long as the whole system is allowed to thermalize. Notice that dynamical properties such as vibrational frequencies, diffusion coefficients, transport properties, are affected by the thermostat, so if you're interested in those, you must switch the thermostat off after equilibration is achieved, and only then start collecting dynamical averages.



    S. Scandolo -- Latin American School on Computational Materials Science, Santiago, Chile, January 19-30, 2009