def plot_dos_from_calc(calc, iatom, ispin=None, sum_atoms=False, average_atoms=False, orbital='total', xfactor=1, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False):
# get number of atoms
if calc.process_class == KkrCalculation or calc.process_class == KkrimpCalculation:
natom = calc.outputs.output_parameters.get_attribute('number_of_atoms_in_unit_cell')
else:
print('ERROR : node input must be a KkrCalculation or KkrimpCalculation node')
# get number of spins
nspin = calc.outputs.output_parameters.get_attribute('nspin')
# get number of atoms
natoms = calc.outputs.output_parameters.get_attribute('number_of_atoms_in_unit_cell')
if not lmdos:
orbital_map = {
'total': 0,'s': 1,'p': 2,'d': 3,'f': 4,'ns': -1
}
else:
orbital_map = {
's': 0,
'p1': 1, 'p2': 2, 'p3': 3,
'd1': 4,'d2': 5,'d3': 6,'d4': 7,'d5': 8,\
'f1': 9,'f2': 10,'f3': 11,'f4': 12,'f5': 13,'f6': 14,'f7': 15
}
if orbital in orbital_map:
iorb = orbital_map[orbital]
else:
if not lmdos:
print('ERROR: chose between total / s / p / d / f / ns')
else:
print('ERROR: Invalid orbital for lmdos')
# get dos_data node from the calc retrieved folder
if calc.process_class == KkrCalculation:
dos_host_dict = parse_dosfiles(calc.outputs.retrieved)
if calc.process_class == KkrimpCalculation:
try:
host_calc = calc.inputs.host_Greenfunction_folder.get_incoming(node_class=KkrCalculation).first().node
parent_calc = host_calc.inputs.parent_folder.get_incoming(node_class=KkrCalculation).first().node
ef = parent_calc.outputs.output_parameters.get_attribute('fermi_energy')
except:
#get fermi energy from parent impurity calculation
if 'parent_calc_folder' in calc.inputs:
host_folder = calc.inputs.parent_calc_folder
parent_calc = host_folder.get_incoming(node_class=CalcJob).all()[-1].node
out_potential_content = parent_calc.outputs.retrieved.get_object_content("out_potential")
else:
potential = calc.get_incoming(node_class=kkr_imp_sub_wc).first().node.inputs.host_imp_startpot
out_potential_content = potential.get_content()
out_potential_content = out_potential_content.split('\n')
ef = float(out_potential_content[3].split()[1])
dos_host_dict = parse_impdosfiles(calc.outputs.retrieved, natom, nspin, ef, lmdos)
dos_data = dos_host_dict['dos_data']
x_data = dos_data.get_x()[1][0]
y_data = dos_data.get_y()[iorb][1]
y_data_atoms = []
if nspin==1:
y_data_atom = np.zeros(len(x_data))
else:
y_data_atom = np.array([np.zeros(len(x_data)), np.zeros(len(x_data))])
plots = []
if type(iatom) == int:
iatom = [iatom]
elif iatom == 'all':
iatom = range(1,natoms+1,1)
elif iatom == 'vac':
vac_atoms = []
core_states = imp_calc.outputs.output_parameters.get_attribute('charge_core_states_per_atom')
for i, core_state in enumerate(core_states):
if core_state == 0:
vac_atoms.append(i+1)
iatom = vac_atoms
for atom in iatom:
if nspin ==1:
y_data_atom = y_data[atom-1]
elif nspin ==2:
y_data_atom = (y_data[(atom-1)*2], y_data[(atom-1)*2+1])
if ispin == 1:
y_data_atom = y_data_atom[0]
elif ispin == 2:
y_data_atom = y_data_atom[1]
elif ispin == None:
y_data_atom = y_data_atom[0] + y_data_atom[1]
else:
print('ERROR : Spin value should be 1 or 2')
if _abs == True:
y_data_atom=abs(y_data_atom)
y_data_atoms.append(y_data_atom)
if not return_data:
if not sum_atoms:
plot, = plt.plot(x_data*xfactor, y_data_atom*yfactor+yoffset)
plots.append(plot)
if sum_atoms:
y_data_sum = np.sum(y_data_atoms, axis=0)
if average_atoms:
plot, = plt.plot(x_data*xfactor, y_data_sum*yfactor/len(iatom)+yoffset)
else:
plot, = plt.plot(x_data*xfactor, y_data_sum*yfactor+yoffset)
if return_data:
return dos_host_dict, x_data*xfactor, y_data_atoms
else:
if len(iatom) == 1 or sum_atoms:
return plot
else:
return plots
def plot_ldos_iatom(DOS_calc, iatom, xfactor=1000, xlim=[-3,3], at_type=None, legend_show=False):
calc = DOS_calc.get_outgoing(node_class=kkr_imp_sub_wc).first().node.get_outgoing(node_class=KkrimpCalculation).first().node
# Plot and store the plot objects
plot1 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='total', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
plot2 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='s', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
plot3 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='p', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
plot4 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='d', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
plot5 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='f', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
plot6 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='ns', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
# Set colors for the plots
plot1.set_color('black')
plot1.set_linewidth(3)
ls = '-'
lw = 2
title_size = 24
label_size = 20
legend_size = 19
legend_loc = 2
plot2.set_color(f'C{0}')
plot2.set_linestyle(ls)
plot2.set_linewidth(lw)
plot3.set_color(f'C{1}')
plot3.set_linestyle(ls)
plot3.set_linewidth(lw)
plot4.set_color(f'C{2}')
plot4.set_linestyle(ls)
plot4.set_linewidth(lw)
plot5.set_color(f'C{3}')
plot5.set_linestyle(ls)
plot5.set_linewidth(lw)
plot6.set_color(f'C{4}')
plot6.set_linestyle(ls)
plot6.set_linewidth(lw)
plt.xlim(xlim)
plt.ylabel('DOS (1/eV)', fontsize=label_size)
plt.xlabel('Energy (meV)', fontsize=label_size)
plt.yticks(fontsize=label_size)
plt.xticks(fontsize=label_size)
if iatom ==1:
if legend_show:
plt.legend(labels=['Total', 's', 'p', 'd', 'f', 'ns'], fontsize=legend_size, loc=legend_loc)
plt.title('Gd', fontsize=title_size)
else:
if legend_show:
plt.legend(labels=['Total', 's', 'p', 'd', 'f', 'ns'], fontsize=legend_size, loc=legend_loc)
if at_type is None:
if iatom ==8:
plt.title('Nb1', fontsize=title_size)
if iatom ==9:
plt.title('Nb2', fontsize=title_size)
if iatom ==14:
plt.title('Nb3', fontsize=title_size)
else:
plt.title(f'{at_type}', fontsize=title_size)
def plot_d_dos_iatom(DOS_calc, iatom, xfactor=1000, xlim=[-3, 3], at_type=None, legend_show=False):
calc = DOS_calc.get_outgoing(node_class=kkr_imp_sub_wc).first().node.get_outgoing(node_class=KkrimpCalculation).first().node
# Plot and store the plot objects
plot1 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='total', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
# plot2 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='d', xfactor=1000, yfactor=1, _abs=False, yoffset=0., lmdos=False, return_data=False)
plot3 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='d1', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=True, return_data=False)
plot4 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='d2', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=True, return_data=False)
plot5 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='d3', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=True, return_data=False)
plot6 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='d4', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=True, return_data=False)
plot7 = plot_dos_from_calc(calc, iatom, ispin=None, orbital='d5', xfactor=xfactor, yfactor=1, _abs=False, yoffset=0., lmdos=True, return_data=False)
# Set colors for the plots
plot1.set_color('black')
plot1.set_linewidth(3)
# plot2.set_color(f'C{0}')
ls = '-'
lw = '2'
title_size = 24
label_size = 20
legend_size = 19
legend_loc = 2
plot3.set_color(f'C{0}')
plot3.set_linestyle(ls)
plot3.set_linewidth(lw)
plot4.set_color(f'C{1}')
plot4.set_linestyle(ls)
plot4.set_linewidth(lw)
plot5.set_color(f'C{2}')
plot5.set_linestyle(ls)
plot5.set_linewidth(lw)
plot6.set_color(f'C{3}')
plot6.set_linestyle(ls)
plot6.set_linewidth(lw)
plot7.set_color(f'C{4}')
plot7.set_linestyle(ls)
plot7.set_linewidth(lw)
plt.xlim(xlim)
plt.ylabel('DOS (1/eV)', fontsize=label_size)
plt.xlabel('Energy (meV)', fontsize=label_size)
plt.yticks(fontsize=label_size)
plt.xticks(fontsize=label_size)
if iatom==1:
if legend_show==True:
plt.legend(labels=['Total', '$d_{x^2-y^2}$', '$d_{xz}$', '$d_{z^2}$', '$d_{yz}$', '$d_{xy}$'], fontsize=legend_size, loc=legend_loc)
plt.title('Gd', fontsize=title_size)
else:
if legend_show==True:
plt.legend(labels=['Total', '$d_{x^2-y^2}$', '$d_{xz}$', '$d_{z^2}$', '$d_{yz}$', '$d_{xy}$'], fontsize=legend_size, loc=legend_loc)
if at_type is None:
if iatom==8:
plt.title('Nb1', fontsize=title_size)
if iatom==9:
plt.title('Nb2', fontsize=title_size)
if iatom==14:
plt.title('Nb3', fontsize=title_size)
else:
plt.title(f'{at_type}', fontsize=title_size)
def plot_gap_edges(gap=1.91, lim=3, color='black', alpha=0.4):
plt.axvspan(-lim, -gap, facecolor=color, alpha=alpha/8)
plt.axvspan(gap, lim, facecolor=color, alpha=alpha/8)
plt.axvline(-gap, color=color, linestyle='--', alpha=alpha)
plt.axvline(gap, color=color, linestyle='--', alpha=alpha)
def get_peak_energy(calc, iatom, iorb, lmdos=True):
if type(iorb) is not list:
iorb = [iorb]
iatom = iatom-1
if lmdos:
x_data = calc.outputs.dos_data_lm.get_x()[1][0]
y_data = calc.outputs.dos_data_lm.get_y()
else:
x_data = calc.outputs.dos_data.get_x()[1][0]
y_data = calc.outputs.dos_data.get_y()
x_peak = []
y_peak = []
for orb in iorb:
y_data_orb = y_data[orb][1]
y_data_orb_atom = y_data_orb[2*iatom] + y_data_orb[2*iatom+1]
y_max_orb_atom = np.max(y_data_orb_atom)
y_peak.append(y_max_orb_atom)
for i, val in enumerate(y_data_orb_atom):
if val == y_max_orb_atom:
x_peak.append(x_data[i]*1000)
return x_peak, y_peak
# Function to download the content of a remote file into memory
def get_remote_file_data(remote_data, relpath):
"""Retrieve the content of a remote file as an in-memory array."""
# Create a temporary file in-memory
with io.BytesIO() as temp_file:
with open('/tmp/temp_remote_file', 'wb') as temp_path:
remote_data.getfile(relpath, temp_path.name)
# Read and return the file data as a numpy array
with open(temp_path.name, 'r') as f:
return np.loadtxt(f)
# Function to store JIJ data into arrays
def get_J_all_energies(calc, pattern, reload_data=False, test=False):
remote_data = calc.outputs.remote_folder
# Get Gd atom indices
Gd_atoms = [int(k.split('=')[1]) + 1 for k, v in calc.inputs.settings_LDAU.items() if 'iatom' in k]
# print('Gd atoms:', Gd_atoms)
# Initialize data storage variables
J_tot_all, J_all = [], []
D_tot_all, D_all = [], []
Dx_tot_all, Dx_all = [], []
Dy_tot_all, Dy_all = [], []
Dz_tot_all, Dz_all = [], []
J_tot = D_tot = Dx_tot = Dy_tot = Dz_tot = 0
# Load data from node extra if exists and reload_data is False, otherwise read from files
if 'uuid_saved_Jij_data' in remote_data.extras and not reload_data:
print('Loading data from node extra')
data_node = load_node(remote_data.extras['uuid_saved_Jij_data'])
data_all = data_node.get_array('data')
else:
print('Loading data from files and storing in node extra')
# List all files in the remote directory
try:
all_files = remote_data.listdir()
except IOError as e:
print(f"Error while listing directory: {e}")
return
# Filter files matching the given pattern
Jij_files = sorted([f for f in all_files if re.match(pattern, f)])
if not Jij_files:
print("No files matching the pattern found.")
return
# Limit files for testing if `test` is set to True
if test:
Jij_files = Jij_files[:2]
data_all = [get_remote_file_data(remote_data, filename) for filename in Jij_files]
data_all = np.array(data_all)
# Save data to a new ArrayData node
data_node = ArrayData()
data_node.set_array('data', data_all)
data_node.store()
remote_data.set_extra('uuid_saved_Jij_data', data_node.uuid)
# Process each file and compute J and D values
for i in range(len(data_all)):
data = data_all[i]
for line in data:
if line[0] == Gd_atoms[0] and line[1] == Gd_atoms[1]:
J_matrix = line[4:].reshape(3, 3) * 1000 * get_Ry2eV() * -1 # Convert to meV and adjust sign
J = np.trace(J_matrix) / 3 # Average diagonal for isotropic J
Dx, Dy, Dz = J_matrix[1, 2], J_matrix[2, 0], J_matrix[0, 1]
D = np.sqrt(Dx**2 + Dy**2 + Dz**2)
# Accumulate values
J_tot += J
D_tot += D
Dx_tot += Dx
Dy_tot += Dy
Dz_tot += Dz
# Append values to arrays
J_all.append(J)
D_all.append(D)
Dx_all.append(Dx)
Dy_all.append(Dy)
Dz_all.append(Dz)
J_tot_all.append(J_tot)
D_tot_all.append(D_tot)
Dx_tot_all.append(Dx_tot)
Dy_tot_all.append(Dy_tot)
Dz_tot_all.append(Dz_tot)
# Print values for verification
# print(f'J = {J}, J_tot = {J_tot}')
# print(f'D = {D}, D_tot = {D_tot}')
break
return J_all, J_tot_all, D_all, D_tot_all, Dx_tot_all, Dy_tot_all, Dz_tot_all
# Function to plot DOS
def plot_dos(ax, DOS_calc, d_ef=0, ylabel_left=True, ylabel_right=True, no_ytick_labels_left=False, no_ytick_labels_right=False):
x_data = DOS_calc.outputs.dos_data_lm.get_x()[1][0]
y_data_all_i, c_all_i = [], []
for i, color in zip(range(4, 9), range(5)):
y_label = DOS_calc.outputs.dos_data_lm.get_y()[i][0]
y_data = DOS_calc.outputs.dos_data_lm.get_y()[i][1][0] + DOS_calc.outputs.dos_data_lm.get_y()[i][1][1] # Sum the two channels
y_data_all_i.append(y_data)
# Calculate cumulative charge
c_all = np.cumsum((x_data[::-1][1:] - x_data[::-1][:-1]) * y_data[::-1][:-1])
c_all_shifted = c_all - c_all[len(c_all) // 2] # Shift to have c=0 at ef
c_all_i.append(c_all_shifted)
#Total DOS
y_data_i_sum = np.sum(y_data_all_i, axis=0)
ax.plot(x_data, y_data_i_sum, label='d tot', color='black', linewidth=linewidth)
# Total cumulative charge
ax2 = plt.twinx(ax)
twin_axes.append(ax2) # Store the twin axis
c_all_i_sum = np.sum(c_all_i, axis=0)
c_all_i_sum = c_all_i_sum + d_ef # Shift to have c=c(ef) at ef
ax2.plot(x_data[::-1][:-1], c_all_i_sum, label='Total charge', color=f'C{3}', linestyle='-', linewidth=linewidth)
# ax2.tick_params(axis='y', colors=darker_C3)
if ylabel_left:
ax.set_ylabel(r'DOS (1/eV)', fontsize=label_fontsize)
ax.yaxis.set_label_coords(label_x_coord, 0.5) # Set fixed x-position for alignment
if no_ytick_labels_left:
ax.set_yticklabels([])
if ylabel_right:
ax2.set_ylabel(r'$\rho_d$ (e)', color=darker_C3, fontsize=label_fontsize)
ax2.yaxis.set_label_coords(1-label_x_coord, 0.5) # Set fixed x-position for alignment
if no_ytick_labels_right:
ax2.set_yticklabels([])
#ax.legend(labels=['$d_{x^2-y^2}$', '$d_{xz}$', '$d_{z^2}$', '$d_{yz}$', '$d_{xy}$', '$d_{tot}$', r'$\Delta\rho_{d}$'], loc=2, prop={'size': legend_size})
ax.axvline(0, color=f'black', linestyle='--', linewidth=linewidt_vline, alpha=alpha_vline)
# ax.axhline(0, color='black', linewidth=linewidt_hline, alpha=alpha_hline)
ax.set_xlim([-1, 1])
ax.set_ylim([0, 2.5])
ax2.set_ylim([0, 2.5])
# Function to plot J_tot and D_tot
def plot_j_tot(ax, J_tot, D_tot, J_offset, D_offset, label_y_right=False,
no_ytick_labels_left=False, no_ytick_labels_right=False):
x = np.linspace(-2, 2, 192)
J_shifted = np.array(J_tot) - (J_tot[95] + J_tot[96]) / 2 + J_offset
ax.plot(x, J_shifted, color='C0', linewidth=linewidth)
ax.axhline(0, color='C0', linewidth=linewidt_hline, alpha=alpha_hline)
ax.axvline(0, color='black', linestyle='--', linewidth=linewidt_vline, alpha=alpha_vline)
ax.set_ylim([-20, 40])
if no_ytick_labels_left:
ax.set_yticklabels([])
ax2 = ax.twinx()
twin_axes.append(ax2)
D_shifted = np.array(D_tot) - (D_tot[95] + D_tot[96]) / 2 + D_offset
ax2.plot(x, D_shifted, color='C1', linewidth=linewidth)
if label_y_right:
ax2.set_ylabel('D (meV)', color='C1', fontsize=label_fontsize)
ax2.axhline(0, color='C1', linewidth=linewidt_hline, alpha=alpha_hline)
ax2.set_ylim([-3.5, 1])
if no_ytick_labels_right:
ax2.set_yticklabels([])
def get_ldec_charges(calc):
# Access the 'retrieved' folder
retrieved_folder = calc.outputs.retrieved
# Open and read 'out_log.000.txt'
with retrieved_folder.open("out_log.000.txt", "r") as file:
lines = file.readlines()
# Find the last occurrence of 'l-dec'
last_occurrence_index = None
for i, line in enumerate(lines):
if 'l-dec' in line:
last_occurrence_index = i
# Initialize dictionary to store orbital information
orbital_data = {}
# If 'l-dec' was found, get the 12 lines following the last occurrence
if last_occurrence_index is not None:
selected_lines = lines[last_occurrence_index + 1:last_occurrence_index + 13]
# Define a regex pattern to capture orbital information with '=' sign
pattern = re.compile(r"(\w+)\s*=\s*([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)")
for line in selected_lines:
# Check if the line matches the pattern with '='
match = pattern.search(line)
if match:
orbital = match.group(1)
spin_dn = float(match.group(2))
spin_up = float(match.group(3))
m_spin = float(match.group(4))
# Store the data in the dictionary
orbital_data[orbital] = {
'spin_dn': spin_dn,
'spin_up': spin_up,
}
elif len(line.strip().split()) == 4: # Handle 'ns' and 'TOT' lines
parts = line.strip().split()
orbital = parts[0]
spin_dn = float(parts[1])
spin_up = float(parts[2])
m_spin = float(parts[3])
# Store the data in the dictionary
orbital_data[orbital] = {
'spin_dn': spin_dn,
'spin_up': spin_up,
}
return orbital_data
else:
print("No 'l-dec' found in the file.")
return None
def get_J_and_D(calc):
impurity1_output_node = impurity2_output_node = Dimer_short_wf.inputs.impurity1_output_node
retrieved = calc.outputs.retrieved
impurity_info = calc.inputs.impurity_info
calc_JiJ_data = parse_Jij(retrieved, impurity_info, impurity1_output_node, impurity2_output_node)['Jijdata'].get_array('JijData')
J = calc_JiJ_data[0][3]
D = calc_JiJ_data[0][4]
return J, D
