514 lines
25 KiB
Python
514 lines
25 KiB
Python
import numpy as np
|
|
from . import tip_types
|
|
from .GerischerMarkus import GM
|
|
import matplotlib.pyplot as plt
|
|
from scipy.interpolate import griddata
|
|
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
|
from pymatgen.io.vasp import Procar
|
|
from echem.core import constants
|
|
from echem.io_data import vasp
|
|
|
|
|
|
class kHET:
|
|
"""
|
|
This class calculates heterogeneous electron transfer rate constant with spatial resolution
|
|
"""
|
|
# TODO update tips types
|
|
|
|
AVAILABLE_TIPS_TYPES = ['oxygen', 'IrCl6', 'RuNH3_6', 'RuNH3_6_NNN_plane', 'RuNH3_6_perpendicular',
|
|
'oxygen_parallel_x', 'oxygen_parallel_y']
|
|
|
|
def __init__(self, working_folder=''):
|
|
if working_folder == '':
|
|
working_folder = '.'
|
|
# TODO think about better than get class objects for outcar and poscar info
|
|
self.outcar = vasp.Outcar.from_file(working_folder + '/OUTCAR')
|
|
self.poscar = vasp.Poscar.from_file(working_folder + '/POSCAR')
|
|
self.procar = Procar(working_folder + '/PROCAR')
|
|
self.working_folder = working_folder
|
|
self.path_to_data = working_folder + '/Saved_data'
|
|
self.wavecar = None
|
|
|
|
self.C_EDL = None
|
|
self.T = None
|
|
self.lambda_ = None
|
|
self.sheet_area = None
|
|
self.V_std = None
|
|
self.overpot = None
|
|
self.dE_Q = None
|
|
self.kb_array = None
|
|
self.E = None
|
|
|
|
def set_parameters(self, T, lambda_, overpot=0, V_std=None, C_EDL=None, dE_Q=None, linear_constant=26, threshold_value=1e-5, shifted_DOS='all'):
|
|
"""
|
|
:param C_EDL: float
|
|
Capacitance of electric double layer (microF/cm^2)
|
|
:param T: int, float
|
|
Temperature. It is used in computing Fermi function and distribution function of redox system states
|
|
:param lambda_: float
|
|
Reorganization energy in eV
|
|
:param V_std: float
|
|
Standart potential of the redox couple (Volts)
|
|
:param overpot: float
|
|
Overpotential (Volts). It shifts the electrode Fermi energy to -|e|*overpot
|
|
:return:
|
|
"""
|
|
#TODO check get_DOS parameters
|
|
self.E, DOS = self.outcar.get_DOS(zero_at_fermi=True, smearing='Gaussian', sigma=0.1, dE = 0.01)
|
|
self.T = T
|
|
self.lambda_ = lambda_
|
|
self.overpot = overpot
|
|
self.sheet_area = self._get_sheet_area() # Area of investigated surface(XY) in cm^2
|
|
self.linear_constant = linear_constant
|
|
if V_std == None or C_EDL == None:
|
|
if dE_Q == None:
|
|
raise ValueError ("Set either V_std and C_EDL or dE_Q parameter")
|
|
else:
|
|
self.dE_Q = dE_Q
|
|
else:
|
|
self.V_std = V_std
|
|
self.C_EDL = C_EDL
|
|
self.dE_Q = self._calculate_dE_Q()
|
|
if shifted_DOS == 'all':
|
|
self.are_frozen_states = False
|
|
self.kb_array = self._get_kb_array(threshold_value)
|
|
else:
|
|
self.are_frozen_states = True
|
|
self.kb_array = self._get_kb_array(threshold_value, shifted_DOS)
|
|
if self.kb_array == []:
|
|
print('Error! kb_array is empty. Try to decrease threshold_value')
|
|
|
|
def load_wavecar(self):
|
|
self.wavecar = vasp.Wavecar.from_file(self.working_folder + '/WAVECAR', self.kb_array)
|
|
|
|
def plot_distributions(self):
|
|
#TODO make it
|
|
pass
|
|
|
|
def _get_atom_localization(self, kpoint, band, target_atom_types):
|
|
for key in self.procar.data.keys():
|
|
procar_data = self.procar.data[key]
|
|
list_of_ions = []
|
|
atomnames = self.poscar.structure.species
|
|
orbital_names = self.procar.orbitals
|
|
for i, name in enumerate(atomnames):
|
|
if name in target_atom_types:
|
|
list_of_ions.append(i)
|
|
list_of_orbitals = [i for i in range(len(orbital_names))]
|
|
|
|
localization = 0
|
|
for ion in list_of_ions:
|
|
for orb in list_of_orbitals:
|
|
localization += procar_data[kpoint - 1][band - 1][ion][orb]
|
|
return localization
|
|
|
|
def _get_kb_array(self, threshold_value, shifted_DOS='all', threshold_for_localization_to_be_shifted=0.4):
|
|
if shifted_DOS=='all':
|
|
fermi_distribution = GM.fermi_func(self.E - self.dE_Q, self.T)
|
|
W_ox = GM.W_ox(self.E - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
E_satisfy_mask = W_ox * fermi_distribution > np.max(W_ox) * np.max(fermi_distribution) * threshold_value
|
|
list_of_E_indices = [i for i, E in enumerate(E_satisfy_mask) if E]
|
|
min_E_ind = min(list_of_E_indices)
|
|
max_E_ind = max(list_of_E_indices)
|
|
Erange = [self.E[min_E_ind], self.E[max_E_ind]]
|
|
|
|
kb_array = []
|
|
for band in range(1, self.outcar.nbands + 1):
|
|
for kpoint in range(1, self.outcar.nkpts + 1):
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
if energy >= Erange[0] and energy < Erange[1]:
|
|
kb_array.append([kpoint, band])
|
|
else:
|
|
fermi_distribution_shifted = GM.fermi_func(self.E - self.dE_Q, self.T)
|
|
W_ox_shifted = GM.W_ox(self.E - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
E_satisfy_mask_shifted = W_ox_shifted * fermi_distribution_shifted > np.max(W_ox_shifted) * np.max(
|
|
fermi_distribution_shifted) * threshold_value
|
|
list_of_E_indices_shifted = [i for i, E in enumerate(E_satisfy_mask_shifted) if E]
|
|
min_E_ind_shifted = min(list_of_E_indices_shifted)
|
|
max_E_ind_shifted = max(list_of_E_indices_shifted)
|
|
Erange_shifted = [self.E[min_E_ind_shifted], self.E[max_E_ind_shifted]]
|
|
|
|
fermi_distribution_frozen = GM.fermi_func(self.E - 0, self.T)
|
|
W_ox_frozen = GM.W_ox(self.E - 0 - self.overpot, self.T, self.lambda_)
|
|
E_satisfy_mask_frozen = W_ox_frozen * fermi_distribution_frozen > np.max(W_ox_frozen) * np.max(
|
|
fermi_distribution_frozen) * threshold_value
|
|
list_of_E_indices_frozen = [i for i, E in enumerate(E_satisfy_mask_frozen) if E]
|
|
min_E_ind_frozen = min(list_of_E_indices_frozen)
|
|
max_E_ind_frozen = max(list_of_E_indices_frozen)
|
|
Erange_frozen = [self.E[min_E_ind_frozen], self.E[max_E_ind_frozen]]
|
|
|
|
kb_array = []
|
|
self.frozen_mask = []
|
|
for band in range(1, self.outcar.nbands + 1):
|
|
for kpoint in range(1, self.outcar.nkpts + 1):
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
if energy >= Erange_frozen[0] and energy < Erange_frozen[1]:
|
|
if self._get_atom_localization(kpoint, band,
|
|
target_atom_types=shifted_DOS) < threshold_for_localization_to_be_shifted:
|
|
kb_array.append([kpoint, band])
|
|
self.frozen_mask.append(1)
|
|
if energy >= Erange_shifted[0] and energy < Erange_shifted[1]:
|
|
if self._get_atom_localization(kpoint, band,
|
|
target_atom_types=shifted_DOS) >= threshold_for_localization_to_be_shifted:
|
|
kb_array.append([kpoint, band])
|
|
self.frozen_mask.append(0)
|
|
return kb_array
|
|
|
|
def calculate_kHET_spatial(self, tip_type='s', z_pos=None, from_center=False, cutoff_in_Angstroms=5, all_z=False, dim='2D', shifted_separately=False):
|
|
xlen, ylen, zlen = np.shape(self.wavecar.wavefunctions[0])
|
|
|
|
if dim == '2D':
|
|
k_HET_ = np.zeros((xlen, ylen))
|
|
if z_pos == None:
|
|
raise ValueError('For dim=2D z_pos is obligatory parameter')
|
|
b1, b2, b3 = self.poscar.structure.lattice
|
|
bn1 = b1 / xlen
|
|
bn2 = b2 / ylen
|
|
bn3 = b3 / zlen
|
|
if b3[0] != 0.0 or b3[1] != 0.0:
|
|
print("WARNING! You z_vector is not perpendicular to XY plane, Check calculate_kHET_spatial_2D function")
|
|
if from_center:
|
|
z = int(zlen // 2 + z_pos // np.linalg.norm(bn3))
|
|
else:
|
|
z = int(z_pos // np.linalg.norm(bn3))
|
|
real_z_position = z_pos // np.linalg.norm(bn3) * np.linalg.norm(bn3)
|
|
print(f"Real z_pos of tip = {real_z_position}")
|
|
|
|
if any(tip_type == kw for kw in self.AVAILABLE_TIPS_TYPES):
|
|
cutoff = int(cutoff_in_Angstroms // np.linalg.norm(bn1))
|
|
if cutoff > xlen // 2 or cutoff > ylen // 2:
|
|
print("ERROR: Cutoff should be less than 1/2 of each dimension of the cell. "
|
|
"Try to reduce cutoff. Otherwise, result could be unpredictible")
|
|
acc_orbitals = self.generate_acceptor_orbitals(tip_type, (xlen, ylen, zlen), z_shift=z)
|
|
if all_z == True:
|
|
zmin = 0
|
|
zmax = zlen
|
|
elif z - cutoff >= 0 and z + cutoff <= zlen:
|
|
zmin = z - cutoff
|
|
zmax = z + cutoff
|
|
else:
|
|
print("Can't reduce integrating area in z dimention. "
|
|
"Will calculate overlapping for all z. You can ignore this message "
|
|
"if you don't care about efficiency")
|
|
zmin = 0
|
|
zmax = zlen
|
|
|
|
for i, kb in enumerate(self.wavecar.kb_array):
|
|
kpoint, band = kb[0], kb[1]
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
weight = self.outcar.weights[0][kpoint - 1]
|
|
f_fermi = GM.fermi_func(energy - self.dE_Q, self.T)
|
|
w_redox = GM.W_ox(energy - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
overlap_integrals_squared = self._get_overlap_integrals_squared(self.wavecar.wavefunctions[i], cutoff,
|
|
acc_orbitals, zmin, zmax)
|
|
# TODO check eq below
|
|
matrix_elements_squared = overlap_integrals_squared * self.linear_constant * np.linalg.norm(bn3) ** 3
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
|
|
elif tip_type == 's':
|
|
for i, kb in enumerate(self.wavecar.kb_array):
|
|
kpoint, band = kb[0], kb[1]
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
weight = self.outcar.weights[0][kpoint - 1]
|
|
f_fermi = GM.fermi_func(energy - self.dE_Q, self.T)
|
|
w_redox = GM.W_ox(energy - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
matrix_elements_squared = np.abs(self.wavecar.wavefunctions[i][:, :, z]) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
|
|
elif tip_type == 'pz':
|
|
for i, kb in enumerate(self.wavecar.kb_array):
|
|
kpoint, band = kb[0], kb[1]
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
weight = self.outcar.weights[0][kpoint - 1]
|
|
f_fermi = GM.fermi_func(energy - self.dE_Q, self.T)
|
|
w_redox = GM.W_ox(energy - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
wf_grad_z = np.gradient(self.wavecar.wavefunctions[i], axis=2)
|
|
matrix_elements_squared = np.abs(wf_grad_z[:, :, z]) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
else:
|
|
print(f"Try another tip type, for now {tip_type} is unavailiable")
|
|
|
|
elif dim == '3D':
|
|
k_HET_ = np.zeros((xlen, ylen, zlen))
|
|
|
|
if tip_type == 's':
|
|
if self.are_frozen_states:
|
|
if shifted_separately:
|
|
k_HET_shifted = np.zeros((xlen, ylen, zlen))
|
|
for i, kb in enumerate(self.wavecar.kb_array):
|
|
kpoint, band = kb[0], kb[1]
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
weight = self.outcar.weights[kpoint - 1]
|
|
if self.frozen_mask[i] == 1:
|
|
f_fermi = GM.fermi_func(energy, self.T)
|
|
w_redox = GM.W_ox(energy - self.overpot, self.T, self.lambda_)
|
|
matrix_elements_squared = np.abs(self.wavecar.wavefunctions[i]) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
else:
|
|
f_fermi = GM.fermi_func(energy - self.dE_Q, self.T)
|
|
w_redox = GM.W_ox(energy - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
matrix_elements_squared = np.abs(self.wavecar.wavefunctions[i]) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
try:
|
|
k_HET_shifted += matrix_elements_squared * f_fermi * w_redox * weight
|
|
except:
|
|
pass
|
|
else:
|
|
for i, kb in enumerate(self.wavecar.kb_array):
|
|
kpoint, band = kb[0], kb[1]
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
weight = self.outcar.weights[0][kpoint - 1]
|
|
f_fermi = GM.fermi_func(energy - self.dE_Q, self.T)
|
|
w_redox = GM.W_ox(energy - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
matrix_elements_squared = np.abs(self.wavecar.wavefunctions[i]) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
|
|
elif tip_type == 'pz':
|
|
if self.are_frozen_states:
|
|
if shifted_separately:
|
|
k_HET_shifted = np.zeros((xlen, ylen, zlen))
|
|
for i, kb in enumerate(self.wavecar.kb_array):
|
|
kpoint, band = kb[0], kb[1]
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
weight = self.outcar.weights[0][kpoint - 1]
|
|
if self.frozen_mask[i] == 1:
|
|
f_fermi = GM.fermi_func(energy, self.T)
|
|
w_redox = GM.W_ox(energy - self.overpot, self.T, self.lambda_)
|
|
wf_grad_z = np.gradient(self.wavecar.wavefunctions[i], axis=2)
|
|
matrix_elements_squared = np.abs(wf_grad_z) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
else:
|
|
f_fermi = GM.fermi_func(energy - self.dE_Q, self.T)
|
|
w_redox = GM.W_ox(energy - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
wf_grad_z = np.gradient(self.wavecar.wavefunctions[i], axis=2)
|
|
matrix_elements_squared = np.abs(wf_grad_z) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
try:
|
|
k_HET_shifted += matrix_elements_squared * f_fermi * w_redox * weight
|
|
except:
|
|
pass
|
|
else:
|
|
for i, kb in enumerate(self.wavecar.kb_array):
|
|
kpoint, band = kb[0], kb[1]
|
|
energy = self.outcar.eigenvalues[0][kpoint - 1][band - 1] - self.outcar.efermi
|
|
weight = self.outcar.weights[0][kpoint - 1]
|
|
f_fermi = GM.fermi_func(energy - self.dE_Q, self.T)
|
|
w_redox = GM.W_ox(energy - self.dE_Q - self.overpot, self.T, self.lambda_)
|
|
wf_grad_z = np.gradient(self.wavecar.wavefunctions[i], axis=2)
|
|
matrix_elements_squared = np.abs(wf_grad_z) ** 2
|
|
k_HET_ += matrix_elements_squared * f_fermi * w_redox * weight
|
|
|
|
else:
|
|
raise ValueError("dim should be 3D or 2D")
|
|
|
|
# TODO: check THIS below
|
|
#k_HET_ *= 2 * np.pi / constants.PLANCK_CONSTANT
|
|
if shifted_separately:
|
|
return k_HET_, k_HET_shifted
|
|
else:
|
|
return k_HET_
|
|
|
|
def plot_2D(self, func, show=True, save=False, filename='fig.png', method='linear', logscale=None):
|
|
"""
|
|
Function for plotting 2D images
|
|
"""
|
|
xlen, ylen = np.shape(func)
|
|
b1, b2, b3 = self.poscar.structure.lattice
|
|
bn1 = b1 / xlen
|
|
bn2 = b2 / ylen
|
|
|
|
R = []
|
|
for j in range(ylen):
|
|
for i in range(xlen):
|
|
R.append(i * bn1 + j * bn2)
|
|
X = np.array([x[0] for x in R])
|
|
Y = np.array([x[1] for x in R])
|
|
|
|
Z = func.transpose().flatten()
|
|
xi = np.linspace(X.min(), X.max(), 1000)
|
|
yi = np.linspace(Y.min(), Y.max(), 1000)
|
|
zi = griddata((X, Y), Z, (xi[None, :], yi[:, None]), method=method)
|
|
if logscale == None:
|
|
plt.contourf(xi, yi, zi, 500, cmap=plt.cm.rainbow)
|
|
else:
|
|
zi = np.log10(zi)
|
|
min, max, step = logscale
|
|
levels = np.arange(min, max, step)
|
|
plt.contourf(xi, yi, zi, 500, cmap=plt.cm.rainbow, levels=levels)
|
|
ax = plt.gca()
|
|
ax.set_aspect('equal')
|
|
divider = make_axes_locatable(ax)
|
|
cax = divider.append_axes("right", size="7%", pad=0.1)
|
|
cbar = plt.colorbar(cax=cax)
|
|
if show == True:
|
|
plt.show()
|
|
if save == True:
|
|
plt.savefig(filename, dpi=300)
|
|
plt.close()
|
|
|
|
def _calculate_dE_Q(self):
|
|
"""
|
|
This function calls GerischerMarcus module with GM class to calculate distribution of redox species
|
|
and Fermi-Dirac distribution according to Gerischer-Marcus formalism
|
|
:return:
|
|
"""
|
|
gm = GM(path_to_data=self.path_to_data)
|
|
gm.set_params(self.C_EDL, self.T, self.lambda_, self.sheet_area)
|
|
return gm.compute_distributions(self.V_std, overpot=self.overpot, add_info=True)[2]
|
|
|
|
def _get_sheet_area(self):
|
|
"""
|
|
Inner function to calculate sheet_area (XY plane) in cm^2
|
|
"""
|
|
b1, b2, b3 = self.poscar.structure.lattice
|
|
return np.linalg.norm(np.cross(b1, b2))*1e-16
|
|
|
|
def generate_acceptor_orbitals(self, orb_type, shape, z_shift=0, x_shift=0, y_shift=0):
|
|
"""
|
|
This is generator of acceptor orbitals using tip_types.py module
|
|
:return:
|
|
"""
|
|
bohr_radius = 0.529177 #TODO Check
|
|
b1, b2, b3 = self.poscar.structure.lattice
|
|
bn1 = b1 / shape[0]
|
|
bn2 = b2 / shape[1]
|
|
bn3 = b3 / shape[2]
|
|
basis = np.array([bn1, bn2, bn3])
|
|
transition_matrix = basis.transpose()
|
|
number_of_orbitals = len(tip_types.orbitals(orb_type))
|
|
acc_orbitals = []
|
|
for i in range(number_of_orbitals):
|
|
acc_orbitals.append(np.zeros(shape))
|
|
for i in range(shape[0]):
|
|
for j in range(shape[1]):
|
|
for k in range(shape[2]):
|
|
if i - x_shift >= shape[0] / 2:
|
|
x = i - shape[0] - x_shift
|
|
else:
|
|
x = i - x_shift
|
|
if j - y_shift >= shape[1] / 2:
|
|
y = j - shape[1] - y_shift
|
|
else:
|
|
y = j - y_shift
|
|
if k - z_shift >= shape[2] / 2:
|
|
z = k - shape[2] - z_shift
|
|
else:
|
|
z = k - z_shift
|
|
r = np.dot(transition_matrix, np.array([x, y, z])) / bohr_radius
|
|
for o, orbital in enumerate(acc_orbitals):
|
|
orbital[i][j][k] = tip_types.orbitals(orb_type)[o](r)
|
|
return acc_orbitals
|
|
|
|
def _get_overlap_integrals_squared(self, wf, cutoff, acc_orbitals, zmin, zmax):
|
|
xlen, ylen, zlen = np.shape(wf)
|
|
overlap_integrals_squared = np.zeros((xlen, ylen))
|
|
for i in range(xlen):
|
|
if i - cutoff < 0:
|
|
wf_rolled_x = np.roll(wf, cutoff, axis=0)
|
|
orb_rolled_x = []
|
|
for orbital in acc_orbitals:
|
|
orb_rolled_x.append(np.roll(orbital, i + cutoff, axis=0))
|
|
xmin = i
|
|
xmax = i + cutoff * 2
|
|
elif i - cutoff >= 0 and i + cutoff <= xlen:
|
|
wf_rolled_x = np.copy(wf)
|
|
orb_rolled_x = []
|
|
for orbital in acc_orbitals:
|
|
orb_rolled_x.append(np.roll(orbital, i, axis = 0))
|
|
xmin = i - cutoff
|
|
xmax = i + cutoff
|
|
elif i + cutoff > xlen:
|
|
wf_rolled_x = np.roll(wf, -cutoff, axis=0)
|
|
orb_rolled_x = []
|
|
for orbital in acc_orbitals:
|
|
orb_rolled_x.append(np.roll(orbital, i - cutoff, axis=0))
|
|
xmin = i - cutoff * 2
|
|
xmax = i
|
|
else:
|
|
print(f"ERROR: for i = {i} something with rolling arrays along x goes wrong")
|
|
for j in range(ylen):
|
|
if j - cutoff < 0:
|
|
wf_rolled = np.roll(wf_rolled_x, cutoff, axis=1)
|
|
orb_rolled = []
|
|
for orbital in orb_rolled_x:
|
|
orb_rolled.append(np.roll(orbital, j + cutoff, axis=1))
|
|
ymin = j
|
|
ymax = j + cutoff * 2
|
|
elif j - cutoff >= 0 and j + cutoff <= ylen:
|
|
wf_rolled = np.copy(wf_rolled_x)
|
|
orb_rolled = []
|
|
for orbital in orb_rolled_x:
|
|
orb_rolled.append(np.roll(orbital, j, axis=1))
|
|
ymin = j - cutoff
|
|
ymax = j + cutoff
|
|
elif j + cutoff > ylen:
|
|
wf_rolled = np.roll(wf_rolled_x, -cutoff, axis=1)
|
|
orb_rolled = []
|
|
for orbital in orb_rolled_x:
|
|
orb_rolled.append(np.roll(orbital, j - cutoff, axis=1))
|
|
ymin = j - cutoff * 2
|
|
ymax = j
|
|
else:
|
|
print(f"ERROR: for i = {i} something with rolling arrays along y goes wrong")
|
|
|
|
integral = []
|
|
for orbital in orb_rolled:
|
|
integral.append(np.linalg.norm(wf_rolled[xmin:xmax, ymin:ymax, zmin:zmax]*\
|
|
orbital[xmin:xmax, ymin:ymax, zmin:zmax]))
|
|
overlap_integrals_squared[i][j] = max(integral) ** 2
|
|
return overlap_integrals_squared
|
|
|
|
def save_as_cube(self, array, name, dir):
|
|
# TODO: rewrite
|
|
import os
|
|
if not os.path.exists(dir):
|
|
print(f"Directory {dir} does not exist. Creating directory")
|
|
os.mkdir(dir)
|
|
|
|
shape = np.shape(array)
|
|
with open(self.working_folder + '/POSCAR') as inf: #TODO get data from poscar class would be better
|
|
lines = inf.readlines()
|
|
natoms = sum(map(int, lines[6].strip().split()))
|
|
atomtypes = lines[5].strip().split()
|
|
numbers_of_atoms = list(map(int, lines[6].strip().split()))
|
|
type_of_i_atom = []
|
|
basis = []
|
|
for i in [2, 3, 4]:
|
|
vector = list(map(float, lines[i].strip().split()))
|
|
basis.append(vector)
|
|
basis = np.array(basis)
|
|
for i, number in enumerate(numbers_of_atoms):
|
|
for j in range(number):
|
|
type_of_i_atom.append(atomtypes[i])
|
|
|
|
with open(dir + '/' + name + '.cube', 'w') as ouf:
|
|
ouf.write(' This file is generated using stm.py module\n')
|
|
ouf.write(' Good luck\n')
|
|
ouf.write(' ' + str(natoms) + '\t0.000\t0.000\t0.000\n')
|
|
ouf.write(' ' + str(-shape[0]) + lines[2])
|
|
ouf.write(' ' + str(-shape[1]) + lines[3])
|
|
ouf.write(' ' + str(-shape[2]) + lines[4])
|
|
for i, line in enumerate(lines[8:natoms + 8]):
|
|
coordinate = np.array(list(map(float, line.strip().split())))
|
|
if lines[7].strip() == 'Direct':
|
|
coordinate = coordinate.dot(basis)
|
|
coordinate *= constants.BOHR_RADIUS
|
|
elif lines[7].strip() == 'Cartesian':
|
|
coordinate *= constants.BOHR_RADIUS
|
|
else:
|
|
print('WARNING!!! Cannot read POSCAR correctly')
|
|
atomtype = type_of_i_atom[i]
|
|
atomnumber = list(constants.ElemNum2Name.keys())[list(constants.ElemNum2Name.values()).index(atomtype)]
|
|
ouf.write(
|
|
' ' + str(atomnumber) + '\t0.00000\t' + str(coordinate[0]) + '\t' + str(coordinate[1]) + '\t' + str(
|
|
coordinate[2]) + '\n')
|
|
counter = 0
|
|
for i in range(shape[0]):
|
|
for j in range(shape[1]):
|
|
for k in range(shape[2]):
|
|
ouf.write(str('%.5E' % array[i][j][k]) + ' ')
|
|
counter += 1
|
|
if counter % 6 == 0:
|
|
ouf.write('\n ')
|
|
ouf.write('\n ')
|
|
print(f"File {name} saved") |