#***********************************************************
# * Copyright (C) 2017-2019 Brendan A. Smith and Alexey V. Akimov
# * This file is distributed under the terms of the
# * GNU General Public License as published by the
# * Free Software Foundation; either version 3 of the
# * License, or (at your option) any later version.
# * http://www.gnu.org/copyleft/gpl.txt
#***********************************************************/
"""
.. module:: compute_hprime
:platform: Unix, Windows
:synopsis: This module computes the transition dipole moments in the basis of KS
orbitals extracted from QE
.. moduleauthor:: Brendan A. Smith and Alexey V. Akimov
"""
import os
import sys
# Fisrt, we add the location of the library to test to the PYTHON path
if sys.platform=="cygwin":
from cyglibra_core import *
elif sys.platform=="linux" or sys.platform=="linux2":
from liblibra_core import *
[docs]def compute_hprime_dia(es, info, filename):
"""Computes the matrix elements of transition dipole moments
This function computes the matrix elements of the dipole operator for the
case without SOC, and prints them to them to the file specificed by the
function parameter "filename"
Args:
es ( dictionary ): Information about electronic structure.
* **es["Coeff_dia"]** ( list of CMATRIX(npw, len(act_space)) objects ) : the
wavefunction coefficients in the planewave basis for the spin-diabatic wavefunctions, such
that res_curr["Coeff_dia"][k] is a matrix for the k-point with index k.
Can be generated by ..seealso::step2.read_wfc_grid
* **es["grid"]** ( list of VECTOR objects ): the grid point vectors [ units: tpiba ]
Can be generated by ..seealso::step2.read_wfc_grid
info ( dictionary ): The basic information regarding the system. ..seealso::QE_methods.read_qe_index
step2.read_info to see how this object can be generated. For the purpose of this calculation, it should include:
* **info["b1"]** ( VECTOR ): reciprocal lattice vector 1 [units: Bohr^-1]
* **info["b2"]** ( VECTOR ): reciprocal lattice vector 2 [units: Bohr^-1]
* **info["b3"]** ( VECTOR ): reciprocal lattice vector 3 [units: Bohr^-1]
filename ( string ): This is the name of the output file where the data will be printed out
Returns:
[CMATRIX(N/2, N/2) , CMATRIX(N/2, N/2), CMATRIX(N/2, N/2) ]: Hprime_x, Hprime_y, Hprime_z:
the matrices with i*hbar*<i|r_alpha|j> matrix elements, where N = len(act_space) for alpha = x, y, z
Here, the act_space includes both the alpha and beta-orbitals
"""
coeff = es["Coeff_dia"]
grid = es["grid"][0]
b1 = info["b1"]
b2 = info["b2"]
b3 = info["b3"]
g = MATRIX(len(grid),3)
for i in range(0,len(grid)):
g.set(i,0,grid[i].x)
g.set(i,1,grid[i].y)
g.set(i,2,grid[i].z)
reci = MATRIX(3,3)
reci.set(0,0,b1.x); reci.set(0,1,b2.x); reci.set(0,2,b3.x)
reci.set(1,0,b1.y); reci.set(1,1,b2.y); reci.set(1,2,b3.y)
reci.set(2,0,b1.z); reci.set(2,1,b2.z); reci.set(2,2,b3.z)
Hprime = compute_Hprime(coeff,g,reci)
Hprime[0].real().show_matrix("%sx_re" % (filename)); Hprime[0].imag().show_matrix("%sx_im" % (filename))
Hprime[1].real().show_matrix("%sy_re" % (filename)); Hprime[1].imag().show_matrix("%sy_im" % (filename))
Hprime[2].real().show_matrix("%sz_re" % (filename)); Hprime[2].imag().show_matrix("%sz_im" % (filename))
return Hprime
[docs]def hprime_py(es, info, filename):
"""Computes the matrix elements of transition dipole moments
This function computes the matrix elements of the dipole operator for the
case without SOC, and prints them to them to the file specificed by the
function parameter "filename"
This is a Python-only version of the ::funct:```compute_hprime_dia```, so it is going
to be slower than that version, but it allows more flexibility
Args:
es ( dictionary ): Information about electronic structure.
* **es["Coeff_dia"]** ( list of CMATRIX(npw, len(act_space)) objects ) : the
wavefunction coefficients in the planewave basis for the spin-diabatic wavefunctions, such
that res_curr["Coeff_dia"][k] is a matrix for the k-point with index k.
Can be generated by ..seealso::step2.read_wfc_grid
* **es["grid"]** ( list of VECTOR objects ): the grid point vectors [ units: tpiba ]
Can be generated by ..seealso::step2.read_wfc_grid
info ( dictionary ): The basic information regarding the system. ..seealso::QE_methods.read_qe_index
step2.read_info to see how this object can be generated. For the purpose of this calculation, it should include:
* **info["b1"]** ( VECTOR ): reciprocal lattice vector 1 [units: Bohr^-1]
* **info["b2"]** ( VECTOR ): reciprocal lattice vector 2 [units: Bohr^-1]
* **info["b3"]** ( VECTOR ): reciprocal lattice vector 3 [units: Bohr^-1]
filename ( string ): This is the name of the output file where the data will be printed out
Returns:
[CMATRIX(N/2, N/2) , CMATRIX(N/2, N/2), CMATRIX(N/2, N/2) ]: Hprime_x, Hprime_y, Hprime_z:
the matrices with i*hbar*<i|r_alpha|j> matrix elements, where N = len(act_space) for alpha = x, y, z
Here, the act_space includes both the alpha and beta-orbitals
"""
coeff = es["Coeff_dia"][0]
grid = es["grid"][0]
b1 = info["b1"]
b2 = info["b2"]
b3 = info["b3"]
g_sz = len(grid)
npw = coeff.num_of_rows
print("g_sz = ", g_sz)
print("npw = ", npw)
scl1 = (1.0+0.0j)
scl2 = 2.0*(0.0+1.0j)
# figure out if we are using a completed wavefunction or not
if g_sz != npw:
g_sz = min(g_sz,npw)
is_compl = 1
else:
is_compl = 0
N = coeff.num_of_cols
sz = int(N/2)
Hx = CMATRIX(sz, sz); Hy = CMATRIX(sz, sz); Hz = CMATRIX(sz, sz)
for g in range(0,g_sz):
gx = scl1*(grid[g].x*b1.x + grid[g].y*b2.x + grid[g].z*b3.x)
gy = scl1*(grid[g].x*b1.y + grid[g].y*b2.y + grid[g].z*b3.y)
gz = scl1*(grid[g].x*b1.z + grid[g].y*b2.z + grid[g].z*b3.z)
for i in range(0, sz):
for j in range(0, sz):
tmp = coeff.H().get(i,g) * coeff.get(g,j)
hx, hy, hz = 0.0, 0.0, 0.0
if(is_compl==0):
hx = tmp.real*gx; hy = tmp.real*gy; hz = tmp.real*gz;
if(is_compl==1):
if(g==0):
hx = tmp.real*gx; hy = tmp.real*gy; hz = tmp.real*gz; # This should give zero!
#Now the Hprime_ matrices are purely imaginary, for the case of gamma-symmetry.
else:
hx = scl2*tmp.real*gx
hy = scl2*tmp.real*gy
hz = scl2*tmp.real*gz
Hx.add(i,j,hx)
Hy.add(i,j,hy)
Hz.add(i,j,hz)
Hx.real().show_matrix("%sx_re" % (filename)); Hx.imag().show_matrix("%sx_im" % (filename))
Hy.real().show_matrix("%sy_re" % (filename)); Hy.imag().show_matrix("%sy_im" % (filename))
Hz.real().show_matrix("%sz_re" % (filename)); Hz.imag().show_matrix("%sz_im" % (filename))
return [Hx, Hy, Hz]