COGITO

Core class for reading VASP output files and projecting DFT wavefunctions onto atomic orbitals to build a chemically interpretable tight-binding model.

class COGITO_dft.COGITO.COGITO(wavecar_dir, readmode=False, spin=0, spin_polar=False)

Bases: object

generate_TBmodel(invariant: bool = True, irreducible_grid=True, verbose=0, tag='', include_excited=1, calc_nrms=False, save_orb_converg_info: bool = True, orbfactor=1.0, num_steps=50, num_outer=4, plot_orbs=False, min_proj=0.01, band_opt=True, orb_opt=True, orb_orth=False, start_from_orbnpy=False, plot_projBS=False, plot_projDOS=False, orbs=None, save_orb_data=False, save_orb_figs=False, minimum_orb_energy: float = -60, min_duplicate_energy: float = -60, file_type: str = 'npy')

Runs all the functions neccessary to generate the TB interpolation. REQUIRES UNIFORM KPT GRID WITH NO SYMMETRY OR FULL SYMMETRY FOR TB MODEL.

Parameters:

• invariant (bool) – Whether to construct the orbitals as eigenvectors of the local environment tensor.

• irreducible_grid (bool) – Whether or not the kpoint grid is irreducible. True for ISYM=1``|``2``|``3; False for ISYM=-1 (does not work with ISYM=0).

• verbose (str) – How much to ouput, includes 0,1,2,3. Higher numbers result in more output.

• tag (bool) – Tag to be appended to the runs output files.

• include_excited (int) – How much to include excited orbital states, includes 0, 1, 2. To do this best avoid using POTCARs with semi-core states. (0) Includes no excited orbital states (only ones which are partially occupied in isolated atom). (1) Includes some excited states (if the POTCAR has the excited states + another higher one of the same l) (generally includes p orbital for d block). (2) Includes maximally recommended states (condition of 1 + if the energy listed is the same as in isolated atom) (generally includes p orbitals for 1&2nd row).

• calc_nrms (bool) – Extra error analysis that requires more runtime to generate. Correlates well with the square root projection spilling. Both spilling and nrmse are output to “error_output”+tag+”.txt” during all runs.

• save_orb_converg_info (bool) – Creates the file ‘orb_converg_info.json’.

• orbfactor (float) – Multipled by the radial part of the pseudo radial orbital to either shrink or grow them.

• num_steps (int) – Maximum number of steps to perform the convergence. Generally convergence is achieved around 5-20 steps. Setting equal to 0 will run the standard direct algorithm where |X_a> = sum_n(c_na |Y_n>). Where |Y_n> are the set of bands (equal to number of orbitals) of the highest projection.

• num_outer (bool) – The number of outer loops in COGITO. Aka, the guassian fit from fit_to_atomic_orb is feed back into the convergence and fitting. Converges very quickly with some possibly of divergence in some systems if num_outer is too large. Values of 2-5 are recommended. Running with 0 creates model from just VASP psuedo orbitals.

• plot_orbs (bool) – Plot the radial converged orbital being fit to Gaussian functions.

• min_proj (bool) – Sets the minimum projectability of a wavefunction to be included Hamiltonian, overlap, etc matrix construction.

• band_opt (bool) – Whether or not to apply Lowdin+Gram+extras band optimization scheme for improved TB interpolation (improves band error from ~40meV to ~2meV).

• orb_opt (bool) – Whether or not enforce zero orbital mixing for improved TB interpolation (ends up having little effect).

• orb_orth (bool) – Make orthogonal tight binding (indirectly creates orthogonal basis) from Lowdin orthogonalized orbital coefficients.

• start_from_orbnpy (bool) – Start from the self.directory+”orbitals”+self.tag+”.npy” of a previous run.

• plot_projBS (bool) – Plot the orbital projected bandstructure (not BS if the kpt grid is uniform).

• plot_projDOS (bool) – Plot the orbital projected density of states (not correct DOS unless kpt grid is uniform).

• orbs (dictionary) – The orbitals to plot the projection of. Defaults to all orbitals. FORMAT: {“element1”:[orbital types]} e.g. {“Si”:[“s”,”p”],”C”:[“s”,”p”]}.

• save_orb_data (bool) – Saves the numerical orbital data for converging basis on line between central and nearest atom. Not generally useful but used to make COGITO variance figure. Best when uncommenting make_fit_wannier() below and running on dense grid.

• save_orb_figs (bool) – Saves the plot of numerical orbital data for the converging basis. Not easily interpreted.

• minimum_orb_energy (float) – The lower limit to add semi-core states in the POTCAR into the COGITO basis. Uses the atomic orbital energy listed in POTCAR.

• min_duplicate_energy (float) – The lower limit to add semi-core states when there is another valence state of the same l quantum number in the POTCAR into the COGITO basis. Uses the atomic orbital energy listed in POTCAR.

• file_type (str) – The lower limit to add semi-core states in the POTCAR into the COGITO basis. Uses the atomic orbital energy listed in POTCAR.

Usage:

new_model = COGITO(“silicon/”) new_model.generate_TBmodel(plot_orbs=True)

test_initialorbs(verbose=0, plot_orbs=False, include_excited=1, low_factor=0.8, high_factor=1.2, num_fac=9, num_outer=4, tag='', min_proj=0.01, num_steps=50, minimum_orb_energy: float = -60, min_duplicate_energy: float = -60)

Tests dependance of orbital radius and quality on the size of initial orbitals. Only runs the convergence of the orbitals without projecting them or generating the TB model.

Parameters:

• verbose (str) – How much to ouput, includes 0,1,2,3. Higher numbers result in more output.

• plot_orbs (bool) – Plot the radial converged orbital being fit to Gaussian functions.

• low_factor (float) – Minimum to multiply the orbital size by; Default is 80%.

• high_factor (float) – Maximum to multiply the orbital size by; Default is 120%.

• num_fac (int) – The number of steps between low_factor and high_factor to try.

Usage:

new_model = COGITO(“silicon/”) new_model.generate_TBmodel(plot_orbs=True)

get_WFdata_fromPOT()

Function which reads the pseudo and ae orbital information from the POTCAR. Creates the global variables self.atmProjectordata, self.atmRecipProjdata, and self.atmPsuedoAeData. FORMATs: self.atmProjectordata [=] {“element”,[float(cutoff_rad),{“orbtypes”:[proj_real_radial_data]}]} e.g. {“Si”,[1.45,{“s”:[],”s_ex”:[],”p”:[],”p_ex”:[]}} self.atmRecipProjdata [=] {“element”,[float(proj_gmax),{“orbtypes”:[proj_recip_radial_data]}]} self.atmPsuedoAeData [=] {“element”,[radial_grid, {“orbtypes”:[pseudo_radial_data]}}, {“orbtypes”:[ae_radial_data]}]} e.g. {“Si,[[0,0.01,…,1.45],{“s”:[],”s_ex”:[],”p”:[],”p_ex”:[]},{“s”:[],”s_ex”:[],”p”:[],”p_ex”:[]}]

get_local_environment()

This function gets the local environments of each atom in terms of local atom theta, phi, and radius. This is used to project onto the spherical harmonics to find the optimal rotation. Return:

get_orbs_fromPOT(orbfactor=1.0, invariant=True)

Creates the 3D initial orbitals in real space from the pseudo radial orbitals. Also sets the amount of orbitals by looping through atoms and their orbitals, thus intializes many orbital-dependent variables.

Parameters:

orbfactor (float) – Multipled by the radial part of the pseudo radial orbital to either shrink or grow them.

get_coefficients(orbitalWF, crystalWF, overlap_matrix, recip=False, band=None)

Find the coefficients for the amount of each pseudo orbital (Φ) in the pseudo wavefunction (Ψ). If the overlap is identity, the coefficients would just be the integral of the orbital-wavefunction overlap: Φa*Ψn = O_an. When the overlap is not identity S_ab, obtaining the coefficients (C_bn) requires solving the linear problem S_ab*C_bn = O_an.

Parameters:

• orbitalWF (dict) – of length M: All the orbital functions in a dictionary. The orbital may be defined in real 3D space on reciprocal G vectors.

• crystalWF (dict) – of length N: All the DFT wavefunctions in a dictionary. The wavefunction is defined on the same space as orbitalWF.

• overlap_matrix (MxM) – matrix of complex float: The overlap of the orbitals. NOTE: The orbitals and their overlaps have a k-dependence.

• recip (bool) – Whether the wavefunctions are defined in real or reciprocal space.

• band (int) – If the coefficents of only one band in the crystalWF dict is needed, pass that band as an integer here.

Returns:

The coefficients C_an of orbital a in band n.

Return type:

unknown

get_aecoefficients(orbitalWF, crystalWF, aeoverlap_matrix, kpt, recip=False, band=None, full_kpt=False, prints=False, gpnts=None, set_gpnts=False)

Find the coefficients for the amount of each ae orbital (Φ) in the ae wavefunction (Ψ). If the overlap is identity, the coefficients would just be the integral of the orbital-wavefunction overlap: Φa*Ψn = O_an. When the overlap is not identity S_ab, obtaining the coefficients (C_bn) requires solving the linear problem S_ab*C_bn = O_an. The orbital-wavefunction overlap is modified from the pseudo overlap using standard PAW methods.

Parameters:

• orbitalWF (dict) – of length M: All the pseudo orbital functions in a dictionary. The orbital may be defined in real 3D space on reciprocal G vectors.

• crystalWF (dict) – of length N: All the pseudo DFT wavefunctions in a dictionary. The wavefunction is defined on the same space as orbitalWF.

• aeoverlap_matrix (MxM) – matrix of complex float: The overlap of the ae orbitals. NOTE: The orbitals and their overlaps have a k-dependence.

• recip (bool) – Whether the wavefunctions are defined in real or reciprocal space.

• band (int) – If the coefficents of only one band in the crystalWF dict is needed, pass that band as an integer here.

Returns:

The coefficients C_na of orbital a in band n.

Return type:

unknown

get_ae_overlap_matrix(orbitalWF, secondWF=None, secondisarray=False, recip=True, kpt=0)

Finds the overlap matrix S_ab = Φa*φb. If secondWF is not defined, φ = Φ.

Parameters:

• orbitalWF (dict) – of length M: All the orbital functions in a dictionary. The orbital may be defined in real 3D space on reciprocal G vectors.

• secondWF (dict) – of length N: The second orbital functions in a dictionary. Defined on the same space as orbitalWF.

• recip (bool) – Whether the orbitals are defined in real or reciprocal space.

Returns:

The overlap matrix S_ab for orbitals a and b.

Return type:

unknown

get_overlap_matrix(orbitalWF, secondWF=None, secondisarray=False, recip=False)

Finds the overlap matrix S_ab = Φa*φb. If secondWF is not defined, φ = Φ.

Parameters:

• orbitalWF (dict) – of length M: All the orbital functions in a dictionary. The orbital may be defined in real 3D space on reciprocal G vectors.

• secondWF (dict) – of length N: The second orbital functions in a dictionary. Defined on the same space as orbitalWF.

• recip (bool) – Whether the orbitals are defined in real or reciprocal space.

Returns:

The overlap matrix S_ab for orbitals a and b.

Return type:

unknown

lowdin_orth(low_orbitals, set_overlap=False, overlap=None, recip=False)

Orthogonalized the orbital based on the Lowdin scheme. The new orbitals Ψ are defined by the original orbitals Φ as Ψ_b = conj(S_ab)^(-1/2)*Φ_a.

Parameters:

• low_orbitals (dict) – of length M: All the orbital functions in a dictionary. The orbital may be defined in real 3D space on reciprocal G vectors.

• set_overlap (bool) – Whether the orbital overlap is being passed to the function (True) or should be calculated (False).

• overlap (MxM) – matrix of complex float: The overlap of the orbitals. NOTE: The orbitals and their overlaps have a k-dependence.

• recip (bool) – Whether the orbitals are defined in real or reciprocal space.

Returns:

Lowdin orthogonalized orbitals.

Return type:

unknown

converge_orbs_recip(num_steps=50)

Converge to the atomic-like Bloch orbitals which best fit the DFT wavefunction. Please reference COGITO paper for theoretical description.

Creates global variable one_orbitalWF which is referenced in the Bloch to atomic orbital fitting.

Parameters:

num_steps (int) – Maximum number of steps to perform the convergence. Setting equal to 0 with run the standard direct algorithm where |X_a> = sum_n(c_na |Y_n>). Where |Y_n> are the set of band (equal to number of orbitals) of the highest projection.

fit_to_atomic_orb(plot_orbs=False, orbfactor=1.0, outer=0)

This function is necessary when using reciprocal integrals because they need an individual orbital rather than the unk = sum_R(y(r-R)). The algorithm used here determines the percent of the center orbital which is exact if the orbital has a precise exponential decay with radius. See COGITO paper for more details.

spher_bessel_trans(orbital_coeffs, rad_grid_size=500)

center_real_orbs(real_orbs, kpt=0, make_real=False)

recip_to_real(recip_orbs, kpt=0, make_real=False)

real_to_recip(real_orbs, kpt=0)

get_kdep_recipprojs(kpt, full_kpt=False, gpnts=None, set_gpnts=False)

get_kdep_reciporbs(kpt, full_kpt=False, gpnts=None, set_gpnts=False, orb_params=[], orb_vecs=[])

proj_all_kpoints(max_bandavg=None, calc_nrms=False)

optimize_band_set(band_opt=True, orb_opt=True, orb_orth=False)

This funciton refines the orbital coefficients for perfect ‘spilling’ in the valence bands and improving spilling and correct orthogonalize to VB in the conduction bands. This process significantly decreases the band error in the TB interpolation while making little to no change to the chemically interpretible metrics (i.e. COHP, COOP, etc).

Returns:

Set of band which minimizes the difference between the orbital states left after and the band set orbital states; maybe later: also minimize overlap between Lowdin orthogonalized band sets.

Return type:

unknown

expand_irred_kgrid()

This function does a couple things: 1. Finds the kpoints of the reducible grid and the coorespond symmetry operations to get them from the irreducible points. 2. Creates the eigenvalues, eigenvectors, and overlaps matrices for the new reducible kpoint grid.

symmetrize_orbs(recip_orbs)

This function is to symmetrize the iterated orbitals to decrease orbital mixing by ensuring orbital has s orbital symmetry.

get_Qab()

get_ae_overlap_info(orbs=None, kpt=0, just_ae_overlap=False, recip=False, test=False, full_kpt=False, gpnts=None, set_gpnts=False)

get_orth_coeffs(coeff, kpoint=None)

get_hamiltonian(kpoints=None)

make_fit_wannier()

get_TBparameter(file_type: str = 'npy')

get_interp_ham(kind, return_truevec=True, return_params=False)

get_neighbors()

orthogonalize_basis()

get_offset(a=2)

write_TBfiles(file_type: str = 'npy')

write_recip_rad_orbs()

read_recip_rad_orbs(tag='')

get_proj_on_aeorb()

write_input_file()

plot_BS()

plot_projectedBS(orbs)

plot_projectedDOS(orbs, xlim=None)

generate_gpnts(kpt)

COGITO_dft.COGITO.func_for_rad(x, a, b, c, d, e, f, g, h, l)

COGITO_dft.COGITO.func_for_rad_fit(x, a, b, con1, con2, con3, c, g, con4, l)

COGITO_dft.COGITO.func_for_rad_exp(x, a, b, c, d, e, f, l)

COGITO_dft.COGITO.lowdin_orth_vectors(vectors)

COGITO_dft.COGITO.lowdin_orth_vectors_orblap(vectors, orblap, energies=None, return_energy=False)

COGITO_dft.COGITO.GS_orth_twoLoworthSets(vectors1, vectors2)

COGITO_dft.COGITO.GS_orth_twoLoworthSets_orblap(vectors1, vectors2, orblap, energy1=None, energy2=None, return_energy=False, energycut=0, ratio=[])

COGITO_dft.COGITO.GS_combine_states_orblap(vectors1, vectors2, orblap, energy1=None, energy2=None, return_energy=False)

COGITO_dft.COGITO.complex128funs(phi, theta, sphharm_key)

COGITO_dft.COGITO.normalize_wf(wavefunc, prim_vec, gridnum, return_integral=False, recip=False)

COGITO_dft.COGITO.periodic_integral_3d(f, prim_vec, n, multiple_wfs=False)

COGITO_dft.COGITO.reciprocal_integral(f)

COGITO_dft.COGITO.smooth(y, box_pts)

COGITO_dft.COGITO.extract_line_data(grid, gridXYZ, p1, p2, tolerance=0.5)

Extracts orbital magnitude along a line passing through p1 and p2 in a non-Cartesian 3D grid.

Parameters:

• grid (numpy.ndarray) – 3D array of orbital magnitudes.

• gridXYZ (numpy.ndarray) – 2D array (3, num_points) of real-space coordinates.

• p1 (tuple) – First point (x1, y1, z1) in real space.

• p2 (tuple) – Second point (x2, y2, z2) in real space.

• tolerance (float) – Distance threshold to include points near the line.

Returns:

Distances along the line.

magnitudes (numpy.ndarray): Orbital magnitudes at corresponding distances.

Return type:

distances (numpy.ndarray)

COGITO_dft.COGITO.combine_and_save_plots(plots, filename='combined_plot.png', layout=None)

Combines multiple Matplotlib plots into a single figure and saves to a file.

Parameters:

• plots (list of matplotlib.figure.Figure) – List of Matplotlib figures.

• filename (str) – Output filename (supports .png, .pdf, .svg, etc.).

• layout (tuple or str) – (rows, cols) for custom layout or “auto” for automatic grid.

Returns:

None

COGITO_dft.COGITO.plot_matrix(matrix, low_center=0.2, high_center=0.85, filename='matrix.png')

COGITO_dft.COGITO.D_real_l3_from_R(R)

Compute the Wigner D matrix for f orbitals (l=3) from rotation matrix R using properly symmetrized and trace-free rank-3 tensors.

Parameters: R : 3x3 rotation matrix

Returns: D_matrix : 7x7 matrix mapping old f orbitals to new f orbitals Order: f_z3, f_xz2, f_yz2, f_xyz, f_z(x2-y2), f_x(x2-3y2), f_y(3x2-y2)

class COGITO_dft.COGITO.Tee(*streams)

Bases: object

write(data)

flush()

COGITO_dft.COGITO.run_cogito(directory, save_metadata: bool = True, invariant=True, irreducible_grid=True, verbose=0, tag='', include_excited=1, calc_nrms=False, save_orb_converg_info: bool = True, orbfactor=1.0, num_steps=50, num_outer=4, plot_orbs=False, min_proj=0.01, band_opt=True, orb_opt=True, orb_orth=False, start_from_orbnpy=False, plot_projBS=False, plot_projDOS=False, orbs=None, save_orb_data=False, save_orb_figs=False, minimum_orb_energy: float = -60, min_duplicate_energy: float = -60)

Runs COGITO. See generate_TBmodel() for a description of all the other arguments.

Parameters:

directory (str) – The directory that has all the input files.

Usage:

from COGITO import run_cogito run_cogito(“silicon”)

COGITO_dft.COGITO.main(argv=None)