1ph IP (and excph eventually) resonant raman scattering

tutorial/exciton-phonon/ip_raman.py

@@ -0,0 +1,116 @@
+"""
+IP one-phonon Stokes Raman with yambopy.
+
+The driver does all the band-window bookkeeping:
+  - loads lattice, electrons, dipoles, ndb.elph
+  - BZ-expands the dipoles via YamboDipolesDB(expand=True) and reorders
+    them into the LetzElPhC k-grid
+  - slices electron energies, dipoles and gkkp to the user-chosen
+    `raman_bands` window (Fortran 1-indexed, inclusive)
+  - passes ready-to-use arrays to ip_resonant_raman_tensor_oneph
+
+The Raman function itself takes only the energies, dipoles, gkkp,
+cell volume and broadening and evaluates the formula.
+
+Assumes yambo wrote velocity-gauge dipoles (DIP_v in ndb.dipoles).
+"""
+
+import numpy as np
+
+from yambopy import (YamboLatticeDB, YamboElectronsDB,
+                     YamboDipolesDB, LetzElphElectronPhononDB)
+from yambopy.exciton_phonon.excph_resonant_raman import ip_resonant_raman_tensor_oneph
+
+
+# ---- User inputs ----------------------------------------------------------
+folder      = '.'
+SAVE_dir    = folder + '/SAVE'
+elph_path   = folder + '/ndb.elph'
+dipole_path = folder + '/dipoles/ndb.dipoles'
+
+omega_range = (1.0, 6.0, 501)         # laser energies (eV): min, max, npts
+broading    = 0.10                    # Lorentzian broadening (eV)
+ph_fre_th   = 5.0                     # acoustic-mode cutoff (cm^-1)
+raman_bands = (44, 56)                # 1-indexed Fortran inclusive
+
+
+# ---- Databases ------------------------------------------------------------
+lattice   = YamboLatticeDB.from_db_file(filename='%s/ns.db1' % SAVE_dir)
+electrons = YamboElectronsDB.from_db_file(folder=SAVE_dir,
+                                          filename='ns.db1', Expand=True)
+dipdb     = YamboDipolesDB.from_db_file(lattice,
+                                        filename=dipole_path,
+                                        dip_type='v',
+                                        project=False,
+                                        expand=True)
+elph      = LetzElphElectronPhononDB(elph_path, read_all=True)
+
+
+# ---- Yambo BZ -> elph BZ permutation --------------------------------------
+diff = elph.kpoints[:, None, :] - lattice.red_kpoints[None, :, :]
+diff = diff - np.round(diff)
+perm = np.argmin(np.linalg.norm(diff, axis=-1), axis=1)
+
+
+# ---- Resolve band range (1-indexed inclusive, Fortran) -------------------
+b_in_F, b_out_F = raman_bands
+nb              = b_out_F - b_in_F + 1
+n_val           = int(electrons.nbandsv) - (b_in_F - 1)
+nc              = nb - n_val
+cell_vol        = float(lattice.lat_vol)
+
+
+# ---- eph_g(q=0) sliced to the raman window --------------------------------
+elph_b_in_F = int(min(elph.bands))
+off         = b_in_F - elph_b_in_F
+
+iq0   = 0
+g_q0  = elph.gkkp[iq0, :, :, 0, :, :]                       # (nk, nmodes, nb_i, nb_f)
+g_q0  = np.swapaxes(g_q0, -1, -2)
+eph_g = np.transpose(g_q0, (1, 0, 2, 3)) / 2.0              # (nmodes, nk, nb_eph, nb_eph), Ha
+eph_g = eph_g[:, :, off : off + nb, off : off + nb]         # (nmodes, nk, nb, nb)
+
+ph_eV = elph.ph_energies[iq0]                               # (nmodes,), eV
+
+
+# ---- KS energies in elph k-order, sliced to the raman window --------------
+el_eV = electrons.eigenvalues[0][perm, b_in_F - 1 : b_out_F]   # (nk, nb)
+
+
+# ---- Velocity dipoles: cv block, reorder, slice to raman cv block --------
+nv_dip       = int(dipdb.nbandsv)
+nc_dip       = int(dipdb.nbandsc)
+elec_dipoles = dipdb.dipoles[:, :, nv_dip : nv_dip + nc_dip, :nv_dip]   # (nk_y, 3, nc_dip, nv_dip)
+elec_dipoles = elec_dipoles[perm]                                       # -> elph k-order
+elec_dipoles = elec_dipoles[:, :, :nc, -n_val:]                         # raman cv block
+elec_dipoles = np.transpose(elec_dipoles, (1, 0, 2, 3))                 # (3, nk, nc, n_val)
+
+
+# ---- Raman calculation ----------------------------------------------------
+laser_eV = np.linspace(omega_range[0], omega_range[1], int(omega_range[2]))
+
+R = ip_resonant_raman_tensor_oneph(
+    laser_energies    = laser_eV,
+    ph_energies       = ph_eV,
+    el_energies       = el_eV,
+    elec_dipoles      = elec_dipoles,
+    eph_g             = eph_g,
+    cell_vol          = cell_vol,
+    broad             = broading,
+    ph_freq_threshold = ph_fre_th,
+)
+
+
+# ---- Save -----------------------------------------------------------------
+raman_inte = np.sum(np.abs(R)**2, axis=(1, 2, 3))
+
+np.savetxt('raman_intensity.dat',
+           np.column_stack([laser_eV, raman_inte]),
+           header='laser_energy(eV)   raman_intensity(arb.u.)')
+
+np.savez('Raman_tensors.npz',
+         laser_energies_eV = laser_eV,
+         ph_freq_cm        = ph_eV * 8065.54,
+         raman_tensor      = R)
+
+print('Wrote raman_intensity.dat and Raman_tensors.npz')

yambopy/dbs/excphondb.py

@@ -68,8 +68,8 @@ def __init__(self,lattice,save_excph="./",read_all=True):
         self.nexc_o = database.variables['EXCITON_SUM'][1].astype(int)
         self.nmodes = database.variables['PHONON_MODES'][0].astype(int)
         self.nqpoints = database.variables['HEAD_R_LATT'][3].astype(int)
-        self.type_exc_i = database.variables['L_kind_in'][...].tostring().decode().strip()
-        self.type_exc_o = database.variables['L_kind_out'][...].tostring().decode().strip()
+        self.type_exc_i = database.variables['L_kind_in'][...].tobytes().decode().strip()
+        self.type_exc_o = database.variables['L_kind_out'][...].tobytes().decode().strip()
         database.close()
 
         #Check how many databases are present
@@ -113,33 +113,65 @@ def read_qpoints(self):
 
         self.car_qpoints = np.array([ q/self.alat for q in self.qpoints ])
 
-    def read_excph(self):
+    def read_excph(self, read_sq=None):
         """
-        Read exciton-phonon matrix elements and their modulus squared
-
+        Read exciton-phonon matrix elements and (optionally) their modulus
+        squared.
+
         NB: EXCPH_GKKP_Q is saved by yambo as (2,mode,exc_out,exc_in), but netCDF stores
             the *transpose* (exc_in,exc_out,mode,2).
             We want to change it to complex (iq,mode,exc_in,exc_out)
-        """    
+
+        Parameters
+        ----------
+        read_sq : None or bool, optional
+            How to obtain the modulus-squared array
+            EXCITON_PH_GKKP_SQUARED_Q. Newer lumen databases no longer store
+            this variable, so ``self.excph_sq`` is always populated for
+            backward compatibility with consumers that expect it:
+
+            * None (default) -- auto: read the stored array when present,
+              otherwise reconstruct it on the fly as ``|excph|**2``.
+            * True  -- read from the DB (errors if the variable is absent).
+            * False -- always reconstruct as ``|excph|**2`` (ignore any
+              stored array).
+
+            The reconstruction ``REAL(g)**2 + AIMAG(g)**2`` is exactly the
+            expression lumen used to fill the stored array, so the result is
+            equivalent up to floating-point precision.
+        """
         var_nm    = "EXCITON_PH_GKKP_Q"
         var_sq_nm = "EXCITON_PH_GKKP_SQUARED_Q"
-            
+
         # excph[q][mode][iexc1][iexc2]
-        excph_full    = np.zeros([self.nfrags,self.nmodes,self.nexc_i,self.nexc_o],dtype=np.complex64)  
-        excph_sq_full = np.zeros([self.nfrags,self.nmodes,self.nexc_i,self.nexc_o])  
-
+        excph_full    = np.zeros([self.nfrags,self.nmodes,self.nexc_i,self.nexc_o],dtype=np.complex64)
+        excph_sq_full = np.zeros([self.nfrags,self.nmodes,self.nexc_i,self.nexc_o])
+
+        # Decide whether to read the stored squared array (absent in newer
+        # lumen databases, which dropped EXCPH_Gkkp_sq).
+        if read_sq is None:
+            with Dataset(self.frag_filename + "1") as db0:
+                read_sq = ('%s1'%var_sq_nm) in db0.variables
+
         for iq in range(self.nfrags):
             fil = self.frag_filename + "%d"%(iq+1)
             database = Dataset(fil)
             excph = database.variables['%s%d'%(var_nm,iq+1)][:]
             excph_full[iq] = np.moveaxis( excph[:,:,:,0]+I*excph[:,:,:,1], -1,0 )
             #excph_full[iq] = np.swapaxes( np.swapaxes(excph[:,:,:,0] + I*excph[:,:,:,1],-1,0), -1,-2)
-
-            excph_sq = database.variables['%s%d'%(var_sq_nm,iq+1)][:]
-            excph_sq_full[iq] = np.moveaxis( excph_sq, -1,0)
-            #excph_sq_full[iq] = np.swapaxes( np.swapaxes(excph_sq[:,:,:],-1,0), -1,-2)
+
+            if read_sq:
+                excph_sq = database.variables['%s%d'%(var_sq_nm,iq+1)][:]
+                excph_sq_full[iq] = np.moveaxis( excph_sq, -1,0)
+                #excph_sq_full[iq] = np.swapaxes( np.swapaxes(excph_sq[:,:,:],-1,0), -1,-2)
             database.close()
-
+
+        # If the stored squared array is not available, reconstruct it on the
+        # fly so self.excph_sq stays a valid array for downstream consumers.
+        # In-place assignment keeps the float64 dtype of the read path.
+        if not read_sq:
+            excph_sq_full[:] = np.abs(excph_full)**2
+
         # Check integrity of elph values
         if np.isnan(excph_full).any(): print('[WARNING] NaN values detected in elph database.')