INSTRUCTIONS FOR EXERCISE
Note that the following color code has been used in this instruction sheet:
Broad headings are in red.
File names are in magenta.
Phrases to be typed into the command line are in blue.
Input parameters are in dark green.
A very brief reminder of useful linux commands:
cp A B copies file A to file B.
mv A B moves file A to file B.
mkdir C creates a directory named C.
ls C lists the names of the files and sub-directories contained in directory C
vi A or gedit A are some of the many ways in which you can open the file A for reading and/or editing.
Adding an ampersand & at the end of your command will make it run in the background, so that you can continue to issue other commands while the original one is still running.
In this exercise, you will first perform simple scf (self-consistent field) calculations on silicon, then plot isosurfaces and contours of the charge density in XCrysDen.
Create the directory exercise3/ and change to it.
Download the following file from Sandro's page:
Si.sample.in this is a sample input file, for a primitive Si cell containing 2 atoms (the same as for exercise1)
Also on the website: Si.rho.in is available as a reference, but first try to write this yourself!
Run the scf calculation:
pw.x < Si.sample.in > Si.sample.out
Extract the charge density and convert it to XCrysDen format
Note that your Si.sample.out can be read in XCrysDen, if you want to see an updated structure (for example, after a relaxation).
To plot the charge density, you will use the utility pp.x
Create an input file, Si.rho.in, using gedit. To add the correct parameters, refer to the Quantum Espresso Wiki (section pp.x):
For convenience, here is a list of keywords that are relevant for this exercise. In &INPUTPP: prefix, filplot, plot_num; and in &PLOT: iflag, output_format, fileout
When the file is created, run the post-processing tool:
pp.x < Si.rho.in
Plot the charge density in XCrysDen:
Start XCrysDen: xcrysden &
Open the file. Note that this is XSF format, not PWSCF format. To open: File > Open Structure > Open XSF and select the file.
Note that the file contains data on a 3D grid (the charge density) so that it is correspondingly larger.
Press c to draw the unit cell.
Examine the charge-density isosurface.
To explore the isosurface, select Tools > Data Grid and press OK.
Put a value in the isovalue box. Be sure to choose a value between minimum grid value and maximum grid value. Press submit. Experiment with different values.
Examine contour sections.
Remove the isosurface by unchecking display isosurface.
Change to Plane #1 tab.
Select an interesting Color Basis.
Select display color-pane.
Press submit.
You can animate the plane using the controls at the bottom.
Experiment with different colors and different plane orientations. You can even add multiple planes by using all three tabs together.
In this part of the exercise, you will calculate the frequency of the Raman-active (zone center) optical phonon of Si using the frozen-phonon method.
The phonon we will analyse is the bond-stretch phonon (Raman).
Run a series of calculations for different bond lengths. Keep the unit cell fixed. You can also keep one of the atoms at the origin, so that just the second can be changed.
Try running a series of calculations with the second atom at (0.22, 0.22, 0.22) to (0.30, 0.30, 0.30) in steps of 0.01.
Plot the energy as a function of the bond-length in Bohr (how do you find the bond length?)
Look at the graph: are we remaining in the linear-response / harmonic regime?
Using the non-linear curve fitting function of XMGrace, fit a quadratic + cubic function (y=A0*(x-A1)^2 + A2*(x-A1)^3 + A3), where 'x' is the bond length in atomic units (bohr), 'y' is the energy in Ry, and A0..A3 are fit parameters). Is the fit good?
The 2*A0 coefficient is the “spring constant” in Ry/bohr^2. Convert this first to rad.s^-1, and then to the conventional cm^-1 units by dividing by the speed of light in cm.s^-1.
How does this compare to the experimental Raman mode (~520 cm^-1)?
How does this mode frequency depend on whether you use the LDA (10.20 Bohr) or experimental (10.26 Bohr) lattice parameters?
Optional Homework: The 2*A0 parameter can also be found using the linear change in force as a function of the bond length. You can print the forces to the output file using "tprnfor=.TRUE." in &CONTROL. Remember to obtain the magnitude of the force directed along the bond -- you will receive the Cartesian components in the output file. The forces are in Ry/Bohr.