===========================================================================================
DAY-6 : The Car-Parrinello method
===========================================================================================

0) COMPILE THE CODES

a) Load the intel module: 
> module load intel/11.1.064

b) Go into the INPUT folder and compile the si_phonon_init and si_random_init programs:
> makeit_phonon
> makeit_random

c) Go back to the day-6 folder and compile the cpmd code:
> makeit 

*******************************************************************************************

1) PHONON_INIT: initial displacement along the optical mode at Gamma

Using the si_phonon_init.in input file, converge to ground state for an initial
configuration in which the 8 atoms of the simulation cell are displaced along the 111
direction:

> si_phonon_init < si_phonon_init.in > si_phonon_init.out &

This calculation produces the unformatted file fort.11, which is read at the next step by cpmd. 
The size of the atomic displacement is controlled by the variable "delta" in the 
file day-6/INPUT/src/si_phonon_init.f90-save   

********************************************************************************************

2) CPMD: Car-Parrinello MD

a) Using the input file cpmd.in run a short Car-Parrinello MD (1000 steps = 1000*5.0*4.84*10^-17 s = 0.24 ps):

> cpmd < cpmd.in > cpmd.out &

b) Grep the "mm" flag from the output file and check the behavior of the Kinetic Energy, Potential Energy and Constant of Motion.
These are the 4th, 5th and 6th fields, respectively:

> grep mm cpmd.out > ene.out

Visualize the behavior of these quantities with (for example) gnuplot:

> gnuplot
> p "ene.out" u 2:4 w l, "ene.out" u 2:5 w l, "ene.out" u 2:6 w l

Q: is the constant of motion constant?
Q: What the frequency of oscillation of the Kinetic (or Potential) Energy?

Remember that we are using a time step of 5.0 Ry a.u. (1 Ry a.u. = 4.84*10^-17 s).
The experimental frequency for this mode is 517 cm^-1 = 15.1 THz. In the original Car-Parrinello paper (PRL 1985) the frequency
they obtained is 20 THz, using the same pseudopotential and XC functional we are using.

****************************************************************************************************

3) RANDOM_INIT: random initial displacement

Using the si_random_init.in input file, converge to ground state for an initial
configuration in which the 8 atoms of the simulation cell are randomly displaced from the equilibrium positions: 

> si_random_init < si_random_init.in > si_random_init.out &

This calculation produces the unformatted file fort.11, which is read at the next step by cpmd.
The size of the atomic displacement is controlled by the variable "delta" in the 
file day-6/INPUT/src/si_random_init.f90-save

********************************************************************************************

4) CPMD: Car-Parrinello MD

a) Using the input file cpmd.in run a short Car-Parrinello MD (1000 steps = 1000*5.0*4.84*10^-17 s = 0.24 ps):

> cpmd < cpmd.in > cpmd.out &

b) Grep the "mm" flag from the output file to check the behavior of the Kinetic Energy, Potential Energy and Constant of Motion.
These are the 4th, 5th and 6th fields, respectively:

> grep mm cpmd.out > ene.out

Visualize the behavior of these quantities with (for example) gnuplot:

> gnuplot
> p "ene.out" u 2:4 w l, "ene.out" u 2:5 w l, "ene.out" u 2:6 w l

Q: Are the oscillations harmonic?

****************************************************************************************************

5) Visualization of the atomic motion with VMD

a) Modify the day-6/src/posver.f90 file to print out the coordinates in xyz format:
around line 237 print to file the number of particles and an empty line (xyz format ...):

  write(91,*) npart
  write(91,*)

around line 264, after "tau" has been updated, print out the coordinates in Angstrom:

 write(91,'(1a4,3f20.10)') 'Si',tau(:,i)*0.529177d0 



b) Run again the CPMD simulation starting from the atomic configuration in which all atoms have been displaced
in the 111 direction along the optical mode at Gamma (i.e. repeat steps 1 and 2). 
This time you will find the file fort.91, containing the positions as a function of time.
Move the file to pos.xyz:

> mv fort.91 pos.xyz

Visualize it with VMD:
> module load vmd
> vmd pos.xyz

Change the representation from LINES to CPK:
Graphics > Representations > CPK

Visualize the distance between two atoms (for example atoms "0" and "4"):

Mouse > Labels > Bonds
and click on two relevant atoms.

Graph the distance as function of time:

Graphics > Labels > Bonds > Graph

You can export the data file in text format or directly in xmgrace for data analysis.

Q: What is the frequency of oscillation of this distance? 
(You can estimate it directly from the graph or you can fit the curve with a sine function in xmgrace)



c)  Run again the CPMD simulation starting from the atomic configuration in which all atoms have been displaced
randomly (i.e. repeat steps 3 and 4).

Move the file to pos.xyz:

> mv fort.91 pos.xyz

Visualize it with VMD:
> vmd pos.xyz 

********************************************************************************************************************

6) CPMD: Car-Parrinello MD with different initial displacements

a) Modify the variable delta that controls the magnitude of the initial displacement from the equilibrium position
in the file day-6/INPUT/src/si_phonon_init.f90-save. For example change it from -0.002 to 0.005.

Minimize the total energy, run the CPMD simulation, recompute the vibrational frequency by plotting the distance between
atoms "0" and "4".

Q: Is the frequency still the same or not?


