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

Open the sample input file "si.in.00". We will focus on the parameters that are specific to molecular dynamics. We will discuss the other parameters tomorrow.

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:

        /usr/local/apps/espresso-5.1/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:

        /usr/local/apps/espresso-5.1/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:

        /usr/local/apps/espresso-5.1/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 the post-processing tool cp2xsf.f90 Compile cp2xsf.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 the post-processing tool gofr.f 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 – MASTANI 2014  (modified by Vardha and Madhura)