Qbox Computing Bulk Modulus of Au

Computing Bulk Modulus of Au

This part of the tutorial is dedicated to a non trivial calculation of Bulk modulus with Qbox. We performed this calculation on a dual-core Ubuntu installation (see tutorial).

We will compute the bulk modulus by computing stress at two different volumes and taking finite differences. The calculation will be done with one atom per cell and many k-points. Qbox does not check for convergence so it is our responsability to verify for this, including energy cutoff, iteration convergence, k-point grid, system size, and pseudopotential quality. Other ways of calculating bulk modulus and convergence tests are left as exercises.

For computing the bulk modulus of gold we will need at least three files (1) a pseudopotential for Gold, (2) input files for defining cell parameters and other calculation specifications (3) a list of k-points for the FCC cell. We will put all of them in one empty directory:

 mkdir au_bulk
 cd au_bulk

We will omit the generation of the pseudopotential, we will only mention that the pseudopotential has to be consistent with the functional choosen (in this case GGA-PBE). The file is directly given here: au_bulk/au-pbe.TM.xml

 wget http://micro.stanford.edu/mediawiki-1.11.0/images/Au_bulk_au-pbe.TM.xml.txt -O au-pbe.TM.xml

Qbox was designed to be an interactive code, so each input line is actually a command given to the Qbox engine to set internal variables or perform a wavefunction optimization. The input file we will use is this au_bulk/au_stress.i . This input file calls another input file with the k-point sampling set au_bulk/kp8881fcc.i.

 wget http://micro.stanford.edu/mediawiki-1.11.0/images/Au_bulk_au_stress.i.txt -O au_stress.i
 wget http://micro.stanford.edu/mediawiki-1.11.0/images/Au_bulk_kp8881fcc.i.txt -O kp8881fcc.i

What follows is a explanation line by line of the input file au_bulk/au_stress.i. Pound sign '#' indicate comment lines.

(begin of au_stress.i)

 ## this input file compute the stress tensor for FCC Au
 # define cell
 set cell  3.9 3.9 0.000  0.000 3.9 3.9  3.9 0.000 3.9

 # define reference cell (always larger than cell),
 #   meant to smooth out discrete changes on pw basis on volume change
 set ref_cell 4 4 0  0 4 4  4 0 4

 # declare a pseudopotential file and give name to such species
 species gold au-pbe.TM.xml

 # add one atom to the cell
 atom Au gold 0 0 0

 # set number of empty states, needed for metals, nempty is 0 by default. 
 set nempty 8
 
 # energy cut off, always check convergece against this variable.
 set ecut 90

 # define kpoints using file kp8881fcc.i
 #  this will call an external file that should be in the current directory
 #  if file is not found Qbox will ignore the line and continue!
 kp8881fcc.i
 
 # activates the calculation of stress tensor.
 set stress ON

 # effective cutoff for stress calculation, improves smoothness relative to change in pw basis set
 set ecuts 85

 # set functional, default is LDA
 set xc PBE

 #set optimization algorithm in non-selfconsistent loop
 set wf_dyn PSD

 # energy cutoff for the wf optimization algorith
 # must be lower than actual ecut,
 # lower values can accelerate converge but increase instabilities  
 set ecutprec 10

 # real calculations start here
 # randomize initial wavefunction 
 randomize_wf

 # updating of charge density self consistent loop, can improve convergence
 set charge_mix_coeff 0.2

 # set smearing of fermi surface, can improve convergence for metals
 set fermi_temp 300

 # run 10 series of 5 self consistent loops with 5 non-self consistent loop 
 run 10 5 5

 # save atom position and wavefunction for later use
 save gs.xml

(end of au_stress.i)

Now we can run the calculation or submit it to a queue.

 mpirun -np 2 ~/soft/qbox-1.45.0/src/qbox au_stress.i | tee au_stress.o

'au_bulk/au_stress.o' will contain all the output from qbox, exploring this file is the only way to check for convergence. The simplest check is to see (and plot) how the total energy converged:

 grep etotal au_stress.o

Other checks can include stress tensor or electron occupancy. Once we checked converge we can look at the final output of the stress tensor. For example in our case we will see

 $ grep '<sigma_xx' au_stress_a3.9.o -A 5 | tail -6
    <sigma_xx>   4.85104411 </sigma_xx>
    <sigma_yy>   4.85104541 </sigma_yy>
    <sigma_zz>   4.85104410 </sigma_zz>
    <sigma_xy>   0.00000018 </sigma_xy>
    <sigma_yz>   0.00000010 </sigma_yz>
    <sigma_xz>  -0.00000064 </sigma_xz>

That means that the resulting stress is diagonal (as it should be) and the magnitude is 4.85104~GPa. The we can rerun an input file with a slightly different cell:

 set cell  3.855 3.855 0.000  0.000 3.855 3.855  3.855 0.000 3.855

Then read the final stress in the same way. In this way we have two volumes 114.58 and 118.64 and obtained two corresponding pressures 11.441950 and 4.85104411 Gpa. The volumes are given in atomic units but we do not need the volume units for this. A numerical approximation of the bulf modulus is simply

 (11.441950-4.85104411)/(114.58 - 118.64)* ((114.58 + 118.64)/2) 

which results numerically in 189.30 in units of GPa.

Possible improvements not shown here for simplicity:

(1) It is possible to write a script (bash or otherwise) to automatize the process. (2) It is possible to reuse the wavefunction of the first calculation in the second by loading the gs.xml wavefunction before resetting the cell size (and removing the 'randomize_wf' line)

 load gs.xml
 set cell 3.855 3.855 0 3.855 0 3.855 0 0 3.855 
 # randomize_wf

(3) In the same way it is possible to save wavefunctions in case the calculation is interrupted

 run 3 5 5
 save gs_checkpoint.xml
 run 3 5 5
 save gs_checkpoint.xml
 run 3 5 5

Further (home)work to convert this tutorial in a serious result:

• Check that the result is converged by increasing the value of ecut.
• Try duplicating the number of iteration 'run 20 5 5' to improve convergence.
• Take finite differences with smaller volume differences and near zero pressure.
• Calculate the bulk modulus using total energies at three different volumes instead of pressures or by fitting to more point near equilibrium volume.
• Change the k-point grid, k-point grids can be generated or just converted from other codes k-point grid outputs. For example from VASP, the only difference is that the k-point weights (see last column in 'au_bulk/kpt8881fcc.i') must be normalized to one.
• Use a different pseudopotential, for example this one (see details inside the file) au_bulk/Au_PBE_11e.xml