Tutorial for Bogoliubov-de Gennes calculations with AiiDA-KKR¶

%load_ext autoreload
%autoreload 2
%matplotlib inline

from aiida import load_profile
profile = load_profile()

import numpy as np
import matplotlib.pyplot as plt
from aiida.orm import load_node, Dict, StructureData, Code
from aiida.engine import submit
from masci_tools.io.common_functions import get_aBohr2Ang
from masci_tools.io.common_functions import search_string
from aiida_kkr import __version__ as aiida_kkr_version
from aiida_kkr.tools import plot_kkr, kkrparams
from aiida_kkr.workflows import kkr_scf_wc, kkr_dos_wc

print('Version of AiiDA-KKR:', aiida_kkr_version, '(needs to be from BdG branch!)')


def load_or_submit(builder, uuid=None):
    """Utility to submit a calculation or load it it the uuid is given"""
    if uuid is None:
        calc = submit(builder)
        print('submitted', calc)
    else:
        calc = load_node(uuid)
    return calc


def get_delta(node, verbose=False):
    """Plot convergence of BdG pairing potential"""
    from aiida_kkr.tools.plot_kkr import get_rms_kkrcalc_from_remote
    try:
        _, _, _, _, res = get_rms_kkrcalc_from_remote(node, **{'l-independent delta:': None})
        d = np.array(res['l-independent delta:'])
    except:# ValueError:
        from masci_tools.io.common_functions import search_string
        d = []
        if node.is_finished:
            with node.outputs.retrieved.open('out_kkr') as f:
                txt = f.readlines()
                itmp = search_string('l-independent delta:', txt)
                while itmp>=0:
                    d.append(float(txt.pop(itmp).split()[-1]))
                    itmp = search_string('l-independent delta:', txt)
                d = np.array(d)
    if verbose:
        print('Delta=', d)
    if len(d)==0:
        print('no Delta values')
        d = []
    return d


# code that uses AMD nodes and has the BdG functionality
kkr_BdG = Code.get_from_string('kkrhost_BdG_AMD@iffslurm')
options_AMD = {
    'withmpi': True,
    'resources': {'num_machines': 1, 'tot_num_mpiprocs': 32},
    'queue_name': 'th1-2020-32', # select AMD nodes that match the AMD codes (see above)
    'max_wallclock_seconds': 3600*4 # 4 hours max runtime
}

if not Nb_scf.is_finished_ok:
    raise ValueError('Wait for scf to finish!')
    
# get the parameters that were used in the scf step
para_kkr = kkrparams(**{k:v for k,v in Nb_scf.inputs.calc_parameters.get_dict().items() if v is not None})

# update parameters to use a semi-circle contour
para_kkr.set_multiple_values(
    use_semi_circle_contour=True, # use semi-circle energy contour
    IM_E_CIRC_MIN= 5e-5, # minimal energy for semi-circle contour
    NPT1=32, # number of points on semi-circle contour
    MAX_NUM_KMESH=4, # limit the number of k-meshes to avoid undersampling in the BZ integration
    DISABLE_CHARGE_NEUTRALITY=True, # fix the Fermi level (otherwise the Fermi energy can jump a lot, this make the calculation stable)
    BZDIVIDE= [100, 100, 100], # larger k-mesh is needed since last energy points is close to real axis (see IM_E_CIRC_MIN)
    RCLUSTZ= 3.5, # increase cluster radius (larger screening cluster is necessary because of last energy point)
    NSTEPS= 200,
    IMIX= 4, 
    # this is for the use of the Chebychev solver (required to run BdG mode)
    RUNOPT=['NEWSOSOL'],
    R_LOG= 0.6,
    NPAN_EQ= 7,
    NPAN_LOG= 18,
    NCHEB= 12,
    DECOUPLE_SPIN_CHEBY=True, # use the Chebychev solver without SOC to speed everything up
)

# create a process builder for the KKR calculation and submit
builder = kkr_BdG.get_builder()
builder.metadata.label = Nb_scf.label + '_semi_circ'
builder.metadata.options = options_AMD
builder.parameters = Dict(dict=para_kkr)
builder.parent_folder = Nb_scf.outputs.last_RemoteData
Nb_semi_circ = load_or_submit(builder, '5a57de5b-02d3-4aaf-9875-25dafab076fb')
# plot_kkr(Nb_semi_circ, silent=True)
    
print('Finished with semi-circle contour?', Nb_semi_circ.is_finished_ok)

with Nb_semi_circ.outputs.retrieved.open('output.000.txt') as _f:
    txt = _f.readlines()
    
ec = []

# read RMS values for pairing potential
iline = search_string('ENERGY =', txt)
while (iline>=0):
    l = txt.pop(iline).split()
    ec.append([l[6], l[7], l[-1]])
    iline = search_string('ENERGY =', txt)
    
ec = np.array(ec, dtype=float)

plt.figure(figsize=(12,6))
plt.scatter(ec[:,0], ec[:,1], c=ec[:,2], )
plt.xlabel('Re(E) (Ry)')
plt.ylabel('Im(E) (Ry)')
cb = plt.colorbar()
plt.title('Energy contour')
plt.ylim(0)
plt.show()

if not Nb_semi_circ.is_finished_ok:
    raise ValueError('Wait for semi-circle calculation to finish!')

# activate BdG mode and change settings for one-shot calculation
# this reuses the settings from the semi-circle calculation (i.e. we reuse the same energy contour etc.)
para_BdG = kkrparams(**{k:v for k,v in Nb_semi_circ.inputs.parameters.get_dict().items() if v is not None})
para_BdG.set_multiple_values(
    # run 1 iteration for the one-shot calcualtion
    NSTEPS= 1,

    # activate BdG mode
    use_BdG=True,
    
    # 'big' starting value (improves convergence)
    Delta_BdG=5e-4,
    
    # use e/h symmetry (uses only half of the mirrored contour, makes everything faster)
    use_e_symm_BdG= True,
    
    # same Delta on all atoms (useful if non-superconducting atoms are in the structure)
    at_scale_BdG = [1.0 for i in range(len(struc.sites))],
)

# create calculation and run it
builder = kkr_BdG.get_builder() # reuse BdG code on AMD node
builder.parameters = Dict(dict=para_BdG)
builder.metadata.label = Nb_semi_circ.label + f'_BdG_start'
builder.metadata.options = options_AMD
# start from scf on the semi-circle contour
builder.parent_folder = Nb_semi_circ.outputs.remote_folder
BdG_scf_start = load_or_submit(builder, 'c4a2f6f8-b863-4f86-ba0c-da7fb8c51c70')

print('Finished with BdG initialization?', BdG_scf_start.is_finished_ok)

if not BdG_scf_start.is_finished_ok:
    raise ValueError('Wait for BdG init to finish!')

# update parmeters for self-consistency of BdG mode
para_BdG_scf = kkrparams(**{k:v for k,v in BdG_scf_start.inputs.parameters.get_dict().items() if v is not None})
para_BdG_scf.set_multiple_values(
    # run many iteration until self-consistency
    NSTEPS= 200,
    
    # coupling constant, this needs to be scaled for the different systems (trial and error) to reproduce the correct gap size
    lambda_BdG= 0.025,
    
    # BdG mixing:
    mixfac_BdG = 0.3, # mixing factor for BdG Broyden mixing
    Ninit_Broyden_BdG = 5, # number of BdG simple mixing steps (should be between 5 and 50 for stability) 

    # control potential mixing
    IMIX = 4, # use Broyden mixing in the potential
    NSIMPLEMIXFIRST = 25, # for the first 25 iterations we just do simple mixing in the potential to first converge the pairing potential and then let the (electronic) potential follow
)

# create new calculation and run it
builder = BdG_scf_start.get_builder_restart() # reuse BdG code on AMD node
builder.parameters = Dict(dict=para_BdG_scf)
builder.metadata.label = Nb_semi_circ.label + f'_BdG_scf'
builder.metadata.options = options_AMD
# start from BdG init (i.e. this uses the output anomalous density from the previous step
builder.parent_folder = BdG_scf_start.outputs.remote_folder
BdG_scf = load_or_submit(builder, '41285105-9951-47a7-86e4-65c5df5448ed')

print('Finished with BdG scf?', BdG_scf.is_finished_ok)

# plot convergence
plot_kkr(BdG_scf, silent=True)

# plot the convergence of the pairing potential
delta = get_delta(BdG_scf)
plt.figure()
plt.plot(delta, 'o-')
plt.ylabel('$\Delta$ (a.u.)')
plt.xlabel('iteration')
# add lines to indicate where mixing changes
plt.axvline(para_BdG_scf.get_value('NSIMPLEMIXFIRST')+0.5, ls=':', color='red', label='start Broyden for potential')
plt.axvline(para_BdG_scf.get_value('<NINIT_BROYDEN_BDG>')+0.5, ls=':', color='blue', label='start Broyden for $\Delta$')
plt.legend(fontsize='large')
plt.show()

if not BdG_scf.is_finished_ok:
    raise ValueError('Wait for BdG scf to finish!')

dos_settings = {'nepts': 32*4, 'tempr': 100, 'emin': -3.0, 'emax': 1.0, 'kmesh': [100, 100, 100]}

builder = kkr_dos_wc.get_builder()
builder.kkr = kkr_BdG
builder.options = Dict(dict=options_AMD)
builder.remote_data = BdG_scf.outputs.remote_folder
builder.wf_parameters = Dict(dict=dos_settings)
builder.metadata.label = BdG_scf.label + '_DOS_overview'

BdG_DOS_overview = load_or_submit(builder, '1cf7a898-4f48-4061-ae2b-b15792d7601e')
plot_kkr(BdG_DOS_overview, silent=True, noshow=True)
plt.title('BdG-DOS in large range around Fermi level')
plt.xlim(-3,1)
plt.show()

if not BdG_scf.is_finished_ok:
    raise ValueError('Wait for BdG scf to finish!')

# needs small temperature and large k-mesh!
dos_settings = {'nepts': 32*4, 'tempr': 0.5, 'emin': -0.01, 'emax': 0.01, 'kmesh': [300, 300, 300], 'KPOIBZ': 600000} # KPOIBZ needs to be big enough to fit the irreducable BZ

builder = kkr_dos_wc.get_builder()
builder.kkr = kkr_BdG
builder.options = Dict(dict={
    'withmpi': True,
    'resources': {'num_machines': 1, 'tot_num_mpiprocs': 64},
    'queue_name': 'th1-2020-64',
    'max_wallclock_seconds': 14400}
)
builder.remote_data = BdG_scf.outputs.remote_folder
builder.wf_parameters = Dict(dict=dos_settings)
builder.metadata.label = BdG_scf.label + '_DOS_zoom'

BdG_DOS_zoom = load_or_submit(builder, '63658dd6-88d2-49ee-bdf0-a5005dfaf6df')
plot_kkr(BdG_DOS_zoom, silent=True, xscale=1000, noshow=True, ptitle='SC gap in BdG-DOS')
plt.xlabel('$E-E_{\mathrm{F}}$ (meV)')
plt.xlim(-10,10)
plt.show()