In [2]:

import numpy.f2py as f2py

fid = open("QE subs/efermig.f90")
source = fid.read()
fid.close()

compilation_out = f2py.compile(source, modulename = "efermig", verbose=False, source_fn="temp.f90")

if compilation_out == 0 :
    import efermig
    print(efermig.efermig_sub.__doc__)
else :
    print("COMPILATION ERROR!")

efermig,nelec_ef = efermig_sub(et,nelec,wk,degauss,ngauss,is,isk,[nbnd,nks])

Wrapper for ``efermig_sub``.

Parameters
----------
et : input rank-2 array('d') with bounds (nbnd,nks)
nelec : input float
wk : input rank-1 array('d') with bounds (nks)
degauss : input float
ngauss : input int
is : input int
isk : input rank-1 array('i') with bounds (nks)

Other Parameters
----------------
nbnd : input int, optional
    Default: shape(et,0)
nks : input int, optional
    Default: shape(et,1)

Returns
-------
efermig : float
nelec_ef : float

In [3]:

import matplotlib.pyplot as plt
import os
import glob

rytoev = np.float(13.605691930242388) #From QE

def Read_Two_Column_File(file_name):
    with open(file_name, 'r') as data:
        x = []
        y = []
        for line in data:
            p = line.split()
            if p[0] != "#" :
                x.append(float(p[0])*rytoev)
                y.append(float(p[1]))

    return x, y

def check_if_string_in_file(file_name, string_to_search):
    """ Check if any line in the file contains given string """
    # Open the file in read only mode
    with open(file_name, 'r') as read_obj:
        # Read all lines in the file one by one
        for line in read_obj:
            # For each line, check if line contains the string
            if string_to_search in line:
                return True
    return False

def Read_Four_Column_File(file_name):
    with open(file_name, 'r') as data:
        x = []
        y = []
        z = []
        w = []
        for line in data:
            p = line.split()
            if p[0] != "#" :
                x.append(float(p[0]))
                y.append(float(p[1]))
                z.append(float(p[2]))
                w.append(float(p[3]))

    return x, y, z, w

<ipython-input-3-2ba4d3685d6a>:5: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  rytoev = np.float(13.605691930242388) #From QE

In [4]:

group_label = 'aiida_3dd_final_scf'

query = QueryBuilder()

PwCalculation = CalculationFactory('quantumespresso.pw')

query.append(Group        , tag='group', filters={'label': group_label})
print("Number of Groups '" + group_label + "':", query.count())

query.append(PwCalculation, with_group='group', tag='pwcalculation', filters={'attributes.exit_status': 0}, project='*' )
print("Number of successiful PwCalculations:", query.count())

query.append(StructureData, with_outgoing='pwcalculation', tag='structures', project='*' )
print("Number of Structures:", query.count())

query.append(BandsData    , with_incoming='pwcalculation', tag='bands', project='*' )
print("Number of Bands:", query.count())

query.append(Dict         , with_incoming='pwcalculation', tag='output_parameters', project='*', edge_filters={'label': 'output_parameters'} )
print("Number of Output Parameters:", query.count())

query.append(Dict         , with_outgoing='pwcalculation', tag='input_parameters', project='*', edge_filters={'label': 'parameters'} )
print("Number of Input Parameters:", query.count())


limit = 100
query.limit(limit)

Number of Groups 'aiida_3dd_final_scf': 1
Number of successiful PwCalculations: 24842
Number of Structures: 24842
Number of Bands: 24842
Number of Output Parameters: 24842
Number of Input Parameters: 24842

Out[4]:

<aiida.orm.querybuilder.QueryBuilder at 0x7f1d211313d0>

In [4]:

def query_by_uuid(uuid):

    # uuid = '03fdf8a9-7ee4-4c85-9bf0-d91fca82a07f' # large efermi diff = -3.87166 
    # uuid = '00f863fb-355c-4100-87f2-3c2dd1786e44' # efermi diff = 1.75443, but olf efermi derivative = 4.0191e-07
    # uuid = '0064e053-14e2-4113-94a2-d83cad1b2f47' # small efermi diff = -0.26632

    group_label = 'aiida_3dd_final_scf'

    query = QueryBuilder()

    PwCalculation = CalculationFactory('quantumespresso.pw')

    query.append(Group        , tag='group', filters={'label': group_label})
    query.append(PwCalculation, with_group='group', tag='pwcalculation', filters={'uuid': uuid}, project='*' )
    query.append(StructureData, with_outgoing='pwcalculation', tag='structures', project='*' )
    query.append(BandsData    , with_incoming='pwcalculation', tag='bands', project='*' )
    query.append(Dict         , with_incoming='pwcalculation', tag='output_parameters', project='*', edge_filters={'label': 'output_parameters'} )
    query.append(Dict         , with_outgoing='pwcalculation', tag='input_parameters', project='*', edge_filters={'label': 'parameters'} )

    limit = 1
    query.limit(limit)
    
    return query

In [5]:

import re
from aiida_quantumespresso_epfl.common.utils import get_bandgap
import json

def compute_new_efermi(query, file_name) :
    
    toprint = False

#     file_name = "fermi_energy_3dd.dat"

    json_file = file_name[0:-4]+'.json'

    output_file = open(file_name, "w") 

    line1 = "                                                                     |---------shapes---------|"
    line2 = "  idx formula              efermi  nks nbnd nelec degau smaer Ngauss  bands   kpoints  weights   ef_diff(eV) nspin  derv-old-ef  derv-new-ef  dir_gap"

    output_file.write( line1 + '\n' + line2 + '\n' )
    if toprint: print( line1 + '\n' + line2 )

    collected_data = {}

    for struc_idx, entry in enumerate( query.iterall() ) :

        node_pwcalculation, node_structure, node_bands, node_output_parameter, node_input_parameter = entry
    #     print('\n', node_pwcalculation)

        ind_bandgap, dir_bandgap = get_bandgap(node_bands, node_output_parameter)

        pwuuid       = node_pwcalculation.uuid 
        formula      = node_structure.get_formula()
        fermi_energy = node_output_parameter.attributes['fermi_energy'       ] 
        nks          = node_output_parameter.attributes['number_of_k_points' ] 
        nbnd         = node_output_parameter.attributes['number_of_bands'    ] 
        nelec        = node_output_parameter.attributes['number_of_electrons'] 
        nspin        = node_output_parameter.attributes['number_of_spin_components']
        degauss      = node_input_parameter.attributes['SYSTEM']['degauss' ]
        smearing     = node_input_parameter.attributes['SYSTEM']['smearing']

        weights = np.array( node_bands.get_array("weights")  , dtype=np.float64)
        kpoints = np.array( node_bands.get_kpoints()         , dtype=np.float64)
        bands   = np.array( node_bands.get_bands()           , dtype=np.float64)/rytoev

        if nspin == 1 :
            isk_array = np.ones(nks, dtype=int)
            is_ = 0 
        elif nspin == 2 :
            nks_half = nks
            nks *= 2

            isk_array = np.ones(nks, dtype=int)
            isk_array[nks_half : nks] = 2 * np.ones(nks_half, dtype=int)

            weights = np.array( 2*list(weights) )

            shape = bands.shape

            bands = bands.reshape([2*shape[1], shape[2]])

            is_ = 0
        else :
            print("Unexpected number of spin. Stopping.")
            exit()

        bands   = np.transpose( np.array(bands, dtype=np.float64) )

        smearing = "methfessel"
#         smearing = "cold"
        
        if   smearing == "cold" :
            Ngauss = -1
        elif smearing == "gauss" :
            Ngauss =  0
        elif smearing == "methfessel":
            Ngauss =  1
                
        new_efermi, nelec_ef = efermig.efermig_sub(       bands, nelec,  weights, degauss, Ngauss,                      is_, isk_array, nbnd, nks)
        derv_old_efermi = rytoev*efermig.sumkg1(          bands,         weights, degauss, Ngauss, fermi_energy/rytoev, is_, isk_array, nbnd, nks)
        derv_new_efermi = rytoev*efermig.sumkg1(          bands,         weights, degauss, Ngauss,   new_efermi/rytoev, is_, isk_array, nbnd, nks)
#         new_efermi_MP, nelec_ef_MP = efermig.efermig_sub( bands, nelec,  weights, degauss,      1,                      is_, isk_array, nbnd, nks)
        
        new_efermi *= rytoev
#         new_efermi_MP *= rytoev
        
        efermi_diff = fermi_energy - new_efermi

        line = '{:5} {:17} {:9.6f}{:5}{:5}{:6.1f}{:6.2f}  {:7}{:2}    {:9} {:9}{:6}   {:8.5f}    {:2}   {:12.4e} {:12.4e} {:9.5f}  '.format(
             struc_idx, formula, fermi_energy, nks, nbnd, nelec, degauss, smearing, Ngauss,
             str(bands.shape), str(kpoints.shape), str(weights.shape), efermi_diff, nspin, derv_old_efermi, derv_new_efermi, dir_bandgap)

        output_file.write( line ) 
        if toprint: print( line, end='\n' ) 

        if abs(efermi_diff) > 0.0001 :

            pattern = 'Final Fermi energy from'

            with open("efermig_out.dat", "r") as file:
                for line in file:
                    if re.search(pattern, line):
                        output_file.write( str(node_pwcalculation) + line )
        else :
            output_file.write('\n')
            
            
        # Determining the Fermi energy determination method
        if check_if_string_in_file('./efermig_out.dat', 'Final Fermi energy from Newton' ) :
            
            if check_if_string_in_file('./efermig_out.dat', "Newton's minimization method" ) :
                ef_method = "Newton's minimization"
            elif check_if_string_in_file('./efermig_out.dat', "Newton's root finding method" ) :
                ef_method = "Newton's root-finding"
                
        else : #If Final Fermi energy is from bisection
            if check_if_string_in_file('./efermig_out.dat', 'Final Fermi energy from Bisection using Ngauss smearing' ) :
                ef_method = "Second bisecion"
            else :
                ef_method = "First bisecion"

        collected_data[struc_idx] = {
            'pwuuid'          : pwuuid,
            'formula'         : formula,
            'old_efermi'      : fermi_energy,
            'new_efermi'      : new_efermi,
            'efermi_diff'     : efermi_diff,
            'nks'             : nks,
            'nbnd'            : nbnd,
            'nelec'           : nelec,
            'degauss'         : degauss,
            'smearing'        : smearing,
            'Ngauss'          : Ngauss,
            'nspin'           : nspin,
            'derv_old_efermi' : derv_old_efermi,
            'derv_new_efermi' : derv_new_efermi,
            'dir_bandgap'     : dir_bandgap,
            'ind_bandgap'     : ind_bandgap,
            'ef_method'       : ef_method,
            'nelec_ef'        : nelec_ef
        }

    output_file.close()
#     print('Compute_new_efermi is done.')

    #Writing into the output file in json format
    with open( json_file, 'w' ) as file :
        result_string = json.dumps(collected_data, indent=2)

        file.write(result_string)
#         if toprint : print("Output written into file: ", json_file)
        
    return

In [6]:

import json
import os
import glob

def load_json(json_file) :

    #Testing if inputed file exists
    if os.path.isfile(json_file) :
        #Try to read and load the input data
        try: 
            with open( json_file, 'r') as file :
             # storing data from inputed file
             collected_data = json.load(file)

        #Handling errors
        except ValueError as e:
            exit("An exception occurred. Probably the json file is invalid. Existing.")
        except :
            exit("Something went wrong at reading the json file. Existing.")

    else : #If inputed file does NOT exist, terminate.
        exit( "Inputted file: '" + json_file + "' does NOT exist.")
        
    return collected_data

In [7]:

# smearing='cold'
def plot_new_efermi() :
    toprint = False
    
    fid = open("efermig_out.dat")
    if toprint : print(fid.read())

    plt.plot([-20,20],[0,0],color='gray', linestyle=':', linewidth=1)
        
    file = glob.glob('./num_electrons_gauss.dat' )
    x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
    plt.plot(x,y,label="N(\u03B5) - nelec   smearing = gauss")

    file = glob.glob('./num_electrons.dat' )
    x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
    plt.plot(x,1.*np.array(y),label='N(\u03B5) - nelec   smearing = '+smearing)

    plt.scatter(fermi_energy,0, marker='x', color='blue', s=100, label='old efermi', zorder=10)
           
#     file = glob.glob('./dev1_num_electrons.dat' )
#     x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
#     plt.plot(x,1.*np.array(y), label='dev1')

#     file = glob.glob('./dev2_num_electrons.dat' )
#     x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
#     plt.plot(x,0.001*np.array(y), label='dev2')
    
    from_Newton = check_if_string_in_file('./efermig_out.dat', 'Final Fermi energy from Newton' )
    
    if from_Newton : #If Final Fermi energy is from Newtow's method

        Newtons_minimization = check_if_string_in_file('./efermig_out.dat', "Newton's minimization method" )
        
        if Newtons_minimization :

            file = glob.glob('./minimum_newton_minimization.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x,y, marker='o', c=colors, s=50, cmap='gnuplot')
            plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="new efermi - Newton's minimization")

#             if smearing == 'methfessel' :
            file = glob.glob('./num_electrons.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            plt.plot(x,np.array(y)**2,label='[N(\u03B5) - nelec]^2 smearing = '+smearing)
                
        else : #If Final Fermi energy is from root-finding method
#             file = glob.glob('./num_electrons.dat' )
#             x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
#             plt.plot(x,abs(np.array(y)), color='gray', label='abs(N(\u03B5) - nelec)   smearing = '+smearing)

            file = glob.glob('./minimum_newton_root.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x,y, marker='o', c=colors, s=50, cmap='gnuplot')
            plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="new efermi - Newton's root-finding")
    
    else : #If Final Fermi energy is from bisection
        
        if check_if_string_in_file('./efermig_out.dat', "Final Fermi energy from Bisection using Ngauss smearing:" ) : 
    
            file = glob.glob('./minimum_bisection.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="new efermi - 2nd Bisection")

        else:
            
            file = glob.glob('./minimum_bisection.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x,y, marker='o', c=colors, s=50, cmap='gnuplot')
            plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="new efermi - 1st Bisection")
        
    extra = 1.0

#     plt.xlim(fermi_energy-dir_bandgap-extra, fermi_energy+dir_bandgap+extra)
#     plt.ylim(-2*extra,2*extra)

    extra = max(ind_bandgap,extra)
    plt.xlim(new_efermi-extra, new_efermi+extra)
    plt.ylim(-0.5,0.5)

#     plt.xlim(0,7.5)
#     plt.ylim(-15,10)
    
    plt.xlabel("Chemical potential (eV)")
    plt.ylabel("Number of electrons")

    plt.legend()
    plt.show()

In [10]:

file_name = "fermi_energy_3dd_" + str(limit) + "struct_M-P.dat"
compute_new_efermi(query, file_name)

In [63]:

# smearing='cold'


params = {'font.size': 14,
            'font.family': 'serif',
            'legend.fontsize': 12,
            'mathtext.fontset': 'dejavuserif'}
            # 'legend.handlelength': 2}

# High resolution for the paper
params = {'font.size': 18,
            'font.family': 'serif',
            'legend.fontsize': 12,
            'mathtext.fontset': 'dejavuserif',
            'figure.figsize': (7, 4.5) ,
            'figure.dpi' : 200,
            'legend.handlelength': 2}

def plot_new_efermi2(xmin=None, xmax=None) :
    toprint = False

    plt.rcParams.update(params)
    
    efermig_out_file = '66990.dat'
    efermig_out_file = './efermig_out.dat'

    fid = open(efermig_out_file)
    if toprint : print(fid.read())

    plt.plot([-20,20],[0,0],color='gray', linestyle=':', linewidth=1)
    
    file = glob.glob('./num_electrons_gauss.dat' )
    x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
    plt.plot(x,y,label=r"$N(\mu)$ -- Gaussian")

    file = glob.glob('./num_electrons.dat' )
    x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
    plt.plot(x,y,label=r'$N(\mu)$ -- Methfessel-Paxton')

    if xmin and xmax:
        plt.plot(x,500*np.array(y)**2, '--',label=r'$[N(\mu)^2]\times 500$ -- Methf.-Pax.')

    # plt.scatter(fermi_energy,0, marker='x', color='blue', s=100, label='old efermi', zorder=10)
    
    from_Newton = check_if_string_in_file(efermig_out_file, 'Final Fermi energy from Newton' )

    if from_Newton : #If Final Fermi energy is from Newtow's method
        
        Newtons_minimization = check_if_string_in_file(efermig_out_file, "Newton's minimization method" )
        print(f'Newtons_minimization={Newtons_minimization}')
        
        if Newtons_minimization :

            file = glob.glob('./minimum_newton_minimization.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x,500*np.array(y), marker='o', c=colors, s=50, cmap='gnuplot', zorder=8)
            plt.scatter(x[-1],500*np.array(y[-1]), marker='o', color='black', s=50, label='correct solution', zorder=9)

            # file = glob.glob('./num_electrons.dat' )
            # x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            # plt.plot(x,np.array(y)**2,label='[N(\u03B5) - nelec]^2 smearing = '+smearing)
        
        else : #If Final Fermi energy is from root-finding method

            file = glob.glob('./minimum_newton_root.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x,y, marker='o', c=colors, s=50, cmap='gnuplot')
            plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="new efermi - Newton's root-finding")
    
    else : #If Final Fermi energy is from bisection
        if check_if_string_in_file(efermig_out_file, "Final Fermi energy from Bisection using Ngauss smearing:" ) : 
    
            file = glob.glob('./minimum_bisection.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="new efermi - 2nd Bisection")

        else:
            
            file = glob.glob('./minimum_bisection.dat' )
            x, y = Read_Two_Column_File(os.path.expanduser(file[0]))
            colors = [-i for i in range(len(x))]
            plt.scatter(x,y, marker='o', c=colors, s=50, cmap='gnuplot', zorder=9 )
            # plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="new efermi - 1st Bisection")
            plt.scatter(x[-1],y[-1], marker='o', color='black', s=50, label="correct solution", zorder=9)


#     plt.xlim(new_efermi-dir_bandgap, new_efermi+dir_bandgap)
#     plt.ylim(-0.5,0.5)
    
    plt.scatter([fermi_energy+0.01, 4.79], [0, 0], marker='x', color='red', s=100, label='incorrect solutions', zorder=10)
    
    if xmin and xmax: # for Cu2S4SnZn
        plt.text(7.95, 0.35, r'Cu$_2$S$_4$SnZn', size = 'large')  #Cu2S4SnZn
    else :
        plt.text(4.7, -0.35, r'P$_8$S$_{28}$Zr$_4$', size = 'large')

    extra = max(ind_bandgap, 1.)
    plt.xlim(new_efermi-extra, new_efermi+extra)
    plt.ylim(-0.5,0.5)

#     plt.xlim(-60,-10)
#     plt.ylim(-280,-220)

    if xmin and xmax:
        plt.xlim(xmin, xmax)

    
    plt.xlabel("Chemical potential (eV)")
    plt.ylabel(r"Number of electrons - $N_0$")

    plt.legend()
    plt.show()