mutable struct tb_k{T}

Tight binding object in k-space. Can be from direct import of DFT band
structure using atomic proj, or from fft'ed tb object. Similar to real-space version.

# Holds
- `Hk::Array{Complex{T},4}` Hamiltonian in k-space
- `K::Array{Float64,2}` K point array. In fractional coordinates of BZ.
- `kweights::Array{Float64,1}` K point weights.
- `k_dict::Dict` Dictionary from k-point like [0,0,0] to index.
- `nwan::Int64` Number of orbitals / generalized wannier functions.
- `nk::Int64` Number of k-points.
- `nspin::Int64` Number of spins (2=magnetic)
- `nonorth::Bool` 
- `Sk::Array{Complex{T},3}` Overlap matrix.
- `scf::Bool` needs self-consistency?
- `h1::Array{T,3} #scf term` holds scf term if present
- `grid::Array{Int64,1}` dimensions of k-point grid, from regular grid like `[8,8,8]`

   mutable struct tb_crys_kspace{T}

Hold k-point tight binding and crystal structure. Similar to `tb_crys`

# Holds
- `tb::tb_k`
- `crys::crystal`
- `nelec::Float64`
- `nspin::Int64`
- `dftenergy::Float64`
- `scf::Bool`
- `gamma::Array{T, 2}`
- `eden::Array{Float64,2}`

function read_tb_crys_kspace(filename; directory=missing)

Reads and returns from filename a tb_crys_kspace object. See write_tb_crys_kspace

If cannot find "filename", will look for "filename.xml", "filename.gz", "filename.xml.gz"

Can read gzipped files directly.

function write_tb_crys_kspace(filename, tbc::tb_crys_kspace)

Save a tb_crys_kspace object to xml format. See read_tb_crys_kspace

    function make_tb_crys(ham::tb,crys::crystal, nelec::Float64, dftenergy::Float64; scf=false, eden = missing, gamma=missing, within_fit=true, screening=1.0, tb_energy=-999, fermi_energy=0.0 )

Constructor function for `tb_crys` object

    function make_tb_crys_kspace(hamk::tb_k,crys::crystal, nelec::Float64, dftenergy::Float64; scf=false, eden = missing, gamma=missing, screening=1.0)

Constructor function for `tb_crys_kspace` object

function make_tb(H, ind_arr, r_dict::Dict; h1=missing)

Constructor function for tb

function make_tb(H, ind_arr, r_dict::Dict, S; h1=missing)

Constructor function for tb with overlaps

function make_tb(H, ind_arr, S; h1=missing)

Constructor function for tb, better programming.

function make_tb_k(Hk, K, kweights, Sk; h1=missing, grid=[0,0,0], nonorth=true)

Constructor for tb_kspace

function load_hr_dat(filename, directory="")

Load a wannier90 hr.dat file Not currently a major part of program, but you can use if you want.

function calc_energy_fft(tbc::tb_crys; grid=missing, smearing=0.01, return_more_info=false)

Get energy using fft.

returns energy

if return_more_info==true then returns etot, efermi, vals, vects

function calc_energy_charge_fft(tbc::tb_crys; grid=missing, smearing=0.01)

Do fft, then calculate energy and charge.

     function calc_energy_charge_fft(tbc::tb_crys_sparse; grid=missing, smearing=0.01)

Sparse matrix non-scf energy/eigenvectors.

function calc_energy(tbc::tb_crys; smearing=0.01, returnk=false)

Calculate energy without fft.

function calc_energy(h::tb_crys, kgrid; smearing=0.01, returnk=false)

Calculate energy no fft

function types_energy(tbc::tb_crys)

Calculate the reference energy of each atom in crystal. This results in the total energy being indexed to seperated non-spin-polarized atoms. This is arbirary.

function types_energy(tbc::tb_crys)

function types_energy(c::crystal)

function types_energy(types)

function make_kgrid(kgrid)

-kgrid is an array of 3 integers like [8,8,8]

returns regular MP Gamma-centered k-point grid and (equal) k-weights.

function trim(h::tb, tol=0.0002)

Remove terms in tb with abs value smaller than tol. (Ryd) Will speed calculations at cost of accuracy. USE WITH CARE. Usually not necessary except for very large or very detailed calculations.

function organizedata(tbc::tb_crys)

Rearrange the data in the tbc as a function of distances for plotting purposes.

return data_onsite, data_arr

Returns two arrays. The first has data on the onsite elements.

• column 1 and 3 have atom indexes
• column 2 and 4 have orbital numbers
• columns 5 and 6 have real and imaginary parts of H
• columns 11 and 12 have real and imaginary parts of S
• column 7 has the closest inter-atomic distance
• column 8 has the index of the closest atom.

The second has intersite data

• column 1 and 3 have atom indexes for the atom pairs
• column 2 and 4 have orbital numbers
• columns 5 and 6 have real and imaginary parts of H
• columns 11 and 12 have real and imaginary parts of S
• column 7 has the inter-atomic distance
• column 8,9,10 have the the direction cosines lmn for the atom pair.

function organizedata(crys::crystal, h::tb)

function myfft_R_to_K(tbc, grid=missing)

Does Fourier Transform R->k (fft) using FFTW for tb_crys

function myfft(crys, nonorth, grid, kpts,ham_kS, Sk=missing)

Does Fourier Transform K->R (ifft) using FFTW.

Arguments

• crys crystal
• nonorth nonorogonal bool
• grid k-point grid size
• kpts the k-points nkpts×3 in the original order, to be rearranged into grid
• ham_kS hamiltonian in k space (nw×nw×nkpts)
• Sk overlaps in k space

function get_sym_R(crys, grid, sss = 1.0)

Figures out the r-space grid using Wigner-Seitz like construction to figure out the best arrangement of r-grid points to keep periodic copies closest to the original atom and take into account symmetry.

return R_grid, R_int_grid, sym_R

returns the R_grid, the integer version, and the symmetry factor of each point.

function tb_indexes(d::dftout)

Figures out mapping between DFT projected hamiltonian orbitals and crystal and the wannier orbitals we want.

return wan, semicore, nwan, nsemi, wan_atom, atom_wan

• wan has the indexes of the wannier orbitals
• semicore has the indexes of semicore states.
• nwan number of wannier orbs
• nsemi number of semicore states
• wan_atom dictionary wannier to atom numbers
• atom_wan dictionary atom numbers to wannier orbitals

function find_vbm_cbm(eigs, fermi)

Find the valence band max and conduction band minimum from eigs, relative to Fermi level.

function ewald_energy(tbc::tb_crys, delta_q=missing)

Return ewald energy term from tbc. If delta_q, the atomic charge density, is missing, loads from tbc.

function ewald_energy(tbc::tb_crys_kspace, delta_q=missing)

function ewald_energy(crys::crystal, gamma, delta_q::Array{Float64,1})

Does the actual calculation.

function get_dq(tbc::tb_crys_kspace)

Get atomic charge density from tb_crys or tb_crys_kspace or crys + eden

function get_dq(tbc::tb_crys)

function get_dq(crys::crystal, chargeden::Array{Float64,1})

function get_energy_electron_density_kspace(tbcK::tb_crys_kspace; smearing = 0.01)

Get energy / charge density from k-space tight binding object.

return bandenergy + etypes + echarge + energy_smear, eden, VECTS, VALS, error_flag

function get_energy_electron_density_kspace(tb_k::tb_k, nelec; smearing = 0.01)

K-space get energy and electron density from tb_k

module AtomicProj

Analyze projections from QE projwfc.x

mutable struct proj_dat

Holds data from projwfc.x

• bs::bandstructure
• natwfc::Int64 number of atomic wavefunctions, from dft projection
• proj::Array{Complex{Float64}, 3} atomc projections nk × natwfc × nbnd , where nbnd is the number of bands in DFT
• overlaps::Array{Complex{Float64}, 3} overlap matrix from dftnk × natwfc × natwfc`
• nspin::Int64 number of psins

This is created from loadXML_proj, which calls make_proj

function create_tb(p::proj_dat, d::dftout; energy_froz=missing, nfroz=0, shift_energy=true)

Does the main creation of TB hamiltonian from DFT projection data in k-space.

Arguments

• p::proj_dat Projection data
• d::dftout DFT scf data
• energy_froz=missing Energy to start freezing eigenvalues
• nfroz=0 number of frozen bands.
• shift_energy=true if true shift eigenvalues so band energy == total energy

function get_ham_r(tbck::tb_crys_kspace)

Do fft to get the real-space ham. Requires tbck to be on a regular grid centered at Gamma with no symmetry.

function get_ham_r_slow(p::proj_dat, d::dftout, grid, ham_k::Array{Complex{Float64}, 3}; nonorth=true, K=missing)

Slow method for doing fourier transform. Uses standard ft formula, not fft.

Use fast version instead.

function loadXML_proj(savedir, B=missing)

Load proj from QE output file. Need the QE save dir like "qe.save". B are reciprocal lattice vectors.

function makeOG(prefix, tmpdir )

opengrid.x inputfile. My workflow doesn't currently use this, as opengrid.x was unreliable, and I've rewritten the code to avoid it and work directly with independent k-points only.

function make_commands(nprocs=1)

Returns a dictionary with command lines to call external programs on the command line nprocs is the number of processors for parallel execution Needs to be changed for you specfic program locations and mpi commands (if any)

This needs to be edited to actually run QE yourself. Running wannier90 is optional, not part of current code.

function make_proj(bs, proj, overlaps)

Constructor for proj_dat.

function make_projwfcx(prefix, tmpdir)

Create the inputfile for running projwfc.x from QE

function makedict_proj(savedir)

Return xml proj data from QE savedir

function prepare_ham_k(p::proj_dat, d::dftout, grid, ham_k::Array{Complex{Float64}, 4}; nonorth=true, K=missing, localized_factor = 0.05, screening=1.0)

Constructs the actual k-space hamiltonain, which involves lots of putting matricies in the correct form and de-orthogonalizing if desired.

#arguments

• p::proj_dat projection data
• d::dftout dft data
• grid kpoint grid spacing
• ham_k::Array{Complex{Float64}, 3} Hamiltonian
• nonorth=true make a non-orthogonal TB
• K=missing k-point array, usually get from p
• localized_factor = 0.15 increase localization of overlaps.
• screening=1.0 mulitply U by this factor. usually not used.

function projwfc_workf(dft::dftout)

This is the main workflow for the creation of TB matrix elements from QE DFT calculations.

Starting from a converged QE scf calculation...

1. Run NSCF calculation with extra empty bands. If you want a real space TB object, these need to be a full k-point grid not using symmtery. k-space only can use irreducible k-points

2. Run projwfc.x

3. Construct the TB hamiltonian in k-space

4. (Optional) FT to real space

5. (Optional) cleanup wavefunctions.

Arguments

• dft::dftout The starting scf calculation
• directory="./"
• nprocs=1 number of processors
• freeze=true Keep occupied eigenvalues fixed to exact DFT values
• writefile="projham.xml" output file for real-space TB
• writefilek="projham_K.xml" output file for k-space TB
• skip_og=true Not used anymore
• skip_proj=true If projections are already run, don't run them again, load from file.
• shift_energy=true Shift energy of eigenvalues s.t. total energy equal band energy.
• cleanup=true Remove the large wavefunction files from disk, keeping nscf/projection files.
• skip_nscf=true If nscf calculation is already done, load results from file.
• localized_factor = 0.15 Adjust extent of overlap matrix. 0.0 uses full atomic wavefunctions, which can be overly delocalized. 1.0 is fully localized.
• only_kspace=false Do not create real-space tb. Usually true in current code, as I can fit directly from k-space tb only.
• screening = 1.0 If use a screening factor to reduce value of U in Ewald calculation. Usually leave at 1.0.

function run_og(filename="og.in";  directory="./", nprocs=1)

runs open_grid.x

function run_projwfcx(projfile="proj.in"; directory="./", nprocs=1)

Run projwfc.x from QE code

function band_summary(tbc, kgrid, fermi=missing)

Produces summary of band structure. See below functions for more specific versions of function that automatically generate the k-points.

Note: gaps are not well-defined for non-magnetic systems with odd numbers of electrons, as they are required to be metals.

Returns direct_gap, indirect_gap, gaptype, bandwidth

-direct_gap: minimum gap at one k-point between nominally filled and empty bands. Can be non-zero in metals. -indirect_gap: LUMO - HOMO. Can be negative if material has a direct gap everywhere, but the conduction band at some k-point is below the valence band at a different k-point. Physically these are indirect gap semimetals. -gaptype : is :metal for all metals, :direct or :indirect for insulators. -bandwidth : HOMO - minimumbandenergy. Included semicore states if they are in the TB calculation.

function band_summary(tbc::tb_crys_kspace)

Will use internal k-points by default.

function band_summary(tbc::tb_crys; kgrid=missing, kpts=missing)

Will automatically generate standard k-grid by default.

function get_kpath(kpath=[0.5 0 0 ; 0 0 0; 0.5 0.5 0.5], names = missing, npts=30)

Construct a k_path for a band structure calculations. Very simple.

• kpath high symmetry k-points in fractional BZ coordinates.
• names names of kpoints like ["Γ", "X"]
• npts number of points between high-symmetry kpoints

function plot_bandstr(h::tb)

Plots using tb

function plot_bandstr(h::tb_crys; kpath, names = missing, proj_types=missing, proj_orbs = missing, proj_nums=missing)

Plot the band structure of a tb_crys object. Can also perform a projected band structure if you specify at least one of proj_types, proj_orbs, proj_nums.

k-path specified by a kpath array and names.

Must do scf calculation before plotting.

Arguments

• h::tb_crys - The tight-biding object we want to plot bands from. Only required argument.
• kpath=[0.5 0 0 ; 0 0 0; 0.5 0.5 0.5; 0 0.5 0.5; 0 0 0 ;0 0 0.5] - nk × 3 array k-point path (high symmetry points).
• npts=30, - number of points between high-symmetry k-points.
• names=missing - nk string array. Names of the high-symmetry k-points
• proj_types=missing - types to project onto. Either proj_types="H" or proj_types=["H", "O"] are valid.
• proj_orbs=missing - orbitals to project onto. either proj_orbs=:s or proj_orbs=[:s, :p].
• proj_nums=missing - atom numbers to project onto. Either proj_nums=1 or proj_nums=[1, 2]
• efermi=missing - allows you to specify fermi energy. Default is to take from h
• color="blue" - specify line color
• MarkerSize=missing" - specify markersize
• yrange=missing" - specify y-range. e.g. yrange=[-0.7, 0.3]
• plot_hk=false - plot things besides the normal band structure. Can be one of :Seig, :Heig, :Hreal, :Himag, :Sreal, :Simag to plot H or S eigvals or components. Primarily for debugging.
• align="vbm" - default or "valence" is to align valence band max to zero energy. Can also be "min", which aligns on the minimum eigenvalue, or "fermi" or "ef", which align on the Fermi level,
• clear_pervious=true - clears the plot before adding new stuff.
• do_display=true - display the plot. If false, can be used with display-less nodes. You can still use savefig from Plots to produce saved images.

function plot_bandstr_dos(c::crystal; sym_prec = 5e-4, npts=30, efermi = missing,color="blue", MarkerSize=missing,yrange=missing,plot_hk=false, align = "fermi",  proj_types = missing, proj_orbs = missing, proj_nums=missing, clear_previous=true, do_display=true, color_spin = ["green", "orange"],  spin = :both, nspin = 1, database=missing, smearing = 0.025)

Run SCF in standardized crystal structure and calculate band structure along symmetry-derived k-lines, as well as the DOS. Then plot them layed out side-by-side.

See plot_bandstr_sym and dos

function plot_bandstr_dos(tbc::tb_crys)

Plot symmetry-derived k-lines and DOS layout out together. This version takes in tb_crys object from a previous SCF calculation. Will need re-run SCF if tbc.crys is not standardized.

See plot_bandstr_sym and dos

function plot_bandstr_sym(c::crystal; sym_prec = 5e-4, npts=30, efermi = missing, color="blue", MarkerSize=missing, yrange=missing, plot_hk=false, align = "vbm", proj_types = missing, proj_orbs = missing, proj_nums=missing, clear_previous=true, do_display=true, color_spin = ["green", "orange"], spin = :both, nspin = 1, database=missing)

Plots the band structure using a set of K-points determined by the space group. Will standardize the crystal structure and run the SCF calculation automatically. Conventions similar to Setyawan and Curtarolo CMS 2010 as well as the jarvis-tools package. Other options similar to plot_bandstr

function plot_bandstr_sym(tbc::tb_crys;  sym_prec = 5e-4, npts=30, efermi = missing, color="blue", MarkerSize=missing, yrange=missing, plot_hk=false, align = "vbm", proj_types = missing, proj_orbs = missing, proj_nums=missing, clear_previous=true, do_display=true, color_spin = ["green", "orange"], spin = :both, nspin = 1, database=missing)

Plots the band structure using a set of K-points determined by the space group. This version takes in a tight-binding crystal object from a pervious SCF calculation. If the crystal is in the standardized unit cell, then it will not repeat the SCF calculation, but otherwise it will.

function plot_compare_dft(tbc, bs; tbc2=missing)

Plots a band structure comparison between a tight-binding crystal object (tb_crys) and a band structure directly from dft (either a dftout or bs object).

The k-points are fixed by the bs object.

tbc2 is an optional second tbc_crys.

function plot_compare_tb(h1::tb, h2::tb; h3=missing)

function plot_compare_tb(h1::tb_crys, h2::tb_crys; h3=missing)

Plot a comparison between different tight binding objects h1, h2, and optionally h3. Options similar to plot_bandstr but more limited.

function run_dft_compare(tbc; nprocs=1, prefix="qe",   outdir="./", kpath=[0.5 0 0 ; 0 0 0; 0.5 0.0 0.5], names = missing, npts=30, efermi = missing, yrange=missing,   align="vbm")

This function will run a new QE dft calculation and band structure calculation and compare the bands with the tight binding model tbc.

outdir is the location the QE files will be stored.

Other options like the other plotting commands

function setup_proj(crys, nwan, proj_types, proj_orbs, proj_nums)

Figure out projection indexes

function dos(tbc::tb_crys; grid=missing, smearing=0.04, npts=missing, proj_type=missing, do_display=true)

DOS, mostly for testing.

• npts is number of energies
• proj_type can be "none", "atomic", or "orbital"
• do_display=false will suppress the actual plot

The combination of smearing and grid are important to get converged results.

See also dos

return energies, dos, projected_dos, pdos_names

function dos_realspace(tbc; direction=3, grid=missing, smearing=0.005, npts=missing, do_display=true, use_sym=false, energy_grid=100, width=missing, energy_lims=missing, scissors_shift=0.0, scissors_shift_atoms=[])

Makes an atom and spatially resoloved DOS-style plot. Along direction=1,2,3 (lattice vector 1,2,3) the atoms are plotted, using transparancy to denote DOS amplitude. Try it out to see. Potentially useful for interfaces. Example usage: c = makecrys([8.0 0 0; 0 8.0 0; 0 0 10.0], [0 0 0; 0 0 0.5], [:Li, :Cl]) #quasi 1D LiCl chain en, tbc, flag = scf_energy(c*[1,1,10]) #ten unit cells dos_realspace(tbc, direction=3)

Adjust the energy_grid and width to adjust plotting parameters. Scissors shift in energy unit to certain atoms will add scissors shift to open band gaps is desired.

function dos_tetra(tbc::tb_crys; grid=missing, npts=missing, proj_type=missing, do_display=true)

DOS, using tetrahedral integration

• grid is the k-point grid. Defaults to 1.6 times the default grid for energy integration.
• npts is number of energies
• proj_type can be "none", "atomic", or "orbital". Defaults to atomic if more than one atom type.
• do_display=false will suppress plotting

return energies, dos, projected_dos, pdos_names

function gaussian_dos(tbc::tb_crys)

Gaussian smearing DOS is now the main DOS. Tetraheral method is still coded, but it is a bit wonky.

function plot_dos(energies, dos, pdos, names; filename=missing, do_display=true)

Does the actual DOS plotting, called by dos or gaussian_dos

function plot_dos_flip(energies, dos, pdos, names; filename=missing, do_display=true, yrange=[-6, 4.0])

For plotting with x and y axes reversed, useful for bandstucture plots

function projection(tbc::tb_crys, vects, sk3, grid; ptype=missing)

Figures out the projections. ptype can be :atomic or :orbs for atom projection or orbital projection (:s,:p,:d) Default is to choose :atomic except for elemental systems.

function setup_tetra(grid)

Setup simple tetrahedron method based on tetra.f90 in QE. Simple tetra method.

module Atomdata

Holds the data for isolated atoms.

atom_prefered_oxidation::Dict()

Prefered oxidation states of atoms, in descending order of preference (approximate).

atom_radius::Dict()

Atomic radius (metallic), in pm

atoms::Dict()

Periodic table information.

cutoff_dist::Dict()

Cutoff distance for Hamiltonian calculation

electronegativity::Dict()

Pauling scale.

min_dimer_dist_dict::Dict()

Approximate minimum distance in current fitting data.

sub_list::Dict()

Most likely substitution list, approximate.

module AtomicMod

Simple module for atom type.

struct atom

Hold basic atomic information

• name::String
• Z::Int64 Atomic number
• row::Float64 Periodic table row
• col::Float64 Periodic table col
• mass::Float64 Mass in amu
• nval::Float64 Number of valence electrons in TB calculation.
• nsemicore::Int64 Number of semicore electrons. Depends on pseudopotential choice
• nwan::Int64 Number of TB orbitals
• orbitals::Array{Symbol,1} Names of the orbitals, like [:s, :p]
• total_energy::Float64 DFT total energy, depends on pseudopotentials.
• eigs::Dict orbital eigenvalues of isolated non-spin-polarized nuetral atom.
• energy_offset::Float64 energy to add to TB calculation to make isolated atoms have zero energy.
• U::Float64 U value for Ewald correction

function makeatom(name, Z, row, col, mass, nval, nsemicore, orbitals, etot, eigs, vac_potential=0.0, U=0.0)

Constructor for atom.

module BandTools

Utility functions for manipulating band structures. These need to be defined early in the code so other Modules have acccess to them.

function band_energy(eigs, weights, nelec, smearing = 0.01; returnk=false, returnocc=false, returnef=false, returnboth=false)

Calculate band energy. Has options for additional return variables. Calculates fermi energy internally.

function calc_fermi(eigs, weights, nelec, smearing = 0.01)

calculate fermi energy using bisection

function smearing_energy(eigs, weights, efermi, smearing = 0.01)

Smearing contribution to total energy. If you don't take this int account in metals or small gap semiconductors, your energy is not variational in the correct way, and things like forces and stresses become wrong.

module CalcTB

Create TB matrix sets from coefficients, or prepare to fit the coefficients.

struct coefs

Hold the TB coefficients for 2body or 3body interactions

• dim::Int64 - dimension (2 or 3)
• datH::Array{Float64,1} Hamiltonian parameters
• datS::Array{Float64,1} Overlap parameters, if any
• sizeH::Int64
• sizeS::Int64
• inds::Dict{Array{Symbol}, Array{Int64,1}} Holds inds that tell which coeffiecients correspond to which atoms/orbitals
• inds_int::Dict{Array{Symbol}, Any} Holds inds that tell which coeffiecients correspond to which atoms/orbitals, but integer array
• ninds_int::Dict{Array{Symbol}, Array{UInt16,2}} Holds number of coefs
• names::Set Names of atoms
• orbs::Array{Any,1} Orbitals
• cutoff::Float64 cutoff distance
• min_dist::Float64 minimum atom-atom distance in fitting data
• maxmin_val_train::Dict max and min value of matrix elements within fitting data
• dist_frontier::Dict dictionary of pareto frontier of shortest fitting distances
• version::Int64 version number

function calc_frontier(crys::crystal, frontier; var_type=Float64, test_frontier=false, diststuff=missing, verbose=true)

Calculate a pareto frontier of short distances. For 2body interactions, this is just a single distance, which is easy. For 3body interactions, there are 3 distances, so multiple points are on the frontier. Useful for deciding if old fitting data applies to a new structure with atoms that are close together.

test_frontier=true is used to check if new structure is in the frontier. othersise, see if new structure changes frontier.

function calc_onsite(t1,s1,s2, database=missing)

Handles atomic matrix els. We do not currently fit onsite matrix els

function calc_tb_LV_sparse(crys::crystal)

Main sparse matrix tight binding Hamiltonian calculator.

calc_tb_fast(crys::crystal; database=missing, use_threebody=true, use_threebody_onsite=true)

Construct tb_crys from crystal stucture, but does not solve. This is usually called internally by functions like scf_energy, but you can use it directly if you want. Until you do a SCF energy calculation, the electron density and Fermi level will be wrong.

Arguments

• crys::crystal - Required crystal structure
• database=missing - Source of coefficients. Will load from default source if not specified.
• use_threebody=true - Use three-body off-site interactions. Only turn off for testing purposes.
• use_threebody_onsite=true - Use three-body on-site interactions. Only turn off for testing purposes.
• verbose=true - set to false for less output.
• var_type=missing - variable type of tb_crys. Default is Float64.

function calc_tb_prepare_fast(reference_tbc::tb_crys; use_threebody=false, use_threebody_onsite=false)

Take a tbc from DFT and rearrange it for use in fitting code. Basically set it up so that it is ready for a least squares linear fit of coefficients.

return twobody_arrays, threebody_arrays, hvec, svec, Rvec, INDvec, h_onsite, ind_conversion, dmin_types, dmin_types3

Where

• twobody_arrays - info for fitting twobody coefficients
• threebody_array - info for fitting twobody coefficients
• hvec - vector of reference TB matrix els, arranged for fitting
• svec - vector of reference overlap matrix els, arranged for fitting
• Rvec - displaments
• INDvec info on which matrix el goes with which row
• h_onsite info for subtracting atomic terms, I think
• ind_conversion - dict to convert between place in hamiltonian and overall counter, which removes duplicates.
• dmin_types - shortest 2body distances
• dmin_types - shortest 3body distances

function calc_threebody(c,ind, t1,t2,t3,orb1,orb2,dist,dist31,dist32,lmn12, lmn31,lmn32; database=missing, memory0=missing, memory1=missing, memory2=missing, memoryV=missing, precalc=false, set_maxmin=false)

Calculate 3body intersite hamiltonian interactions.

memory1, etc have preallocated memory. This function is important for performance.

function calc_threebody_onsite(t1,t2,t3,orb1,dist12,dist13,dist23, database; set_maxmin=false, memory=missing)

calculate threebody onsite interactions

function calc_twobody(t1,t2,orb1,orb2,dist,lmn, database)

Two body intersite Hamiltonian and overlap matrix els.

function calc_twobody(t1,t2,orb1,orb2,dist,lmn, database)

Two body intersite Hamiltonian and overlap matrix els.

function calc_twobody_onsite(t1,t2,orb1,orb2, dist,lmn, database)

Calculate 2body onsite interactions.

function fit_threebody(t1,t2,t3,orb1,orb2,dist,dist31,dist32,lmn12, lmn31,lmn32)

Fit threebody intersite interactions

function fit_threebody_onsite(t1,t2,t3,orb1,dist12,dist13,dist23)

Fit three body onsite interactions.

function fit_twobody(orb1,orb2,dist,lmn)

Fit Two body intersite Hamiltonian and overlap matrix els.

function fit_twobody_onsite(t1,t2,orb1,orb2, dist,lmn)

Fit 2body onsite interactions.

function get_data_info(at_set, dim)

Figure out the arrangement of data in a coefs file.

Loops over various combinations of orbitals and atoms and assigns them places in datH and datS, depending on the terms included in the model and the dimensionaly.

function get_data_info(at_set, dim)

Figure out the arrangement of data in a coefs file.

Loops over various combinations of orbitals and atoms and assigns them places in datH and datS, depending on the terms included in the model and the dimensionaly.

function laguerre(dist, ind=missing; nmax=6, memory=missing)

Calculate laguerre polynomials up to order nmax

function make_coefs(at_list, dim; datH=missing, datS=missing, cutoff=18.01, min_dist = 3.0, fillzeros=false, maxmin_val_train=missing, dist_frontier=missing)

Constructor for coefs. Can create coefs filled with ones for testing purposes.

See coefs to understand arguments.

function plot_database(database, entry, t=missing)

Plot some data from coefs. Needs to be updated to work with Plots, probably doesn't work now.

function read_coefs(filename, directory = missing)

Read coefs from filename. Can read gzip directly.

function renormalize_S(tbc, database, cutoff=17.99)

Function for changing S but keeping the band structure fixed. Not currently used.

function symmetry_factor_fit(s1,s2,lmn)

All the Slater-Koster factors, for fitting.

function symmetry_factor(s1,s2,lmn, dat)

All of the spd Slater-Koster matrix elements. dat is preallocated memory.

function two_body_H(dist, ind=missing)

get 3body hamiltonian terms together.

function two_body_O(dist, ind=missing)

three body onsite.

function trim_dist(tbc, cutoff=18.0001)

Reduce the atom-atom hamiltonian terms longer than cutoff.

Not used in typical code, but can make tbc run faster / reduce memory.

function two_body_H(dist, ind=missing)

Two body intersite Hamiltonian.

function two_body_O(dist, ind=missing)

Two body onsite

function two_body_H(dist, ind=missing)

Two body intersite overlap.

function write_coefs(filename, co::coefs; compress=true)

Write coefs to a file. Compress uses gzip. See read_coefs

function parsePOSCAR(lines)

Parse a POSCAR from VASP

Called by makecrys, doesn't need to be called directly.

function parseQEinput(lines)

Parse a quantum espresso inputfile. Can only handle simple cases with explict CELL_PARAMETERS Cannot handle nonzero ibrav. or celldm

Called by makecrys, doesn't need to be called directly.

function generate_supercell(crys, cell)

Generate supercell. cell is [1,1,2], etc.

Note, perfered notation is to use syntax: c * [1,1,2] , where c is a crystal, thus using overloading of the * operator, rather than calling directly.

function generate_random_distortion(crys, amag, strain_mag)

Randomly distort a crystal. amag is the atom distance, strain_mag is the strain magnitude

function generate_optimum_supercell(c::crystal, dist)

Generate a supercell of crystal c where the periodic copies of all atoms are at least dist apart. Will consider (some) linear combinations of the initial lattice vectors and look for the "best" cell.

function write_poscar(crys, filename)

Write a crystal to a POSCAR.

function function write_efs(crys, energy, forces, stress, filename)

Write crystal, energy, force, stress to fake quantum espresso output file. This is for testing purposes only.

function get_grid(c, kden=55.0)

Get a default k-point grid size with kden density.

function orbital_index(c::crystal)

Get correspondence between crystal and the TB orbital numbers.

return ind2orb, orb2ind, etotal, nval

• ind2orb dictionary which gives [atom_number,atom_type, :orbital_symbol] from TB index integer.
• orb2ind dictionary which gives TB index integer from [atom_number,atom_type, :orbital_symbol]
• etotal total DFT energy of atoms, for calculation atomization energy.
• nval number of valence orbitals.

module DFToutMod

Deal with energy, force, stress, band structure, etc, from a DFT calculation

Generic to different DFT codes.

mutable struct bandstructure

Band structure. Has

• nbnd::Int Number of bands
• nks::Int Number of k-points
• nelec::Float64 Number of electrons
• nspin::Int64 Number of spins (2 for spin-polarized (magnetic), 1 non-sp)
• efermi::Float64 Fermi energy
• kpts::Array{Float64,2} List of k-points, nks × 3 array in BZ crystal units.
• kweights::Array{Float64,1} Weights of k-points (nks).
• kgrid::Array{Int,1} Equivalent gamma centered Monkhorst-Pack k-grid dimensions, if applies.
• eigs::Array{Float64,3} Eigenvalues. nks × nbnd × nspin

mutable struct dftout

DFT output struct. Has

• crys::crystal crystal structure
• energy::Float64 the actual DFT energy, depends on pseudopotentials.
• energy_smear::Float64 smearing energy
• forces::Array{Float64,2} forces Ryd / a.u.
• stress::Array{Float64,2} stress Ryd / (a.u.)^3
• bandstruct::bandstructure See bandstructure struct
• hasband::Bool does this object have a band structure, usually true
• hasham::Bool not used, always false
• prefix::String A string has a name used to find output files.
• outdir::String Directly loaded from
• tot_charge::Float64 If charge of unit cell is nonzero
• atomize_energy::Float64 Atomization energy, relative to non-spin-polarized atoms.
• nspin::Int64 number of spins (1 = non-sp, 2= spin-polarized)
• mag_tot::Float64 Total magnetization = sumi magi
• mag_abs::Float64 Absolute magnetization = sumi | magi |

function get_atomize_energy(d::dftout)

Return atomization energy in current energy units

function makebs(nelec::Number, efermi::Number,  kpoints, kweights,kgrid, vals)

Constructor for bandstructure.

function makedftout(A, pos, types, energy::Number,energy_smear::Number,  forces, stress, bandstruct=missing; prefix="PREFIX", outdir="TMPDIR", tot_charge=0.0)

function makedftout(crys::crystal, energy::Number, energy_smear::Number,  forces, stress, bandstruct=missing; prefix="PREFIX", outdir="TMPDIR", tot_charge=0.0)

Constructor for dftout. Usually called by function that loads DFT output files, not called directly.

module DFT

This is the generic DFT interface. Only QE is currently implemented however.

function runSCF(crys::crystal; inputstr=missing, prefix=missing, tmpdir="./", directory="./", functional="PBESOL", wannier=0, nprocs=1, code="QE", skip=false, calculation="scf", dofree="all", tot_charge = 0.0, smearing = missing, magnetic=false, cleanup=false, use_backup=false)

Workflow for generic DFT SCF calculation. code can only by "QE"

module QE

Module for running Quantum Espresso. Generic DFT version below.

function doclean(d)

Clean up wavefunctions in directory d

function loadXML(savedir)

Load a QE SCF DFT calculation into a dftout object.

function loadXML_bs(savedir::String)

Load bandstructure from QE xml file that was converted to a dict already

function loadXML_bs(savedir::String)

Load bandstructure from QE xml file.

function make_commands(nprocs=1)

Returns a dictionary with command lines to call external programs on the command line nprocs is the number of processors for parallel execution Needs to be changed for you specfic program locations and mpi commands (if any)

This needs to be edited to actually run QE yourself. Running wannier90 is optional, not part of current code.

function makedict(savedir)

Load a data-file-schema.xml as julia dictionary.

function runSCF(crys::crystal, inputstr=missing, prefix=missing, tmpdir="./", directory="./", functional="PBESOL", wannier=0, nprocs=1, skip=false, calculation="scf", dofree="all", tot_charge = 0.0, smearing = 0.01, magnetic=false, cleanup=false, use_backup=false)

Workflow for doing SCF DFT calculation on crys

Return dftout

function run_pwscf(inputstr, outputstr, nprocs=1, directory="./", use_backup=false)

Run the pw.x code from QE on inputstr

module Ewald

Module for electrostatic preperation.

function electrostatics_getgamma(crys::crystal;  kappa=missing, noU=false, onlyU=false, screening = 1.0)

Main function. Does Ewald calculation on crys. Returns gamma, which is used in scf calculation. This is only run once for a given tb_crys object and stored.

• kappa is the splitting parameter between real/k-space in Ewald calculation. Will estimate best one if not provided.
• noU=false for testing only
• onlyU=false for testing only
• screening=1.0 Not used. Purpose is to reduce U values for values < 1.

function estimate_best_kappa(A)

Estimate best value of kappa for Ewald sum. Shouldn't effect final value, only calculation speed. There is probably a better way to do this.

function getU(types)

Get values of U from Atomdata.

function get_onsite(crys::crystal, U::Array{Float64,1})

return Onsite terms, Coulumb U.

function k_space(crys::crystal, kappa, starting_size_kspace=2)

K-space Ewald sum.

function k_space(crys::crystal, kappa, starting_size_kspace=2)

K-space Ewald sum.

function real_space(crys::crystal, kappa::Float64, U::Array{Float64}, starting_size_rspace=2)

Real-space Ewald sum.

function do_fitting(list_of_tbcs; kpoints = missing,  atoms_to_fit=missing, fit_threebody=true, fit_threebody_onsite=true, do_plot = true)

Used for simple linear fitting of coefficients. Interface for more complicated fitting.

• list_of_tbcs - List of tbc_crys or tbccrysk` objects to fit to.
• dft_list - for kspace fitting, use
• fit_threebody=true - Fit threebody coefficients. Sometimes false for testing, but true for production.
• fit_threebody_onsite=true - Fit threebody onsite coefficients. See above.
• do_plot=false - show simple plot comparing coefficients to tbc reference.

function do_fitting_linear(list_of_tbcs; kpoints = missing, dft_list = missing,  fit_threebody=true, fit_threebody_onsite=true, do_plot = false, starting_database=missing, mode=:kspace, return_database=true, NLIM=100, refit_database=missing)

Linear fitting (not recursive). Used as starting point of recursive fitting.

Arguments.

• list_of_tbcs The main data to fit to. Consists of a list of tb_crys or tb_crys_kspace objects.
• kpoints = missing Kpoints to do fitting to in k-space. Usually, this is not used, see below.
• dft_list = missing List of dftout objects. Normally, we get symmetry reduced k-grids from these objects.
• fit_threebody=true Fit three body coefficients. Yes for production runs.
• fit_threebody_onsite=true Fit 3body onsite coefficients. Yes for production runs.
• do_plot = false Make a plot to assess fitting.
• starting_database=missing Use a database dict with some of the coefficents already fit and fixed.
• mode=:kspace Fit in either :kspace or :rspace. Can only use r-space if using only tbc_crys real-space objects to fit to. :kspace is normal.
• return_database=true Return the final database. For use when called by other functions.
• NLIM=100 Largest number of k-points per structure. Set to smaller numbers to make code go faster / reduce memory, but may be less accurate.
• refit_database=missing starting point for coefficients we are fitting. Usually not used, as it doesn't always speed things up in practice. Something may not work about this option.

function do_fitting_recursive(list_of_tbcs ; weights_list = missing, dft_list=missing, X_cv = missing, kpoints = [0 0 0; 0 0 0.5; 0 0.5 0.5; 0.5 0.5 0.5], starting_database = missing,  update_all = false, fit_threebody=true, fit_threebody_onsite=true, do_plot = false, energy_weight = missing, rs_weight=missing,ks_weight=missing, niters=50, lambda=0.0, leave_one_out=false, prepare_data = missing, RW_PARAM=0.0, NLIM = 100, refit_database = missing, start_small = false)

This is the primary function for fitting. Uses the self-consistent linear fitting algorithm.

Arguments

• list_of_tbcs List of tbc_crys or tbc_crys_kspace object to fit to.
• weights_list = missing relative weights of different tbc objects in fitting code.
• dft_list=missing List of dftout objects used to get symmetry-reduced kpoint lists / weights to use in kspace fitting.
• kpoints = [0 0 0; 0 0 0.5; 0 0.5 0.5; 0.5 0.5 0.5] alternate way of setting k-points, not usually used.
• starting_database = missing If using already fit coefficents for some of the atoms, from another calculation, include that database dict here.
• update_all = false If updateall is true, then we refit the starting coefficients from `startingdatabase. Normallyfalse` except for testing.
• fit_threebody=true Fit 3body intersite coefficents. true for production runs.
• fit_threebody_onsite=true Fit 3body onsite coefficents. true for production runs.
• do_plot = false make plot for the linear fit.
• energy_weight = missing Weighting for the total energy terms in the fit.
• rs_weight=missing Real space hamiltonian matrix els weighting. zero for pure k-space fit.
• ks_weight=missing K-space hamiltonian matrix els weighting. Set to zero to ignore hamiltonian matrix els and only fit to band structure.
• niters=50 Maximum number of iterations.
• lambda=0.0 If greater than zero, include a simple ridge regression with this lambda value. Usually zero.
• leave_one_out=false Leave-one-out cross-validation. Too slow to be very useful.
• prepare_data = missing Rarely used option to reuse previous linear fitting.
• RW_PARAM=0.0 Weighting of non-occupied bands in fit.
• NLIM = 100 Maximum k-points per structure. Smaller for faster but less accurate fit that uses less memory.
• refit_database = missing Option to include starting data. Rarely used.
• start_small = false When fitting only 3body data, setting this to true will start the 3body terms with very small values, which can improve convergence. Not useful if also fitting 2body terms.

function extract_database(database_old,nh,ns, KEYS, HIND, SIND)

If we are using fixed prefit coefficients, we have to get them in the same form as our current fitting.

function fourierspace(tbc, kpoints, X_H, X_S, Y_H, Y_S, Xhc, Xsc, rind, Rvec, INDVec, h_on, ind_convert)

Do analytic fourier transform of real space fitting matrices into kspace.

function get_k(dft_list, ncalc; NLIM = 100)

Decide which k-points to include in fitting, as we limit the total number to NLIM or less per structure

Uses some randomness, but puts high symmetry points at front of line.

function hermetian_index(i::Int64,j::Int64,nwan::Int64)

This is used to reduce memory by only keeping track of independet coefficients of Hermetian matrices, which is nearly a factor of 2 reduction.

function make_database(ch, cs,  KEYS, HIND, SIND, DMIN_TYPES, DMIN_TYPES3; scf=false, starting_database=missing, tbc_list=missing)

Construct the coefs and database from final results of fitting.

function prepare_for_fitting(list_of_tbcs; kpoints = missing, dft_list = missing, fit_threebody=false, fit_threebody_onsite=false, starting_database=missing, refit_database=missing)

Make lots of preperations for fitting. Moves things around, put stuff in materices, etc.

module Force_Stress

Module for calculating force and stress

function finite_diff(crys::crystal, database, ind1, ind2; stress_mode=false, step = 0.0002, smearing = 0.01, grid = missing)

Finite differences force/stress, for testing.

Arguments

• crys::crystal Crystal structure
• database Database of fitting coefficents.
• ind1 atom index for first stress index
• ind2 cartesian index or second stress index
• stress_mode=false true for stress, otherwise force.
• step = 0.0002 step_size for finite steps.
• smearing = 0.01 smearing energy
• grid = missing kpoint grid

function get_energy_force_stress(crys::crystal, database; smearing = 0.01, grid = missing)

Get force and stress, non-fft algorithm. Generally use the fft algorithm.

return energy_tot, f_cart, stress

Returns Ryd units. Generally users should use the scf_energy_force_stress function

Uses automatic differentation for gradient.

function get_energy_force_stress(crys::crystal, database; smearing = 0.01, grid = missing)

Get force and stress, non-fft algorithm

function get_energy_force_stress_fft(tbc::tb_crys, database; do_scf=false, smearing = 0.01, grid = missing, e_den0=missing, vv = missing)

Calculate energy/force/stress using fft algorithm. Users should use scf_energy_force_stress, which calls this. Uses automatic differentation for jacobian.

function get_energy_force_stress_fft(tbc::tb_crys, database; do_scf=false, smearing = 0.01, grid = missing, e_den0=missing, vv = missing)

Calculate energy/force/stress using fft algorithm. Users should use scf_energy_force_stress, which calls this. Uses automatic differentation for jacobian.

function get_energy_force_stress_fft(tbc::tb_crys, database; do_scf=false, smearing = 0.01, grid = missing, e_den0=missing, vv = missing)

Calculate energy/force/stress using fft algorithm. Users should use scf_energy_force_stress, which calls this. Uses automatic differentation for jacobian.

function get_energy_force_stress_fft_LV_sym_SINGLE(tbc::tb_crys_sparse)

Main sparse matrix force/stress calculator. Should be faster for large systems, due to faster matrix/vector and matrix/matrix multiply (assuming sparscity). Also requires much less memory.

function safe_mode_energy(crys::crystal, database; var_type=Float64)

Relaxation can accidently lead to very small atom-atom distances during the relaxation precedure if too large of step is taken. This function is a repulsive energy function at short range to make sure the relaxtion doesn't get stuck at very short distances where the fitting doesn't apply.

function sparsify_jac(jac, INDH, INDS, nwan)

A pain of using ForwardDiff.jacobian on a function that returns a sparse matrix Hamiltonian is that ForwardDiff wants a dense return. So we have to do a lot of work returning a dense matrix version of the sparse matrix and the index of the non-zero entries. We then sparsify it here.

module ManageDatabase

Module for reading coefs from files and making database as needed for calculations

function add_to_database(s::Set)

Load elements or twobody terms from precalcuated coefs from files.

function clear_database()

Remove loaded databasea

function prepare_database(at_list)

function prepare_database(c::crystal)

Get ready database of precalculated coefs for crystal

module SCF

Module for self-consistent field calculations for TB objects.

function remove_scf_from_tbc(hk3, sk3, tbc; smearing=0.01, e_den = missing)

This function takes a hk3, sk3 set of hamiltonian / overlap that does not require scf and adjusts it so that it does require scf, but gives the same energy and band structure.

function remove_scf_from_tbc(tbcK::tb_crys_kspace; smearing=0.01, e_den = missing)

This function takes a tbc_crys_kspace object that does not require scf and adjusts it so that it does require scf, but gives the same energy and band structure.

function remove_scf_from_tbc(tbc::tb_crys; smearing=0.01, grid = missing, e_den = missing)

This function takes a tbc_crys object that does not require scf and adjusts it so that it does require scf, but gives the same energy and band structure.

function scf_energy(tbc::tb_crys_sparse)

Sparse matrix implmentation of self-consistent field. Uses dense eigen routines, but sparse matrix mulitplication. Repeats a lot of code currently. TODO : refactor with less waste.

function scf_energy(c::crystal, database::Dict; smearing=0.01, grid = missing, conv_thr = 1e-5, iters = 75, mix = -1.0, mixing_mode=:pulay, verbose=true)

Run scf calculation of c::crystal, using database of coefs. The main user version is scf_energy in ThreeBodyTB, which calls this one.

• smearing is smearing energy in Ryd.
• grid is k-point grid (gamma centered MP), will use default.
• conv_thr convergence threshold in Ryd.
• iters maximum iterations for first attempt
• mix=-1.0 default is choose mixing for you. Otherwise, set between 0.0 and 1.0
• mixing_mode=:DIIS default using Pulay mixing (DIIS). Any other input uses simple mixing.
• verbose=true verbosity level.

function scf_energy(tbc::tb_crys; smearing=0.01, grid = missing, e_den0 = missing, conv_thr = 1e-5, iters = 100, mix = -1.0, mixing_mode=:pulay, verbose=true)

module Utility

Some useful functions, mostly for converting stuff and loading files and reshaping stuff.

function inv_reshape_vec(x, strain, nat; strain_mode=true)

The force and relax algorithms from outside codes take in vectors, not crystal 's. So we have to reshape vectors to and from crystals

function reshape_vec(x, nat; strain_mode=false)

The force and relax algorithms from outside codes take in vectors, not crystal 's. So we have to reshape vectors to and from crystals

function get_kpath_sym(c::crystal; sym_prec = 5e-4, magmoms = missing)

The group-theory based k-point path. Uses Spglib.jl for symmetry and largely follows the conventions of 'High-throughput electronic band structure calculations: Challenges and tools' - Setyawan and Curtarolo Comp Mater Sci 2010

Implementation based on jarvis-tools by C. Choudhary (https://github.com/usnistgov/jarvis/blob/master/jarvis/core/kpoints.py and https://github.com/usnistgov/jarvis/blob/master/jarvis/analysis/structure/spacegroup.py )

function get_standard_crys(c::crystal; sym_prec = 5e-4, magmoms=missing, to_primitive=true)

Return standardized crystal structure, per spglib convention. Check if space-group is correct.

function get_symmetry(c::crystal; verbose=true, sym_prec = 5e-4, magmoms=missing)

Return space group number, print other symmetry information if verbose=true. Based on Spglib.jl

function kassemble(kdict, names)

Helper function, converts names to k-points.

module Classical

Run a CLASSICAL model based on similar computational machinery to the TB3 model.

struct coefs_cl

Hold the classical coefficients for 2body or 3body interactions and EAM interactions.

function calc_energy_cl(crys::crystal

Main function for classical energy calculation from crystal structure.

Four body terms are experimental, may not be supported

function energy_force_stress_cl(crys::crystal)

Main function for running classical model. Returns energy/force/stress from crystal structure.

function make_coefs_cl(at_list, dim)

Constructor for coefs. Can create coefs filled with ones for testing purposes.

See coefs_cl to understand arguments.

function read_coefs_cl(filename, directory = missing)

Read Classical coefs_cl from filename. Can read gzipped files directly.

function write_coefs_cl(filename, co::coefs; compress=true)

Write coefs_cl to an file xml. Compress uses gzip. See read_coefs_cl

module CalcTB

Create TB matrix sets from coefficients, or prepare to fit the coefficients.

function do_fit_cl(DFT)

For fitting the classical model to a specific set of DFT calculations in the list DFT

function do_fit_cl(CRYS::Array{crystal,1})

Does the main fitting for CLASSICAL MODEL

function do_fit_cl_RECURSIVE(DFT_start)

Main CLASSICAL fitting routine. Runs DFT in recursive active-learning styel framework. For fitting classical energy model with basis of Laguerre polynomials times exponentials similar to the TB3 model. Can consider 2 and 3 body terms as well as an EAM (embedded atom model) style term. Experimental.

function kfold(L, n)

Helper function for kfold cross-validation