Ground state analysis for NiO

[mpmdean]

We should probably re-do the Ground state analysis for NiO

import edrixs
import numpy as np
import matplotlib.pyplot as plt

def get_fortran_eigvec(eigvec_num):
    '''
    read eigenvector for nth (n = eigvec_num) eigenvalue from files produced by the Fortran solver.
    '''
    f = open('files/eigvec.'+str(eigvec_num), 'rb')
    dt = np.dtype(np.complex128)
    ffile  = np.fromfile( f, dtype=dt, offset = 4)
    eigval = ffile[0:1].real[0] # first complex  number is the eigenvalue
    v_for  = ffile[1:]          # the rest of complex number are the eigenvector values
    f.close()

    return v_for
    
def decimalToBinary(n, norb):
    '''
    convert a decimal number to its binary form.
    norb is the number of orbitals (i.e., the number of digits).
    '''
    # call python method "bin" and remove the prefix "0b"
    binstr = bin(int(n)).replace("0b", "")
    # convert the string to a list using the edrixs convention (see: https://nsls-ii.github.io/edrixs/reference/fock_basis.html#edrixs.fock_basis.get_fock_full_N)
    binlen = len(binstr)
    binlist = []
    for i in range(norb):
        if i < binlen:
            binlist.append(int(binstr[-1-i]))
        else:
            binlist.append(0)
    return binlist

def get_fortran_fock_i(norb):
    '''
    read fock basis from file "fock_i.in" produced by the Fortran solver, and convert to the binary form.
    '''
    dec_arr = np.loadtxt('files/fock_i.in', dtype='int')[1:]
    fock_i = []
    for dec in dec_arr:
        fock_i.append(decimalToBinary(dec, norb))
    return fock_i

indices, e = np.loadtxt('files/eigvals.dat', unpack=True)

# note different order scipy.linalg.eigh
v = np.array([get_fortran_eigvec(eigvec_num)
              for eigvec_num in range(1, 190+1)])



norb_d = 10

basis = np.array(get_fortran_fock_i(190))
num_d_electrons = basis[:, :norb_d].sum(1)

alphas = np.sum(np.abs(v[:, num_d_electrons==8])**2, axis=1)
betas = np.sum(np.abs(v[:, num_d_electrons==9])**2, axis=1)
gammas = np.sum(np.abs(v[:, num_d_electrons==10])**2, axis=1)

fig, ax = plt.subplots()

ax.plot(e, alphas, label=r'$\alpha$ $d^8L^{10}$')
ax.plot(e, betas, label=r'$\beta$ $d^9L^{9}$')
ax.plot(e, gammas, label=r'$\gamma$ $d^{10}L^{8}$')

ax.set_xlabel('Energy (eV)')
ax.set_ylabel('Population')
ax.set_title('NiO')
ax.legend()
plt.show()