Physical Properties of Boron Nitride

by Wahib Aggoune & Pasquale Pavone for exciting nitrogen

Purpose: In this tutorial, you will learn how to set up and execute exciting calculations for three phases of boron nitride, i.e., the wurtzite (w-BN), hexagonal (h-BN), and cubic (c-BN) phases. Examples of electronic band-structure calculations and full structure optimization will be also shown. Finally, we show how to investigate the optical excitations in h-BN, as an example for the calculation of excited-states properties.


0. Preliminary notes

Read the following paragraphs before starting with the rest of this tutorial!

Here is a list of the scripts which are relevant for this tutorial with a short description.

  • SETUP-excitingroot.sh: Python script for replacing the environment variable in the input.xml.
  • EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
  • OPTIMIZE-lattice.py: Python script. A manager program which calls the setup and analyze scripts
  • OPTIMIZE-submit.sh: (Bash) shell script for running a series of exciting calculations.
  • LATTICE-parameters.sh: (Bash) shell script for extracting the values of the lattice parameters out of any exciting input file.
  • SETUP-volume-optimization.py: Python script for generating structures at different volumes.
  • EXECUTE-volume-optimization.sh: (Bash) shell script for running a series of exciting calculations.
  • PLOT-birch.py: Python script for fitting energy-vs-strain curves using the Birch-Murnaghan equation of state.
  • PLOT-newbirch.py: Python script for fitting energy-vs-volume curves using the Birch-Murnaghan equation of state in polynomial form.
  • PLOT-exciton-weights.py: Python visualization scripts for plotting the excitonic weights along a band structure path.

Furthermore, this tutorial requires the usage of the following visualization packages:

Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.

  • From now on the symbol $ will indicate the shell prompt.

1. Three different phases of boron nitride

As shown in the following figure, boron nitride exhibits three different phases: The hexagonal (honeycomb layers similar to graphite), wurtzite, and cubic (similar to diamond) ones.

str-BN.png

Similar to graphite, bulk h-BN has a layered structure formed by alternating h-BN monolayers, in the AA' stacking arrangement, which are held together by weak van der Waals forces. The unit cell contains four atoms distributed in two different layers, each one can be seen as a triangular lattice with a basis of two atoms (N and B) per unit cell. The lattice vectors in this case can be written as

(1)
\begin{align} {\bf a}_1= a \left(\frac{1}{2}, \frac{\sqrt{3}}{2}, 0\right) \qquad {\rm and} \qquad {\bf a}_2= a \left(-\frac{1}{2}, \frac{\sqrt{3}}{2}, 0\right) \qquad {\rm and} \qquad {\bf a}_3= a \left(0, 0, \frac{c}{a}\right)\;, \end{align}

where $a$ = 2.55 Å and $c/a$ = 2.66 are the experimental lattice constants. The crystal structure of the wurtzite phase of BN is similar to that of h-BN. However, the N and B atoms of different layers in the wurzite structure are bound by strong covalent bonds. The lattice vectors have the same form as for h-BN, with different values of the lattice parameters, being $a$ = 2.55 Å and $c/a$ = 1.64.

For the cubic phase (zincblende structure) the fcc lattice parameter is $a$ = 3.61 Å. The lattice vectors are defined as follow

(2)
\begin{align} {\bf a}_1= a \left(\frac{1}{2}, \frac{1}{2}, 0 \right) \qquad {\rm and} \qquad {\bf a}_2= a \left(0, \frac{1}{2}, \frac{1}{2}\right) \qquad {\rm and} \qquad {\bf a}_3= a \left(\frac{1}{2}, 0, \frac{1}{2}\right)\;, \end{align}

2. Kohn-Sham electronic properties of boron nitride


2. 1 Wurtzite boron nitride (w-BN)

i) Preparation of the input file

The first step before starting an exciting calculations is to create a directory for the example that you want to investigate. Thus, we create a directory $IMPRSWORK/w-BN and we move inside it.

$ mkdir $IMPRSWORK/w-BN
$ cd $IMPRSWORK/w-BN

The first example of an exciting input file (input.xml) for w-BN is given in the following.

<input>
 
   <title> w-BN </title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="4.8188">
         <basevect>  0.500000000  0.86602540  0.000000000 </basevect>
         <basevect> -0.500000000  0.86602540  0.000000000 </basevect>
         <basevect>  0.000000000  0.00000000  1.635000000 </basevect>
      </crystal>
      <species speciesfile="B.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.0000000000"/>
         <atom coord="0.333333333  0.333333333  0.5000000000"/>
      </species>
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.3748987557"/>
         <atom coord="0.333333333  0.333333333  0.8748987557"/>
      </species>
 
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE"
      ngridk="6 6 4"
      rgkmax="5">
   </groundstate>
 
</input>

Input files for exciting are written in the XML format and are typically named input.xml. The XML format allows your data to be written in a structured way. Figuratively speaking, an exciting input is pretty much like an article with its sections and subsections. In case of XML data, sections and subsections are called elements.

Let us examine this example bit-by-bit. The first thing to be said is that an XML is not sensitive to line indentation. However, for the sake of clarity, line indentation is used in all examples of this tutorial to illustrate the hierarchy among elements. Now, let's come back to the input.xml shown above. The first and the last line indicate the beginning and the end of the input.

The element title contains some freely chosen text simply to describe the calculation. Keywords <title> and </title> indicate the beginning and the end of the element (note in the second keyword the presence and position of the slash /).

The next element, structure, describes the geometry and the chemical composition of a studied system. Notice that the declaration of the structure section contains as additional information the parameter speciespath.

    <structure speciespath="$EXCITINGROOT/species/">

Such parameters in the XML language are called attributes, and their values are always given in quotes regardless whether it is numerical, symbolic, or boolean information. In particular, the attribute speciespath defines the location, where the files with the data about chemical elements are stored.

Please, note that exciting is always using atomic units (Hartree, Bohr, etc.)! Thus, you can notice that in the example above the value of the experimental lattice constant has been considered (in units of bohr!).

...
    <crystal scale="4.8188">
...

Let us now investigate in more detail the other instructions given in the input file. The lattice vectors of w-BN are chosen according to Eq. (1):

...
      <basevect>  0.500000000  0.86602540  0.000000000 </basevect>
      <basevect> -0.500000000  0.86602540  0.000000000 </basevect>
      <basevect>  0.000000000  0.00000000  1.635000000 </basevect>
...

notice that the $c$ parameter in w-BN is 4.17 Å, which gives 1.635 for the $c/a$ parameter.

The position of the four atoms (2B and 2N) inside the unit cell of w-BN is chosen accordingly to the figure shown above.

...
      <species speciesfile="B.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.0000000000"/>
         <atom coord="0.333333333  0.333333333  0.5000000000"/>
      </species>
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.3748987557"/>
         <atom coord="0.333333333  0.333333333  0.8748987557"/>
      </species>
...

The last part of the input file contains a list of the parameters which are set for this calculation:

...
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE"
      ngridk="6 6 4"
      rgkmax="6">
   </groundstate>
...

In particular, according to this part:

  • The GGA_PBE functional is chosen for the exchange-correlation energy.
  • Brillouin-zone integration are performed by summing contributions from a 6$\times$6$\times$4 k-points grid. Notice that the number of k-points along the out-of-plane direction is set to 4 which is approximatly ${6}(c/a)^{-1}$.
  • The number of APW basis functions used in the calculation is set by assigning a value to the attribute rgkmax.
ii) Execute calculations

To perform the first example calculation, copy the complete input listed above inside a file input.xml which you should create in the current directory.

N.B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT using the command

$ SETUP-excitingroot.sh

You can now run the script EXECUTE-single.sh in the directory where the input file input.xml has been created. If you want to save the results of this run in the subdirectory first-example you should run as in the following.

$ EXECUTE-single.sh first-example

===> Output directory is "first-example" <===

Running exciting for file input.xml -------------------------------------

   Elapsed time = 0m14s

Run completed for file input.xml ----------------------------------------

$

This (preliminary) calculation solve self-consistently the Kohn-Sham equations and allows to determine the ground-state total energy of w-BN in the configuration specified in the input file. More information can be found inside the file first-example/INFO.OUT.

iii) Calculate density of states

After you have completed the ground-state run and have obtained the corresponding total energy, you can go for more properties of the system. One of the most fundamental ones is the density of states (DOS). The DOS gives you information on the energy levels in your system, or — more precisely — about how many electronic states there are at any given energy. To calculate it, move inside the directory first-example.

$ cd first-example

Then, you need to do the following simple modifications in the file input.xml inside this directory (for more details, see Input Reference):

  1. modify the value of the attribute do to "skip" to the element groundstate;
  2. add the element properties after the groundstate element;
  3. insert the subelement dos into the element properties;
  4. add the attribute nsmdos = "3" to the element dos.

The corresponding part of the input.xml should now look like this:

...
   <groundstate
      do="skip"
      ...
   </groundstate>
 
   <properties>
      <dos nsmdos="3"/>
   </properties>
...

All other instructions inside input.xml remain unchanged.

Now, execute exciting again by typing on the command line:

$ time excitingser

The run should take only a few seconds. Then, to visualize the DOS, execute

$ DOS w-BN_dos.agr

This produces the file w-BN_dos.agr which can be visualized using the program xmgrace. Open it with the command

$ xmgrace w-BN_dos.agr >/dev/null 2>&1 &
iv) Calculate electronic band structure

Now, we are ready for a more detailed view on the Kohn-Sham electronic structure: The band structure. In addition to the energy of each state, the band structure shows the dependence of the energy eigenvalues on the coordinates in k-space.

To calculate the band structure of w-BN, insert the subelement bandstructure in the element properties with the following specifications:

...
   <properties>
 
      <bandstructure>
         <plot1d>
            <path steps="100">
               <point coord="0.33333333  -0.33333333  0.0" label="K"    />
               <point coord="0.0          0.0         0.0" label="GAMMA"/>
               <point coord="0.5          0.0         0.0" label="M"    />
               <point coord="0.33333333  -0.33333333  0.0" label="K"    />
               <point coord="0.33333333  -0.33333333  0.5" label="H"    />
               <point coord="0.0          0.0         0.5" label="A"    /> 
               <point coord="0.5          0.0         0.5" label="L"    />
               <point coord="0.33333333  -0.33333333  0.5" label="H"    />
            </path>
         </plot1d>
      </bandstructure>
 
   </properties>
...

Now execute excitingser again on the command line:

$ time excitingser

This makes the code produce the band structure, which is written to bandstructure.xml. To visualize the band structure, execute the following command:

$ BAND

Now, you have produced the xmgrace file w-BN_bandstructure.agr, which you can open, visualize and manipulate with xmgrace:

$ xmgrace w-BN_bandstructure.agr >/dev/null 2>&1 &

In order to visualize the character of these bands (i.e., to identify which atom and kind of orbital generate a given band), we add to the subelement bandstructure in the element properties the following attribute:

...
      <bandstructure character="true">
...

Now execute excitingser again on the command line:

$ time excitingser

This will takes few seconds and the code produce the band structure with character, which is written to BAND_Sss_Aaaaa.OUT for all atoms. To visualize these results, execute the command:

$ BAND

The command generates several files of the type w-BN_band_atom_n_orbital.agr related to each atoms (e.g., boron_1, nitrogen_2, etc.) and orbitals (e.g., l0 for l=0, l1 for l=1, etc.). These files can be visualized using the program xmgrace. For instance, the contribution of the p orbitals of the first N atom can be visualized by typing

$ xmgrace w-BN_band_nitrogen_1_l1.agr >/dev/null 2>&1 &

To improve the quality of the visualization, you can follow the following steps. In the menu bar of xmgrace on the top, press Plot -> Set appearance…. Then, click inside the Select set: window and type Ctrl A to select all the bands. In the Symbol properties window set the Type to Circle. By selecting Accept, black circles appear, indicating the different bands and regions in the Brillouin zone that are dominated by p orbital of first N atom in the input. You can do the same for other orbital and atoms to see their contribution. In the case of w-BN the p orbitals of the N and B atoms dominate the valence band maximum and conduction band minimum, respectively.

Notice that the different files contain information about contributions of each atom and orbital with respect to the same electronic band structure.


2. 2 Cubic boron nitride (c-BN)

In order to study the cubic structure c-BN, we create now a directory $IMPRSWORK/c-BN and proceed in a similar way to the previous phase. The new input file, input.xml, for c-BN is given in the following. Please, notice that executing this input file will perform in a "single shot" the total energy, density of states, and electronic band-structure calculations.

<input>
 
   <title> c-BN </title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="6.84">
         <basevect> 0.5  0.5  0.0 </basevect>
         <basevect> 0.5  0.0  0.5 </basevect>
         <basevect> 0.0  0.5  0.5 </basevect>
      </crystal>
      <species speciesfile="B.xml" rmt="1.25">
         <atom coord="0.00000000  0.00000000  0.00000000"/>
      </species>
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.25000000  0.25000000  0.25000000"/>
      </species>
 
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE"
      ngridk="4 4 4"
      rgkmax="6">
   </groundstate>
 
   <properties>
 
      <dos nsmdos="3"/>
 
      <bandstructure character="true">
         <plot1d>
            <path steps="100">
               <point coord="0.50  0.00  0.50" label="X"    />
               <point coord="0.50  0.25  0.75" label="W"    />
               <point coord="0.50  0.50  0.50" label="L"    />
               <point coord="0.00  0.00  0.00" label="GAMMA"/>
               <point coord="0.50  0.00  0.50" label="X"    />
            </path>
         </plot1d>
      </bandstructure>
 
   </properties>
 
</input>

Following all instructions you got for the previous phase, you should be able tu running successfully exciting as in the following example.

$ EXECUTE-single.sh first-example

===> Output directory is "first-example" <===

Running exciting for file input.xml -------------------------------------

...

   Elapsed time = 0m7s

Run completed for file input.xml ----------------------------------------

$

You can now plot the density of states and the Kohn-Sham electronic band structure as explained above.


2. 3 Hexagonal boron nitride (h-BN)

In order to study the third structure h-BN, we create now a directory $IMPRSWORK/h-BN and proceed in a similar way to the previous phases. The "single shot" input file, input.xml, for h-BN is given in the following.

<input>
 
   <title> Bulk_h-BN  </title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="4.73">
         <basevect>  0.5  0.866025404  0.00 </basevect>
         <basevect> -0.5  0.866025404  0.00 </basevect>
         <basevect>  0.0  0.000000000  2.66 </basevect>
      </crystal>
      <species speciesfile="B.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.0"/>
         <atom coord="0.333333333  0.333333333  0.5"/>
      </species>
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.333333333  0.333333333  0.0"/>
         <atom coord="0.000000000  0.000000000  0.5"/>
      </species>
 
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE"
      ngridk="6 6 2"
      rgkmax="6">
   </groundstate>
 
   <properties>
 
      <dos nsmdos="3"/>
 
      <bandstructure character="true">
         <plot1d>
            <path steps="100">
               <point coord="0.33333333  -0.33333333  0.0" label="K"    />
               <point coord="0.0          0.0         0.0" label="GAMMA"/>
               <point coord="0.5          0.0         0.0" label="M"    />
               <point coord="0.33333333  -0.33333333  0.0" label="K"    />
               <point coord="0.33333333  -0.33333333  0.5" label="H"    />
               <point coord="0.0          0.0         0.5" label="A"    /> 
               <point coord="0.5          0.0         0.5" label="L"    />
               <point coord="0.33333333  -0.33333333  0.5" label="H"    />
            </path>
         </plot1d>
      </bandstructure>
 
   </properties>
 
</input>

After performing the calculation for this phase, you can try to compare the results for the different phases.

Question: What is the character of the valence-band minimum and conduction-band maximum in the three phases of boron nitride? Why?


3. Full structure optimization

Up to here, all calculations have been performed considering BN phases in their experimental equilibrium configurations. However, exciting is also able to predict the equilibrium configuration, i.e., the equilibrium lattice constant of the different phases. In this section, it will be shown how to proceed in order to find the equilibrium geometry for each system.

Please note:

  • The values of the parameter ngridk and rgkmax in the input file have been chosen in order to allow for a fast calculation. Of course, the complete optimization procedure presented in this tutorial should be repeated for more accurate values of all computational parameters up to the desired accuracy in the resulting lattice constants.

3. 1 w-BN

We take as a starting point the input file used in the last section. We create a new working directory $IMPRSWORK/h-BN/structure-optimization and we move inside it.

$ mkdir $IMPRSWORK/w-BN/structure-optimization
$ cd $IMPRSWORK/w-BN/structure-optimization

Copy inside this directory the input file you used in the first example above.

$ cp $IMPRSWORK/w-BN/input.xml ./w-BN_opt.xml

Before performing the structure optimization, some attention should be paid at the input file: The do attribute of the groundstate element should be restored to "fromscratch". Furthermore, the possibly present properties block could be omitted. Finally, a new element, relax, must be added just after the element groundstate. See Input Reference for the meaning of the different attributes of relax. The presence of the element relax guarantees the optimization of the relative position of the atoms in the unit cell at any value of the general lattice parameters of the crystal.

...
   <relax method="harmonic" taunewton="1.0"/>
 
</input>

To optimize the general crystal structure of w-BN, you will set up and execute a series of calculations for different volumes (at constant c/a ratio) and for differentc/a ratios (at constant volume) as shown in General lattice optimization. The final goal is to find the set of lattice constants which minimize the total energy of the system. The values of the initial lattice parameters, $a, b, c, \alpha, \beta,$ and $\gamma$, corresponding to the initial input file w-BN_opt.xml can be extracted using the following command

$ LATTICE-parameters.sh w-BN_opt.xml

As an output on the screen you should get

     Bravais lattice: Hexagonal
     a          b          c          alpha      beta       gamma
     4.81880    4.81880    7.87874   90.00000   90.00000  120.00000

STEP1) Optimizing the volume

At the first step, we optimize the energy with respect to the volume. In order to generate input files for a series of volumes you have to use the script OPTIMIZE-lattice.py.

$ OPTIMIZE-lattice.py 

>>>> Please enter the exciting input file name: w-BN_opt.xml

     Number and name of space group: 186 (P 63 m c)     
     Hexagonal I structure in the Laue classification.

     Which parameter would you like to optimize?
     1 ... volume                               
     2 ... c/a ratio with constant volume       
>>>> Please choose '1' or '2': 1

>>>> Please enter the maximum physical strain value  
     The suggested value is between 0.001 and 0.050: 0.03
     The maximum physical strain is 0.03

>>>> Please enter the number of the distorted structures [odd number > 4]: 5
     The number of the distorted structures is 5

$

Entry values must be typed on the screen when requested. In this case, entries are the following.

  1. w-BN_opt.xml, the name of the input file you have created.
  2. 1, the type of optimization you choose.
  3. 0.03, the absolute value of the maximum strain, ϵmax, for which we want to perform the calculation. Notice that in this case the physical strain ϵ is defined by the relationship V=Vinitial·(1+ϵ)3, where Vinitial is the volume of the unit cell defined in w-BN_opt.xml.
  4. 5, the number of deformed structures equally spaced in strain, which are generated between the maximum negative strain and the maximum positive one.

Now, you must move to the directory VOL and execute the calculation there.

$ cd VOL

To execute the calculations, you have to run the script OPTIMIZE-submit.sh. If you do so, the screen output will be similar to the following.

$ OPTIMIZE-submit.sh 

SCF calculation of "vol_01" starts --------------------------------
   Elapsed time = 0m25s

SCF calculation of "vol_02" starts --------------------------------
   Elapsed time = 0m25s

...

SCF calculation of "vol_05" starts --------------------------------
   Elapsed time = 0m27s

$

When calculations finished to run, results can be analyzed. In order to do this, you have to run again the OPTIMIZE-lattice.py python script in the current directory.

$ OPTIMIZE-lattice.py

>>>> Murnaghan or Birch-Murnaghan EOS: [M/B]

At this point, the script is asking whether you desire to use a Murnaghan (M) or Birch-Murnaghan (B) equation of state for extracting equilibrium parameters such as the equilibrium volume and bulk modulus.

Therefore, if you type B or b, you will see as a result that the script has generated some new files. Among them, e.g., the BM_eos.out file, for instance, contains all final fit parameters and pressures which are calculated from equation of state fitting procedure. Also, the script generates a plot, PostScript (BM_eos.eps) or PNG (BM_eos.png) file, where you can see the calculated energies and the fit curve. On this plot, you can also find the optimized values of the parameters appearing in the equation of state (minimum energy, equilibrium volume, bulk modulus, and bulk modulus pressure derivative).

A file corresponding to an exciting input file for the optimized geometry is created with the name BM-optimized.xml. The values lattice parameters $a, b, c, \alpha, \beta,$ and $\gamma$ at this step can be visualized using the command

$ LATTICE-parameters.sh BM-optimized.xml

As an output on the screen you should get

     Bravais lattice: Hexagonal
     a          b          c          alpha      beta       gamma
     4.86883    4.86883    7.96053   90.00000   90.00000  120.00000

Due to the fact that in this step the volume only has been optimized, the relative variation with respect to the initial values of the three lattice parameters $a, b,$ and $c$ is the same.

At this point, you have performed the first optimization step by varying only the volume. In order to be prepared for the next step, you should move now to the parent directory and rename the VOL directory to 1-VOL (first step, optimizing only the volume). Then, you should copy the BM-optimized.xml file to the current directory with the new name 1-VOL.xml. This file will be used as the input file in the next step.

$ cd ..
$ mv VOL 1-VOL
$ cp 1-VOL/BM-optimized.xml 1-VOL.xml

STEP2) Optimizing the $c/a$ ratio

In order to perform the next optimization step, you have to run again OPTIMIZE-lattice.py in the /home/tutorials/w-BN/structure-optimization directory, typing entries like in the following.

$ OPTIMIZE-lattice.py

>>>> Please enter the exciting input file name: 1-VOL.xml

     Number and name of space group: 186 (P 63 m c)      
     Hexagonal I structure in the Laue classification.

     Which parameter would you like to optimize?
     1 ... volume                               
     2 ... c/a ratio with constant volume       
>>>> Please choose '1' or '2': 2

>>>> Please enter the maximum physical strain value  
     The suggested value is between 0.001 and 0.050: 0.03
     The maximum physical strain is 0.03

>>>> Please enter the number of the distorted structures [odd number > 4]: 5
     The number of the distorted structures is 5

$

Notice that in this case the physical strain ϵ is defined by the relationship (c/a)=(c/a)initial·(1+ϵ), where (c/a)initial corresponds to the unit cell defined in 1-VOL.xml. Now, move to the COA directory and run OPTIMIZE-submit.sh.

$ cd COA
$ OPTIMIZE-submit.sh

When calculations are finished, run again OPTIMIZE-lattice.py.

$ OPTIMIZE-lattice.py

 =====================================================================
 Optimized lattice parameter saved into the file: "coa-optimized.xml".
 =====================================================================

$

In this case, the optimization is performed using a fourth-order polynomial fit for calculating the minimum energy and the corresponding strain. The resulting plot is also available as PostScript file as coa.eps).

At this point, you have completed the second optimization step (by varying only the $c/a$ ratio). The optimized structure is saved in the coa-optimized.xml file. The current values of the lattice parameter are obtained once more by running the script LATTICE-parameters.sh

$ LATTICE-parameters.sh coa-optimized.xml

     Bravais lattice: Hexagonal
     a          b          c          alpha      beta       gamma
     4.86314    4.86314    7.97918   90.00000   90.00000  120.00000
$

Here, you can notice that by changing the $c/a$ ratio at fixed $V$, the values of $a, b$ and $c$ differ from the ones at the previous step. Now, you should move out to the parent directory, rename COA to 2-COA (second optimization step, varying only $c/a$). Then, copy the coa-optimized.xml file to the current directory with the name 2-COA.xml.

$ cd ..
$ mv COA 2-COA
$ cp 2-COA/coa-optimized.xml 2-COA.xml

STEP3) Optimizing again the volume

Repeat now the procedure already explained in STEP1, running the script OPTIMIZE-lattice.py and using as entries values 2-COA.xml, 1, 0.03, and 5 in the given order. After having performed the calculation (running the script OPTIMIZE-submit.sh inside the directory VOL), you run OPTIMIZE-lattice.py.

At this point, you have optimized the volume for the second time. Follow the last part of STEP1 and copy the file BM-optimized.xml to the parent directory under the name 3-VOL.xml.

$ cd ..
$ mv VOL 3-VOL
$ cp 3-VOL/BM-optimized.xml 3-VOL.xml

STEP4) Optimizing again the $c/a$ ratio

Proceed in a similar way to STEP2. Run the script OPTIMIZE-lattice.py using as entries values 3-VOL.xml, 2, 0.03, and 5 in the given order. Using the same procedure as in the previous steps, check the resulting lattice parameters.

At this point, the second optimization of the $c/a$ ratio has been completed and the new optimized structure has been saved in coa-optimized.xml.


Reaching convergence

In the following table you find a summary of the results for the lattice parameters of the first 4 optimization steps.

Step OPT $a_{\hspace{0.1mm}}$min [Bohr] $b_{\hspace{0.1mm}}$min [Bohr] $c_{\hspace{0.1mm}}$min [Bohr] $\alpha_{\hspace{0.1mm}}$min [degrees] $\beta_{\hspace{0.1mm}}$min [degrees] $\gamma_{\hspace{0.1mm}}$min [degrees]
0 4.81880 4.81880 7.87874 90.00000 90.00000 120.00000
1 $V$ 4.86883 4.86883 7.96053 90.00000 90.00000 120.00000
2 $c/a$ 4.86314 4.86314 7.97918 90.00000 90.00000 120.00000
3 $V$ 4.86661 4.86661 7.98488 90.00000 90.00000 120.00000
4 $c/a$ 4.86319 4.86319 7.99610 90.00000 90.00000 120.00000

As you can see from the previous table, the iterative procedure which we used allows for a reduction of the variation of the lattice constants from one step to the next. Residual oscillations can be attributed to the non-optimal choice of the computational parameters such as ngridk, rgkmax, epsengy, and the number of calculations at each step. Improving the "quality" of these parameter will improve the quality of the results.


3. 2 c-BN

Since the lattice parameters $a, b, c$ in the cubic BN phase are equal, in this system we optimize only the volume of the unit cell as shown in Volume optimization for cubic systems. We create a new working directory $IMPRSWORK/c-BN/structure-optimization and we move inside it.

Copy inside this directory the input file you used in the first example above

$ cp $IMPRSWORK/c-BN/input.xml ./input.xml

and then remove the element properties after the groundstate element. Due to the cubic symmetry, there is here no extra internal relaxation. Therefore, the relax element can be omitted.

In order to generate input files for a series of volumes you have to run the script SETUP-volume-optimization.py. A directory name can be specified by adding the name in the command line.

$ SETUP-volume-optimization.py c-BN-vol-opt

The script SETUP-volume-optimization.py produces the following output on the screen.

$ SETUP-volume-optimization.py c-BN-vol-opt

Enter the number of volume values >>>> 9

$

The working directory is named c-BN-vol-opt, the (on screen) input entry is preceded by the symbol ">>>>". The entry value must be typed on the screen when requested. This entry (in our example 9) is the number of structures which are generated with equally spaced intervals of the lattice constant with a variation of between -5% and +5% from the reference lattice constant.

To execute the series of calculation with input files created by SETUP-volume-optimization.py you have to run the script EXECUTE-volume-optimization.sh. If a name for the working directory has been specified, then you must give it here, too.

$ EXECUTE-volume-optimization.sh c-BN-vol-opt

===> Output directory is " c-BN-vol-opt" <===

Running exciting for file input-01.xml ----------------------------------

...

Run completed for file input-09.xml ------------------------------------

$

After the complete run, the results of the calculation for the input file input-i.xml are contained in the subdirectory (of the working directory c-BN-vol-opt) rundir-i where i is running from 01 to the total number of volumes. The data for the energy-vs-volume curve are contained in the file energy-vs-volume inside c-BN-vol-opt.

In order to analyse the results of the calculations you first have to move to the working directory.

$ cd c-BN-vol-opt

At this point, you can use the python script PLOT-newbirch.py for performing the fit using the Birch-Murnaghan equation of state in order to extract the equilibrium volume and the corresponding lattice parameter.

$ PLOT-newbirch.py 

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     V0        B0         Bp        a-sc       a-bcc      a-fcc     log(chi)
   80.44644   374.444     4.114     4.3169     5.4389     6.8526      -4.56
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

$

The resulting fcc lattice parameter of c-BN is now 6.8526 (Bohr). This script generates also the output file PLOT.png.


3. 3 h-BN

In order to optimize the crystal structure of h-BN, we can follow the same procedure used for the w-BN structure shown above. In addition, since van der Waals forces play an important role in h-BN, we use an empirical energy correction, as given by Tkatchenko-Scheffler, to the standard DFT total energy by setting the attribute vdWcorrection="TSvdW" in the groundstate element. This correction will increase the execution time. In addition, the values of the mesh of k-points ngridk and the size of the basis set rgkmax should be increased to allow for a better convergence of the resulting lattice parameters.

As usual, we create a new working directory $IMPRSWORK/h-BN/structure-optimization and we move inside it. We copy inside this directory the input file you used above.

$ cp $IMPRSWORK/c-BN/input.xml ./input.xml

and then remove the element properties after the groundstate element. Other modifications to the input file are shown below.
   <groundstate
      ...
      vdWcorrection="TSvdW">
   </groundstate>
 
   <relax method="harmonic" taunewton="1.0"/>
 
</input>

Following now the same procedure used for w-BN, you can obtain the theoretical equilibrium lattice constants of h-BN. Notice that in the latter case, the
choice of 0.05 for the maximum strain in the input of OPTIMIZE-lattice.py could improve the quality of the results.


4. "Quality" of the calculations

In order to be able to rely on your calculation, you need to understand how the choice of the computational parameters change of the values of the physical quantities which are relevant for you. For most of the calculations performed in this tutorial, the most relevant computational parameters are listed below:

  • The mesh of k-points (groundstate attribute: ngridk).
  • The size of the basis set for expanding the wave function (groundstate attribute: rgkmax).
  • The maximum strain amplitude and the number of calculated structures at each optimization step.

Try to get a feeling on how the variation of these parameters changes the physical properties which are relevant in this tutorial.


5. Excited-states properties: The optical spectrum and electron-hole pairs

In this section you will learn how to perform a basic Bethe-Salpeter-equation (BSE) calculation in order to compute the optical spectrum and analyze the electron-hole pairs in the h-BN. The resulting imaginary part of the microscopic dielectric function is the quantity which is directly measurable as optical absorption. In order to analyze the resulting full BSE spectra (singlet, including the electron-hole interaction), we compute also the corresponding independent-particle spectra (IP, neglecting the electron-hole interaction). Further details can be found in Excited states from BSE.

Please note: that the values of the parameter in the input file have been chosen in order to speed up the calculations. The converged results for structural, electronic, and optical properties of h-BN can be found in the paper by Aggoune et al..

Now, we start by creating a new working directory $IMPRSWORK/h-BN/BSE and moving inside it.

In order to calculate excited-states properties, we compute first the ground-state eigenvalues and eigenvectors. The corresponding input.xml is shown in the following. Notice that we modify the lattice vectors to the optimized values, by setting $c/a=2.65$ and scale to 4.737 Bohr. In addition, we set the mesh of k-points (groundstate attribute: ngridk) to 9 9 2 in order to improve the quality of the resulting spectra.

<input>
 
   <title> Bulk_h-BN </title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="4.737">
         <basevect>  0.5  0.866025404  0.00   </basevect>
         <basevect> -0.5  0.866025404  0.00   </basevect>
         <basevect>  0.0  0.000000000  2.6578 </basevect>
      </crystal>
      <species speciesfile="B.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.0"/>
         <atom coord="0.333333333  0.333333333  0.5"/>
      </species>
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.333333333  0.333333333  0.0"/>
         <atom coord="0.000000000  0.000000000  0.5"/>
      </species>
 
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE"
      ngridk="9 9 2"
      rgkmax="5">
   </groundstate>
 
   <properties>
      <bandstructure>
         <plot1d>
            <path steps="100">
               <point coord="0.33333333  -0.33333333  0.0" label="K"    />
               <point coord="0.0          0.0         0.0" label="GAMMA"/>
               <point coord="0.5          0.0         0.0" label="M"    />
               <point coord="0.33333333  -0.33333333  0.0" label="K"    />
               <point coord="0.33333333  -0.33333333  0.5" label="H"    />
               <point coord="0.0          0.0         0.5" label="A"    /> 
               <point coord="0.5          0.0         0.5" label="L"    />
               <point coord="0.33333333  -0.33333333  0.5" label="H"    />
            </path>
         </plot1d>
      </bandstructure>
   </properties>
 
</input>

i) Execute the ground-state calculations using the optimized structure.

To perform the ground-state calculation, copy the complete input listed above inside a file input.xml which you should create in the current directory.

N.B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT using the command

$ SETUP-excitingroot.sh

You can now run the calculation, which will takes few seconds, as in the following.

$ time excitingser

This calculation solve self-consistently the Kohn-Sham equations and allows to determine the ground-state eigenvalues and eigenvectors as well as the the bandstructure.dat file where the band structure is printed. This data are needed in order to to proceed with the BSE calculations.


ii) Start the BSE calculation

You may start by copying the input file to input-BSE.xml. This file must be modified as foffows:

  1. set the value of the attribute do to "skip" inside the element groundstate;
  2. delete the full properties block;
  3. insert the element xs just after groundstate as in the next example (see Input Reference and Excited states from BSE for more information about the meaning of the different attributes):
...
   <xs 
      xstype="BSE" 
      ngridk="9 9 2"
      rgkmax="5"
      vkloff="0.001 0.02 0.003"
      ngridq="9 9 2"
      nempty="80"
      gqmax="2.0"
      broad="0.0036"
      scissor="0.0"
      tappinfo="true"
      tevout="true">
 
      <energywindow intv="0.0 1.0" points="1200"/>
 
      <screening screentype="full" nempty="100"/>
 
      <BSE bsetype="singlet" nstlbse="6 8 1 3" />
 
      <qpointset>
         <qpoint>0.0 0.0 0.0</qpoint>
      </qpointset>
 
      <storeexcitons selectenergy="false" 
         MinNumberExcitons="1" 
         MaxNumberExcitons="5"/>
 
   </xs>
...

We are now ready to run the BSE calculation:

$ time excitingser input-BSE.xml &

While the calculation is running, some strings will be printed on screen. A successfully completed BSE run will produce the following lines:

Calculating RPA Dielectric Function:                100.000 
Calculating Screened Coulomb Potential:             100.000 
Calculating Screened Coulomb Matrix Elements:       100.000 
Calculating Plane-wave matrix elements:             100.000 
Calculating Exchange Interaction Matrix Elements:   100.000 
Solving BSE Eigenvalue Problem:                     100.000 
   Elapsed time = 8m3s

A number of files and folders are now present in the working directory. Most of them contain technical information about the calculation that is not strictly related to the physical interpretation of the results. Since we are interested in the optical absorption spectrum of h-BN, given by the imaginary part of the macroscopic dielectric function, we focus on the folder EPSILON. In addition, information about the excitons is found inside the folder EXCITON. We notice that the subelement storeexcitons allows to store the information about the spatial distribution of the first 5 excitons both in real and reciprocal space.

In order to perform a complete analysis of the results, we also calculate the independent-particle spectrum where electron-hole interactions are neglected. First, you have to copy the input file input-BSE.xml to input-BSE-IP.xml. Then, the latter file must be modified by setting the attribute bsetype to "IP" instead of "singlet", and by deleting the subelement storeexcitons. We add also the subelement plan to avoid repeating other calculations. The corresponding part of the input-BSE-IP.xml should now look like this:

...
      <BSE bsetype="IP" ... />
...
      <plan>
         <doonly task="bse"/>
      </plan>
...

Run again the calculation, which this time is very fast:

$ time excitingser input-BSE-IP.xml &

All the interesting results are stored In the folder EPSILON.


iii) Analysis of the spectrum

Since we are dealing with a layered hexagonal system, we focus only on the in-plane optical component, corresponding to the label OC11.
In order to plot the full BSE spectrum for the OC11 component, execute the following commands:

$ cd EPSILON
$ cp EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT.xml tmp.xml
$ xsltproc --stringparam filename dielectricfunction $EXCITINGVISUAL/xmldielectricfunction2agr.xsl tmp.xml

Now, three xmgrace files are generated:

  1. dielectricfunction_Re.agr
  2. dielectricfunction_Im.agr
  3. dielectricfunction_ReKK.agr

which contain the real part, the imaginary part, and the real part from a Kramers-Kronig transform of the macroscopic dielectric function, respectively. This information is contained inside the EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT file in the second, third, and fourth column, respectively. To visualize the plot of the imaginary part with xmgrace (see also Xmgrace: A Quickstart), type:

$ xmgrace dielectricfunction_Im.agr >/dev/null 2>&1 &
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License