Absorption spectra by Fermi golden Rule
There is a possibility to evaluate rough estimation of UV-VIS absorption spectra by Fermi golden rule between Kohn-Sham orbitals. For a given range of photon energy are selected all pairs of occupied and unoccupied ground-state Kohn-Sham orbital with appropriate energy difference. This two orbitals are than projected onto real space grid and transition dipole moment is computed. The Oscillator strength (or probability of transition) is than proportional to square of length of transition dipole moment. The Absorption spectra can be than plotted by summing up contributions of oscillator strength for occupied and unoccupied orbital pairs with energy difference falling into particular energy bins (like weighted histogram).
Note that this method ( Fermi golden rule between ground state single particle states obtained from LDA-DFT ) is a very crude approximation of the true optical transitions. No beyond-DFT nor beyond-Hratree-Fock correction (such as RPA, TD-DFT, GW, BSE … ) are used. For this reason the transition energy ( position of peak ) might not be realistic. The method is, however, useful to explore the dependence of electronic transition probability on symmetry breaking and the changes in shape and localization of electronic states.
Availability
this is implemented only in Prokop's personal versions of Fireball located in
/data/home/hapala/Fireball_dev/src_1.0-BSfinal /data/home/hapala/Fireball_dev/progs_Jellium_mod
Usage
The computation of optical transitions is activated by setting keyword iwrtexcit = 1 in firebal.in. However, currently it is also necessary to activate iwrtewf = 1 in order to initialize a grid, even thought .xsf files written by iwrtewf option are not used at all and have nothing to do with the computation of transitions ( this should be corrected in future). So the typical fireball.in looks like:
&OPTION basisfile = answer.bas nstepf = 1 icluster = 1 ifixcharge = 1 &END &OUTPUT iwrtewf = 1 iwrtexcit=1 &END &MESH iewform = 3 npbands = 1 pbands = 1 &END
Than excitations.optional file should be provided to
0 ! write_out_transition_map ( default 0, moslty just for debuging ) 0.0 ! dEmin ... minimum energy difference ( mostly 0) 3.5 ! dEmax ... maximim energy difference
if write_out_transition_map is set to '1', an '.xsf' file containing spatial distribution of transition dipole moment is outputted for each pair of orbitals. This may be useful to identify parts of molecule which contribute mostly to the electronic transition. However, the number of transitions grow quadratically with the energy range and each .xsf file needs typically more than 100MB of disk space. So in most of cases you would quickly consume HUGE amount of disk space ( like 1 terrabyte ) if you switch this option on.
It is also recommended to chose reasonable small energy interval (dEmin,dEmax) otherwise the computation can take very long time.
Outputs
Excitations.dat file is written out from the calculation. The file contains several columns of following meaning:
j,jindex of occupied and unoccupied orbitaldelta_Eijenergy difference|D|length of transition dipole momentD_x,D_y,D_zcartesian components of transition dipole momentEi,Ejenergy eigenvalues corresponding to the Khon-Sham orbitals
i j delta_Eij |D| D_x D_y D_z Ei Ej 749 879 3.4862376070 0.9824594517 0.2759973643 -0.9336816482 0.1314937604 -3.0970434024 -6.5832810094 750 879 3.4791174668 0.9753364513 -0.4413174043 0.2404091902 0.8358968616 -3.0970434024 -6.5761608692 751 879 3.4582785159 0.5954896260 0.3594301771 -0.1688715641 -0.4437344220 -3.0970434024 -6.5553219183 752 879 3.4453548553 0.8057311791 -0.1487142116 -0.4625862731 -0.6427291469 -3.0970434024 -6.5423982578 753 879 3.4410779875 0.7684644226 0.0080712081 -0.6875508088 -0.3431418216 -3.0970434024 -6.5381213900 754 879 3.4257564875 0.4202213606 0.0892087920 -0.1778280360 -0.3701418282 -3.0970434024 -6.5227998899
Making spectra
For transformation of pairwise table of transition dipole moment into optical spectra for plotting can be used for example following bash/awk scripts.
First script spectranice.sh accumulate (sum up) pairwise contributions of transition dipole moment into bins of given size and then smear them with gaussian filter (to make nice smooth spectra). The most important parameter is step size. For example with step size 0.001 eV:
spectranice.sh 0 3.5 0.001
spectranice.sh code:
min=$1
max=$2
step=$3
awk -v min=$min -v max=$max -v step=$step '
BEGIN{for(i=0;min+i*step<max;i++){hist[i]=0}; minx=10000000; maxx=-10000000 }
(NR>0){
x=$3; y=$4; yy=y*y;
ei=int((x-min)/step)
if(x>maxx){maxx=x}
if(x<minx){minx=x}
hist[ei-5]+=yy*0.00098
hist[ei-4]+=yy*0.00977
hist[ei-3]+=yy*0.04395
hist[ei-2]+=yy*0.11719
hist[ei-1]+=yy*0.20508
hist[ei+0]+=yy*0.24609
hist[ei+1]+=yy*0.20508
hist[ei+2]+=yy*0.11719
hist[ei+3]+=yy*0.04395
hist[ei+4]+=yy*0.00977
hist[ei+5]+=yy*0.00098
}
END{
#print maxx,minx
for(i=0;min+i*step<max;i++){
if(((min+i*step)>minx) && ((min+i*step)<maxx) ) {print min+i*step,hist[i]}}
}
' Excitations.dat > spectranice.dat
in some cases you are not interested in the sum of all contributions ( this would correspond to absorption spectra) but rather in the transition with highest oscillator strength in given interval of energy (energy bin). For this purpose script 'spectramax.sh can be used.
spectramax.sh 0.001
#! /bin/bash
step=$1
min=`awk '
BEGIN{Emin=100000}
(NR>0){if($3<Emin)Emin=$3}
END{print Emin}
' Excitations.dat`
max=`awk '
BEGIN{Emax=0}
(NR>0){if($3>Emax)Emax=$3}
END{print Emax}
' Excitations.dat`
echo $min $max
awk -v min=$min -v max=$max -v step=$step '
BEGIN{for(i=0;min+i*step<max;i++){hist[i]=0;num[i]=0;}}
(NR>0){
ei=int(($3-min)/step)
if(hist[ei]<($4*$4)){hist[ei]=($4*$4)}
}
END{for(i=0;min+i*step<max;i++){print min+i*step,hist[i]}}
' Excitations.dat > spectramax.dat
