Comparing relax calculation with old and new Fermi energy¶

def query(group_label_new, group_label_old) :
    from aiida import orm, load_profile
    load_profile('sph')

    qb_new = orm.QueryBuilder().append(
        orm.Group        , tag='group_new', filters={'label': group_label_new}).append(
        orm.WorkChainNode, tag='new_relax_workchains' , with_group='group_new'                  , project='*', filters={'attributes.process_label': 'PwRelaxWorkChain', 'attributes.exit_status': 0}).append(
        orm.WorkChainNode, tag='new_base_workchains' , with_incoming='new_relax_workchains', project='*', filters={'attributes.process_label': 'PwBaseWorkChain'}, edge_filters={'label': 'final_scf'}).append(
        orm.BandsData    , tag='new_bands'           , with_incoming='new_base_workchains' , project='*').append(
        orm.Dict         , tag='new_out_parameters'  , with_incoming='new_base_workchains' , project='*', edge_filters={'label': 'output_parameters'}).append(
        orm.StructureData, tag='structures'          , with_outgoing='new_relax_workchains' , project='*')

    # qb_new.limit(10)

    print('Total number of entries in qb_new:', qb_new.count() )

    qb_old = orm.QueryBuilder().append(
        orm.Group        , tag='group_old', filters={'label': group_label_old}).append(
        orm.WorkChainNode, tag='new_relax_workchains' , with_group='group_old'                  , project='*', filters={'attributes.process_label': 'PwRelaxWorkChain', 'attributes.exit_status': 0}).append(
        orm.WorkChainNode, tag='new_base_workchains' , with_incoming='new_relax_workchains', project='*', filters={'attributes.process_label': 'PwBaseWorkChain'}, edge_filters={'label': 'final_scf'}).append(
        orm.BandsData    , tag='new_bands'           , with_incoming='new_base_workchains' , project='*').append(
        orm.Dict         , tag='new_out_parameters'  , with_incoming='new_base_workchains' , project='*', edge_filters={'label': 'output_parameters'}).append(
        orm.StructureData, tag='structures'          , with_outgoing='new_base_workchains' , project='*' )
      
    # qb_old.limit(1000)

    print('Total number of entries in qb_old:', qb_old.count() )

    qb_list = []

    qb_old_all = qb_old.all()

    for entry_new in qb_new.all() :
        for entry_old in qb_old_all :
            if entry_new[4].uuid == entry_old[4].uuid :
                qb_list.append( [entry_new, entry_old]  )
    
    return qb_new, qb_old, qb_list

def collect_data_HT(qb_list, json_file) :
    from aiida_quantumespresso_epfl.common.utils import get_bandgap
    import json

    # qb.limit(10)
    # json_file = 'analysis_selected_structure.json'

    line0 = '                         |----- differences (eV) ------|   |-- rel-errors (%) --|  (ev / angs)\n'
    line1 = 'count formula             efermi    ind_bgap  dir_bgap       ind_bgap  dir_bgap     new_force      old_pwbase                            new_pwbase'
    # print(line0 + line1)

    data = {}

    for counter, entry in enumerate( qb_list ) :

        qb_new, qb_old = entry
  
        new_pwrelax, new_pwbase, new_bands, new_out_parameters, old_structure = qb_new
        old_pwrelax, old_pwbase, old_bands, old_out_parameters, old_structure = qb_old

        # Getting the band gaps
        new_ind_bandgap, new_dir_bandgap = get_bandgap(new_bands, new_out_parameters)
        old_ind_bandgap, old_dir_bandgap = get_bandgap(old_bands, old_out_parameters)

        # Getting the Fermi energies
        new_fermi_energy = new_out_parameters.attributes['fermi_energy'] 
        old_fermi_energy = old_out_parameters.attributes['fermi_energy'] 

        # Getting the Fermi energies
        new_volume = new_out_parameters.attributes['volume'] 
        old_volume = old_out_parameters.attributes['volume'] 

        # Getting the forces
        new_force = new_out_parameters.attributes['total_force'] 
        old_force = old_out_parameters.attributes['total_force'] 

        # Getting formula
        formula= old_structure.get_formula()

        # Getting number of steps in the new calculation
        # num_steps = new_pwrelax.outputs.output_trajectory.get_array('steps')
        # print('num_steps', num_steps)

        total_steps = 0
        for wc in new_pwrelax.called:
            td = wc.outputs.output_trajectory
            total_steps += len(td.get_array('steps'))
        # print('total_steps', total_steps)

        # Computing the differences between new and old calculations
        efermi_diff      = new_fermi_energy - old_fermi_energy
        volume_diff      = new_volume - old_volume
        ind_bandgap_diff = new_ind_bandgap - old_ind_bandgap
        dir_bandgap_diff = new_dir_bandgap - old_dir_bandgap

        # Computing the relative errors between the new and old values
        ind_relative_error    = 100.*ind_bandgap_diff/old_ind_bandgap
        dir_relative_error    = 100.*dir_bandgap_diff/old_dir_bandgap
        volume_relative_error = 100.*volume_diff     /old_volume

        line = '{:5} {:17} {:10.5f}{:10.5f}{:10.5f}  {:10.2f}{:10.2f}      {:10.5f}       {:38}{:38}'.format(
                 counter+1, formula, efermi_diff, ind_bandgap_diff, dir_bandgap_diff, ind_relative_error, dir_relative_error, new_force, old_pwbase.uuid, new_pwbase.uuid)

    #     print(line)

        data[counter] = {
            'formula'               : formula ,
            'efermi_diff'           : efermi_diff ,
            'old_volume'            : old_volume ,
            'new_volume'            : new_volume ,
            'volume_diff'           : volume_diff ,
            'volume_relative_error' : volume_relative_error ,
            'total_steps'           : total_steps ,
            'old_ind_bandgap'       : old_ind_bandgap ,
            'new_ind_bandgap'       : new_ind_bandgap ,
            'ind_bandgap_diff'      : ind_bandgap_diff ,
            'dir_bandgap_diff'      : dir_bandgap_diff ,
            'ind_relative_error'    : ind_relative_error ,
            'dir_relative_error'    : dir_relative_error ,
            'new_force'             : new_force ,
            'old_pwbase_uuid'       : old_pwbase.uuid ,
            'new_pwbase_uuid'       : new_pwbase.uuid ,
        }

    with open( json_file, 'w' ) as file :
        result_string = json.dumps(data, indent=2)
        file.write(result_string)
    

json_file = 'analysis_selected_structure_HT3.json'
collect_data_HT(qb_list, json_file)

def sort(filter_key, toprint, toplot, json_file, value_key = None, zoomin=''):
    
    if value_key == None :
        value_key = filter_key

    data = load_json(json_file)
    
    line0 = '                         |----- differences (eV) ------|   |-- rel-errors (%) --|  (ev / angs)\n'
    line1 = 'count formula             efermi    ind_bgap  dir_bgap       ind_bgap  dir_bgap     new_force      old_pwbase                            new_pwbase'
    if toprint : print(line0 + line1)

    # filter_key = 'ind_relative_error'
    sorted_data = sorted(data.items(), key=lambda x: abs(x[1][filter_key]), reverse=True)    

    filter_column = []
    xs = []

    for counter, entry in enumerate(sorted_data) :
        entry_dict = entry[1]

        formula            = entry_dict['formula'           ]
        efermi_diff        = entry_dict['efermi_diff'       ]
        try : volume_diff           = entry_dict['volume_diff'          ]
        except : pass
        try : volume_relative_error = entry_dict['volume_relative_error']
        except : pass
        try : total_steps = entry_dict['total_steps']
        except : pass
        ind_bandgap_diff   = entry_dict['ind_bandgap_diff'  ]
        dir_bandgap_diff   = entry_dict['dir_bandgap_diff'  ]
        ind_relative_error = entry_dict['ind_relative_error']
        dir_relative_error = entry_dict['dir_relative_error']
        new_force          = entry_dict['new_force'         ]
        old_pwbase_uuid    = entry_dict['old_pwbase_uuid'   ]
        new_pwbase_uuid    = entry_dict['new_pwbase_uuid'   ]
        old_ind_bandgap    = entry_dict['old_ind_bandgap'   ]

        line = '{:5} {:17} {:10.5f}{:10.5f}{:10.5f}  {:10.2f}{:10.2f}      {:10.5f}       {:38}{:38}'.format(
                 counter+1, formula, efermi_diff, ind_bandgap_diff, dir_bandgap_diff, ind_relative_error, dir_relative_error, new_force, old_pwbase_uuid, new_pwbase_uuid)

        if toprint : print(line)

        filter_column.append(entry_dict[value_key])
        
        xs.append(old_ind_bandgap)
    
    if toplot : 
        # new_plot(filter_column,value_key,xs,zoomin) 
        plot(filter_column,value_key,filter_key)

import os
import matplotlib.pyplot as plt
import numpy as np
import json

def load_json(json_file) :

    #Testing if inputted 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
                collect_data = json.load(file)
                
            return collect_data
        #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.")
        
        return
    else : #If inputed file does NOT exist, terminate.
        exit( "Inputted file: '" + json_file + "' does NOT exist.")
        
        return


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

params3 = {'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(filter_column, value_key,filter_key='') :

    plt.rcParams.update(params3)

    size = len(filter_column)
    factor = 100./float(size)
    # xs = [i*factor for i in range(size)]
    xs = [i for i in range(size)]
    
    # plt.plot( xs, size*[0], color='gray', linestyle='--', linewidth=1)
    # plt.plot( xs, size*[0.001], color='red', linestyle='-', linewidth=1)
    plt.plot( xs, size*[0], color='black', linestyle='-', linewidth=1)
    # plt.plot( xs, size*[0.3], color='black', linestyle='-', linewidth=1, label='0.3')
    # plt.plot( xs, size*[0.5], color='black', linestyle='-', linewidth=1, label='0.5')
    # plt.plot( xs, size*[  1], color='black', linestyle='-', linewidth=1, label='1.0')
    
    # plt.bar( xs, filter_column,width=100./len(xs), label='value', color='blue')
    plt.bar( xs, filter_column,width=0.75, label='value', color='blue')
# width=0.8, bottom=None, *, align='center', data=None
    # if value_key not in  ['ind_bandgap_diff', 'volume_relative_error', 'volume_diff'] : 
    #     plt.plot( xs, np.abs(np.array(filter_column)), color='red', label='abs value', linewidth=2 )

    # plt.xlim(-1, 100)
    # plt.xlim(-1, 60)
    plt.xlim(-1, 285)
    if value_key == filter_key and value_key != 'old_ind_bandgap' :
        plt.xlim(-1, 50)
    plt.xlabel("Structure index")

    if value_key == 'old_ind_bandgap' :
        label = 'Old bandgap (eV)' 
    elif value_key == 'new_ind_bandgap' :
        label = 'New bandgap (eV)' 
    elif value_key == 'ind_bandgap_diff' :
        label = 'Bandgap difference (eV)'
        plt.ylim(-0.20, 0.25)
    elif value_key == 'ind_relative_error' :
        label = 'Bandgap relative error (%)'
        # plt.ylim(-5, 25)
        # plt.ylim(-2, 10)
        # plt.xlim(0, 30)
    elif value_key == 'new_force' :
        label = 'Force (eV/Å) '
    elif value_key == 'volume_diff' :
        label = 'Volume difference (Å^3) '
        # plt.ylim(-3, 0.20)
    elif value_key == 'volume_relative_error' :
        label = 'Volume relative error (%)'
        # plt.ylim(-1, 5)
        # plt.ylim(-1, 4)
        plt.ylim(-1, 5)
        # plt.xlim(0, 30)
    else :
        label = value_key

    # plt.yscale('log')
    print('value_key', value_key, label )

    plt.ylabel(label)
    # plt.legend()
    plt.tight_layout()
    print(value_key,filter_key)
    plt.savefig(value_key+'_sorted_by_'+filter_key+'_v2.png')
    plt.show()
    # if value_key == 'ind_relative_error' :
    #     plt.savefig(value_key+'.png')
    # elif value_key == 'volume_relative_error' :
    #     plt.savefig(value_key+'.png')

def new_plot(filter_column, filter_key, xs, zoomin='') :

    plt.rcParams.update(params3)

    size = len(filter_column)
    print("size", size)
    # factor = 100./float(size)
    # xs = [i*factor for i in range(size)]
    
    plt.plot( xs, size*[0], color='gray', linestyle='-', linewidth=1)
    
    # plt.bar( xs, filter_column,width=max(xs)/len(xs), label='value', color='blue')
    plt.plot( xs, filter_column, 'o', markersize=4, label='value', color='blue')
# width=0.8, bottom=None, *, align='center', data=None
    # plt.plot( xs, np.abs(np.array(filter_column)), color='red', label='abs value', linewidth=2 )
    
    # plt.yscale('log')
    
    # plt.xlim(-1, 100)
    # plt.xlim(0, 100)
    plt.xlabel("Bandgap (eV)")

    if filter_key == 'old_ind_bandgap' :
        label = 'Bandgap with old protocol (eV)' 
    if filter_key == 'new_ind_bandgap' :
        label = 'Bandgap with new protocol (eV)' 
    if filter_key == 'ind_bandgap_diff' :
        label = 'Bandgap difference (eV)'
        # plt.ylim(-5, 20)
    elif filter_key == 'ind_relative_error' :
        label = 'Bandgap relative error (%)'
        # plt.ylim(-5, 15)
        # plt.ylim(-2, 10)
        if zoomin == '_zoomin' : plt.ylim(-5, 15)
        plt.xlim(min(xs), max(xs))
    elif filter_key == 'new_force' :
        label = 'Force (eV/Å) '
    elif filter_key == 'volume_diff' :
        label = 'Volume difference (Å^3) '
        plt.ylim(-3, 20)
    elif filter_key == 'volume_relative_error' :
        label = 'Volume relative error (%)'
        # plt.ylim(-1, 5)
        # plt.ylim(-1, 4)
        # plt.xlim(0, 30)
        if zoomin == '_zoomin' : plt.ylim(-1, 4)
        plt.xlim(min(xs), max(xs))
    else :
        label = filter_key

    # plt.yscale('log')
    print('filter_key', filter_key, label )

    plt.ylabel(label)
    # plt.legend()
    plt.tight_layout()
    # # plt.savefig('force_kpoint.png')
    # if filter_key == 'ind_relative_error' :
    #     plt.savefig(filter_key+'_bandgap'+zoomin+'.png')
    # elif filter_key == 'volume_relative_error' :
    #     plt.savefig(filter_key+'_bandgap'+zoomin+'.png')
    plt.show()