gen_meff.py

#!/usr/bin/env python3

import sys
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import hbar, m_e, e
import re 


my_path = os.path.abspath("/work/i254/i254/gentles/binaries_4D_2")
if my_path not in sys.path:
    sys.path.insert(0, my_path)
my_path = os.path.abspath("/work/i254/i254/gentles/tools/pp")
if my_path not in sys.path:
    sys.path.insert(0, my_path)

from get_eigs_xml import get_eigs_xml
from bandgap_qescf import qe_scf_lattice_param_2

def lattice_constant(crysfile):
    p1=re.compile(r"alat=([0-9.]+)")
    with open(crysfile, 'r') as fi:
        s=fi.read()
        q1=p1.search(s)
        if q1:
            a=float(q1.group(1))*0.529177
        else:
            a=None 
    return a

def effective_mass(delta_E_eV, delta_k_units, lattice_constant_m):
    """
    Calculate the effective mass from energy and k-point difference.
    
    Parameters:
        delta_E_eV (float): Energy difference (E2 - E1) in electronvolts (eV)
        delta_k_units (float): Difference in k-points (unitless, in units of 2π/a)
        lattice_constant_m (float): Lattice constant 'a' in meters

    Returns:
        float: Effective mass in units of the electron mass (m*/m_e)
    """
    # Convert energy from eV to Joules
    delta_E_J = delta_E_eV * e
    
    # Convert k difference to actual k-value in 1/m
    delta_k = delta_k_units * (2 * np.pi / lattice_constant_m)
    
    # Effective mass formula
    m_star = hbar**2 / (delta_E_J / delta_k**2)
    
    # Return in units of electron mass
    return 0.5*m_star / m_e


def specific_meff(base):
    ks,eigs,_=get_eigs_xml(f"dir_{base}/tmp_{base}/{base}.xml")
    a=lattice_constant(f"dir_{base}/{base}.out.crys")
    return effective_mass(eigs[1,448]-eigs[0,448],np.linalg.norm(ks[1]),a*10**(-10))

def specific_meff_calc(base,i):
    ks,eigs,_=get_eigs_xml(f"dir_{base}/calc_{i}/meff_{base}.{i}/{base}.{i}.xml")
    a=lattice_constant(f"dir_{base}/calc_{i}/{base}.{i}.out.crys")
    return effective_mass(eigs[1,448]-eigs[0,448],np.linalg.norm(ks[1]),a*10**(-10))