View source for Bethe-Salpeter equation tutorial. Optical absorption (BN)

In this tutorial we will deal with different aspects of running a Bethe-Salpeter (BSE) calculation for optical absorption spectra using yambopy. We continue using hexagonal boron nitride to illustrate how to perform convergence tests, and how to compare and analyse the results. 

Before you start this tutorial, make sure you have run the scf and nscf runs. If you have not, you can calculate the relax <code>-r</code>, scf <code>-s</code> and nscf <code>-n</code> using the <code>gs_bn.py</code> file:
 $ python gs_bn.py -r -s -n
Once that is over, you can start the tutorial.

__FORCETOC__

=== BSE convergence of optical spectra ===

In this section of the tutorial we will use the <code>bse_conv_bn.py</code> file. To evaluate the Bethe-Salpeter Kernel we need to first calculate the static dielectric screening, and then the screened Coulomb interaction matrix elements. This will be done with the file <code>bse_conv_bn.py</code>.

'''a. Static dielectric function'''

We begin by converging the static screening. In principle, all parameters can be converged independently one-by-one, and then the you can choose the best values for the final, fully converged result. To converge the static screening you will need:

<blockquote>

<code>FFTGvecs</code>: number of planewaves to include. Usually we need fewer planewaves than the total used in the self-consistent cycle that generated the ground state in the <code>scf</code> run. Reducing the total number of planewaves used diminishes the amount of memory per process that you are going to need. 

<code>BndsRnXs</code>: number of bands to calculate the screening. Typically you need many bands (hundreds or even more) for the screening to converge. 

<code>NGsBlkXs</code>: number of components for the local fields. Averages the value of the dielectric screening over a number of periodic copies of the unit cell. This parameter greatly increases the computational cost of the calculation and thus should be increased slowly.

</blockquote>

We create a dictionary with different values for each variable. The python script (<code>bse_conv_bn.py</code>) will then create a reference input file with the first value in each parameter's list. It will then create input files with the other parameters changing independently according to the values specified on the list:

 #list of variables to optimize the dielectric screening
 conv = { 'FFTGvecs': [[10,10,15,20,30],'Ry'],
          'NGsBlkXs': [[1,1,2,3,5,6], 'Ry'],
          'BndsRnXs': [[1,10],[1,10],[1,20],[1,30],[1,40]] }

To run a calculation for each of the variables in the dictionary, we first define the <code>run</code> function

    def run(filename):
        """
        Function to be called by the optimize function
        """
        path = filename.split('.')[0]
        print(filename, path)
        shell = scheduler()
        shell.add_command('cd %s'%folder)
        shell.add_command('%s mpirun -np %d %s -F %s -J %s -C %s 2> %s.log'%(nohup,threads,yambo,filename,path,path,path))
        shell.add_command('touch %s/done'%path)
        if not os.path.isfile("%s/%s/done"%(folder,path)):
            shell.run()

which together with the optimize function in the class <code>YamboIn</code>, 
 y.optimize(conv,folder='bse_conv',run=run,ref_run=False) 
will run and manage the calculations as defined in the dictionary. All calculation files and outputs will be in their respective directories inside the folder <code>bse_conv</code>. 

To check which options are available in the script <code>bse_conv_bn.py</code>, run 
 $ python bse_conv_bn.py
which will give you the following list:

 usage: bse_conv_bn.py [-h] [-r] [-a] [-p] [-e] [-b] [-u] [-t THREADS]
 
 Test the yambopy script.
 
 optional arguments:
   -h, --help            show this help message and exit
   -r, --run             run BSE convergence calculation
   -a, --analyse         analyse results data
   -p, --plot            plot the results
   -e, --epsilon         converge epsilon parameters
   -b, --bse             converge bse parameters
   -u, --nohup           run the commands with nohup
   -t THREADS, --threads THREADS
                         number of threads to use

So, in order to converge the screening, you will need to run

 $ python bse_conv_bn.py -r -e

By default we have set up the total number of threads to be equal to 2. If you have access to more, you can change this using the <code>-t</code> option shown on the list, followed by the number of threads you want to use. 

As you can see, the python script is running all the calculations changing the value of the input variables. You are free to open the <code>bse_conv_bn.py</code> file and modify it according to your own needs. Using the optimal parameters, you can run a calculation and save the dielectric screening databases <code>ndb.em1s*</code> to re-use them in the subsequent calculations. For that you can copy these files to the SAVE folder. <code>yambo</code> will only re-calculate any database if it does not find it or some parameter has changed.

Once the calculations are done you can plot resulting optical spectra by running 
 $ python bse_conv_bn.py -a 
followed by 
 $ python bse_conv_bn.py -p -e 
It will search the folder <code>bse_conv</code> for all calculations you performed previously and plot the results, using the class <code>YamboAnalyser</code> class.

[[File:Bse-hbn-conv-2.png|x750px|Yambo tutorial image]]

Alternatively, you can obtain the same plots by running the the following script: 
 $ python plot-bse-conv.py 
If you want to, you can now add new entries to the lists in the dictionary, specially the <code>BndsRnXs</code>, which usually requires a lot of bands to achieve convergence. Keep its number less than or equal to the number of bands in the nscf calculation, re-run the script and plot the results again. You can also test if some variables can converge with smaller values than the ones you used up to now. 

'''b. Screened Coulomb interaction'''

Up to this point we have been focused on the variables which control the static screening. We still need to deal with the variables which setup the Bethe-Salpeter auxiliary Hamiltonian matrix. Once this matrix is set up, <code>yambo</code> will use a diagonalisation algorithm to obtain the excitonic states and energies. The relevant variables for this process are:

<blockquote><code>BSEBands</code>: number of bands to generate the transitions. Should be as small as possible as the size of the BSE auxiliary hamiltonian has (in the resonant approximation) dimensions <code>Nk*Nv*Nc</code>. 

Another way to converge the number of transitions is using <code>BSEEhEny</code>. This variable selects the number of transitions based on the electron-hole energy difference (it's the one we are going to use).

<code>BSENGBlk</code> is the number of blocks for the dielectric screening average over the unit cells. This uses the static screening computed previously and controlled by the variable <code>NGsBlkXs</code>. So you <code>BSENGBlk</code> cannot be larger than <code>NGsBlkXs</code>.

<code>BSENGexx</code> in the number of exchange components. Relatively cheap to calculate, but should be kept as small as possible to save memory.

<code>KfnQP_E</code> is the scissor operator for the BSE. The first value is the rigid scissor, the second and third are the stretching for the conduction and valence respectively. The optical absorption spectrum is obtained in a range of energies given by <code>BEnRange</code> and the number of frequencies in the interval is <code>BEnSteps</code>. Alternatively a quasiparticle energy database can be included from a previous GW calculation.
</blockquote>

For the last variable you will be using a predefined scissor operator. If you want to know how to compute a scissor operator you can go and follow the GW tutorial here [[GW tutorial. Convergence and approximations (BN)]]. 

The dictionary of convergence in this case is:

 #list of variables to optimize the BSE
         conv = { 'BSEEhEny': [[[1,10],[1,10],[1,12],[1,14]],'eV'],
                  'BSENGBlk': [[0,0,1,2], 'Ry'],
                  'BSENGexx': [[10,10,15,20],'Ry']}

All these variables do not change the dielectric screening, so optionally you can calculate it once (using the previous section to determine converged parameters) and put the database in the <code>SAVE</code> folder to make the calculations faster. Otherwise the dielectric screening will be computed for each run (in the case of this tutorial, this is still quite fast). This is a slightly advanced use of <code>yambo</code>, so it is better to leave till you are more familiar with the code and its flow. 

To run the convergence calculations for the BSE Hamiltonian write:

 $ python bse_conv_bn.py -r -b
Once the calculations are done you can plot the optical absorption spectra:

 $ python bse_conv_bn.py -a
 $ python bse_conv_bn.py -p -b

[[File:Bse-hbn-conv-1.png|x750px|Yambo tutorial image]]

Alternatively, you can obtain the same plots by running the the following script: <source>$ python plot-bse-conv.py </source>
Again, we invite you to change the values in the dictionary to check if the calculation would converge better using larger or smaller parameters.

=== Coulomb truncation convergence ===

'''a. Running for different layer separation distances'''

Here we will check how the dielectric screening changes with vacuum spacing between layers and including a coulomb truncation technique. For that we define a loop where we do a self-consistent ground state calculation, non self-consistent calculation, create the databases and run a <code>yambo</code> BSE calculation for different vacuum spacings.

To analyze the data we will plot the dielectric screening and check how the different values of the screening change the absorption spectra.
In the folder <code>tutorials/bn/</code> you find the python script <code>bse_cutoff_bn.py</code>. This script takes some time to be executed, you can run both variants without the cutoff and with the cutoff <code>-c</code> simultaneously to save time. You can run this script with:

 python bse_cutoff_bn.py -r -t4    # without coulomb cutoff
 python bse_cutoff_bn.py -r -c -t4 # with coulomb cutoff
where <code>-t</code> specifies the number of MPI threads to use. The main loop changes the <code>layer_separation</code> variable using values from a list in the header of the file. In the script you can find how the functions <code>scf</code>, <code>ncf</code> and <code>database</code> are defined.

'''b. Plot the dielectric function'''

In a similar way as what was done before we can now plot the dielectric function for different layer separations:

 yambopy plotem1s bse_cutoff/*/*     # without coulomb cutoff  
 yambopy plotem1s bse_cutoff_cut/*/* # with coulomb cutoff

[[File:Bn em1s cutoff cut.png|Yambo tutorial image]]

[[File:Bn em1s cutoff.png|Yambo tutorial image]]

In these figures it is clear that the long-range part of the coulomb interaction (q=0 in reciprocal space) is truncated, i. e. it is forced to go to zero.

'''c. Plot the absorption spectrum'''

You can also plot how the absorption spectra changes with the cutoff using:

 python bse_cutoff_bn.py -p
 python bse_cutoff_bn.py -p -c

[[File:Bn bse cutoff.png|bn bse cutoff]]

[[File:Bn bse cutoff cut.png|bn bse cutoff cut]]

As you can see, the spectra is still changing with the vacuum spacing, you should increase the vacuum until convergence. For that you can add larger values to the <code>layer_separations</code> list and run the calculations and analysis again.

=== Excitonic wavefunctions in reciprocal space ===

In this example we show how to use <code>yambopy</code> to plot the excitonic wavefunction weights that result from a BSE calculation (specifically from the diagonalisation of the excitonic Hamiltonian). The script we will use this time is: <code>bse_bn.py</code>. Be aware the parameters specified for the calculation are not high enough to obtain a converged result. 

To run the BSE calculation do:

 python bse_bn.py -r
Alternatively, you can run  
 python bse_bn.py -r -c  
to add the Coulomb cutoff.

Then, the absorption spectrum is readily plotted by running the auxiliary script
 python plot-bse.py
 
[[File:bse_abs_plt.png|Yambo tutorial image]]

Afterwards you can run a basic analysis of the excitonic states and store the wavefunctions of the ones that are the most optically active, plotting them in reciprocal space. Plots in real space (not covered here) are also possible using yambopy (by calling ypp). 

In order to proceed we will use the script <code>plot-excitondb.py</code>.
After running it as 

 python plot-excitondb.py
 
it will produce three plots.

The first one is a representation of the excitonic weights (or rather their modulus) in reciprocal space, band- and kpoint-resolved, for the lowest-bound exciton.
You can see that the electronic transitions contributing to this state originate mostly at symmetry point K (direct band gap) and along the KM line (the bands here have pi-type orbital character).

[[File:exc_bands.png|Yambo tutorial image]]

The second figure is an interpolated version of the first, with the bands and excitonic contributions appearing smoother.

[[File:exc_bands_interp.png|Yambo tutorial image]]


The final figure is a representation of the excitonic weights, summed over all bands, within the Brillouin zone. The brighter colors indicate where the most important electronic transitions for this excitonic state are located.

[[File:exc_BZ.png|Yambo tutorial image]]

You can obtain plots for different excitonic states by changing the line

 # List of states to be merged
 states = [1,2]
 
in <code>plot-excitondb.py</code>.

You may also want to test the script with a more converged calculation: to this end you can increase the number of kpoints in the nscf calculation by going to <code>gs_bn.py</code> and changing the line

 kpoints_nscf = [6,6,1]

to your liking (for example to <code>12x12x1</code>).

Then, run again the following commands (nscf calculation, regeneration of yambo database, new BSE calculation).

 python gs_bn.py -n
 rm -r database
 python bse_bn.py -r -t 4

As for the plotting, you may want to increase the value of <code>lpratio</code> in <code>plot-excitondb.py</code> to ensure that the interpolation works, as well as reduce the value of <code>scale</code> for a nicer plot in the hexagonal BN.
The lines to be changed are

 exc_bands_inter = yexc.interpolate(save,path,states,lpratio=5,f=None,size=0.5,verbose=True)

and

 yexc.plot_exciton_2D_ax(ax,states,mode='hexagon',limfactor=0.8,scale=160)

Then, you can run <code>plot-excitondb.py</code> again.

Again, be aware that these figures serve only to show the kind of representation that can be obtained with <code>yambo</code>, <code>ypp</code> and <code>yambopy</code>. Further convergence tests need to be performed to obtain accurate results, but that is left to the user.

<!--
You can now visualise these wavefunctions in real space using our online tool: [http://henriquemiranda.github.io/excitonwebsite/ http://henriquemiranda.github.io/excitonwebsite/]

For that, go to the website, and in the <code>Excitons</code> section select <code>absorptionspectra.json</code> file using the <code>Custom File</code>. You should see on the right part the absorption spectra and on the left the representation of the wavefunction in real space. Alternatively you can vizualise the individually generated <code>.xsf</code> files using xcrysden.
-->

<!--
=== Parallel static screening ===

In this tutorial we will show how you can split the calculation of the dielectric function in different jobs using <code>yambopy</code>. The dielectric function can then be used to calculate the excitonic states using the BSE.

The idea is that in certain clusters it is advantageous to split the jobs as much as possible. The dielectric function is calculated for different momentum transfer (q-points) over the brillouin zone. Each calculation is independent and can run at the same time. Using the <code>yambo</code> parallelization you can separate the dielectric function calculation among many cpus using the variable <code>q</code> in <code>X_all_q_CPU</code> and <code>X_all_q_ROLEs</code>. The issue is that you still need to make a big reservation and in some cases there is load imbalance (some nodes end up waiting for others). Splitting in smaller jobs can help your jobs to get ahead in the queue and avoid the load imbalance. If there are many free nodes you might end up running all the q-points at the same time.

'''The idea is quite simple:''' you create an individual input file for each q-point, submit each job separately, collect the results and do the final BSE step (this method should also apply for a GW calculation).

'''2. Parallel Dielectric function'''

To run the dielectric function in parallel do:

<source lang="bash">python bse_par_bn.py -r -t2</source>
Here we tell <code>yambo</code> to calculate the dielectric function. We read the number of q-points the system has and generate one input file per q-point. Next we tell <code>yambo</code> to calculate the first q-point. <code>yambo</code> will calculate the dipoles and the dielectric function at the first q-point. Once the calculation is done we copy the dipoles to the SAVE directory. After that we run each q-point calculation as a separate job. Here the user can decide to submit one job per q-point on a cluster or use the python <code>multiprocessing</code> module to submit the jobs in parallel. In this example we use the second option.

<source lang="python">from yambopy import *
import os
import multiprocessing

yambo = "yambo";
folder = "bse_par";
nthreads = 2 #create two simultaneous jobs

#create the yambo input file
y = YamboIn('yambo -r -b -o b -V all',folder=folder)

y['FFTGvecs'] = [30,'Ry']
y['NGsBlkXs'] = [1,'Ry']
y['BndsRnXs'] = [[1,30],'']
y.write('%s/yambo_run.in'%folder)

#get the number of q-points
startk,endk = map(int,y['QpntsRXs'][0])

#prepare the q-points input files
jobs = []
for nk in xrange(1,endk+1):
    y['QpntsRXs'] = [[nk,nk],'']
    y.write('%s/yambo_q%d.in'%(folder,nk))
    if nk != 1:
        jobs.append('cd %s; %s -F yambo_q%d.in -J yambo_q%d -C yambo_q%d 2&gt; log%d'%(folder,yambo,nk,nk,nk,nk))

#calculate first q-point and dipoles
os.system('cd %s; %s -F yambo_q1.in -J yambo_q1 -C yambo_q1'%(folder,yambo))
#copy dipoles to save
os.system('cp %s/yambo_q1/ndb.dip* %s/SAVE'%(folder,folder))

p = multiprocessing.Pool(nthreads)
p.map(run_job, jobs)</source>
'''3. BSE'''

Once the dielectric function is calculated, it is time to collect the data in one folder and do the last step of the calculation: generate the BSE Hamiltonian, diagonalize it and calculate the absorption.

<source lang="python">#gather all the files
if not os.path.isdir('%s/yambo'%folder):
    os.mkdir('%s/yambo'%folder)
os.system('cp %s/yambo_q1/ndb.em* %s/yambo'%(folder,folder))
os.system('cp %s/*/ndb.em*_fragment* %s/yambo'%(folder,folder))

y = YamboIn('yambo -r -b -o b -k sex -y d -V all',folder=folder)
y['FFTGvecs'] = [30,'Ry']
y['NGsBlkXs'] = [1,'Ry']
y['BndsRnXs'] = [[1,30],'']
y['BSEBands'] = [[3,6],'']
y['BEnSteps'] = [500,'']
y['BEnRange'] = [[0.0,10.0],'eV']
y['KfnQP_E']  = [2.91355133,1.0,1.0] #some scissor shift
y.arguments.append('WRbsWF')
y.write('%s/yambo_run.in'%folder)

print('running yambo')
os.system('cd %s; %s -F yambo_run.in -J yambo'%(folder,yambo))</source>
'''3. Collect and plot the results'''

You can then plot the data as before:

<source lang="bash">python bse_par_bn.py -p</source>
This will execute the following code:

<source lang="python">#collect the data
pack_files_in_folder('bse_par')

#plot the results using yambo analyser
y = YamboAnalyser()
print y
y.plot_bse(['eps','diago'])</source>
You should obtain a plot like this:

[[File:Bethe Salpeter in BN.png|Yambo tutorial image]]
-->