pynao package

Submodules

pynao.ao_log module

class pynao.ao_log.ao_log(**kw)

Bases: pynao.log_mesh.log_mesh

Holder of radial orbitals on logarithmic grid. Args:

ionslist of ion structures (read from ion files from siesta)

or

gtogaussian type of orbitals object from pySCF

or

gpaw : ??

Returns:

ao_log:

sp2ion (ion structure from m_siesta_ion or m_siesta_ion_xml):

List of structure composed of several field read from the ions file.

nr (int): number of radial point rmin (float) kmax (float) rmax (float) rr pp psi_log psi_log_rl sp2rcut (array, float): array containing the rcutoff of each specie sp_mu2rcut (array, float) interp_rr instance of log_interp_c to interpolate along real-space axis interp_pp instance of log_interp_c to interpolate along momentum-space axis

Examples:

ao_eval(ra, isp, coords)

comp_moments()

Computes the scalar and dipole moments of the product functions Args:

argument can be prod_log or ao_log

get_aoneo()

Packs the data into one array for a later transfer to the library

init_ao_log_gpaw(**kw)

Reads radial orbitals from a previous GPAW calculation.

init_ao_log_gto(**kw)

supposed to be private

init_ao_log_ion(**kw)

Reads data from a previous SIESTA calculation, interpolates the Pseudo-Atomic Orbitals on a single log mesh.

siesta_ion_interp(sp2ion, fname='paos')

view()

Shows a plot of all radial orbitals

pynao.bse_iter module

class pynao.bse_iter.bse_iter(**kw)

Bases: pynao.gw.gw, pynao.tddft_iter.tddft_iter

apply_kernel4p(ddm)

This applies the 4-point interaction kernel. This operator can be arbitrarily complex, therefore we formulate it as a procedure.

apply_kernel4p_nspin1(ddm)

apply_kernel4p_nspin2(ddm)

apply_l(sab, comega=0j)

This applies the interacting, two-point Green’s function to a suitable vector (e.g. dipole matrix elements)

apply_l0(sab, comega=0j)

This applies the non-interacting four point Green’s function to a suitable vector (e.g. dipole matrix elements)

apply_l0_ref(sab, comega=0j)

This applies the non-interacting four point Green’s function to a suitable vector (e.g. dipole matrix elements)

apply_l0_spin_non_diag(sab, comega=0j)

The definition of L0 which is not spin-diagonal even for parallel-spin references

comp_polariz_inter_ave(comegas)

Compute a direction-averaged interacting polarizability, i.e. <d |L| d>

comp_polariz_nonin_ave(comegas)

Non-interacting average polarizability, i.e. <d |L0| d>

define_e_x_l0(mo_energy, mo_coeff, **kw)

Define eigenvalues and eigenvecs for the two-particle Green’s function L0(1,2,3,4)

seff(sext, comega=0j)

This computes an effective two point field (scalar non-local potential) given an external two point field.

L = L0 (1 - K L0)^-1 We want therefore an effective X_eff for a given X_ext X_eff = (1 - K L0)^-1 X_ext or we need to solve linear equation (1 - K L0) X_eff = X_ext

The operator (1 - K L0) is named self.sext2seff_matvec

sext2seff_matvec(sab)

This is operator which we effectively want to have inverted (1 - K L0) and find the action of it’s inverse by solving a linear equation with a GMRES method. See the method seff(…)

pynao.chi0_matvec module

class pynao.chi0_matvec.chi0_matvec(**kw)

Bases: pynao.mf.mf

A class to organize the application of non-interacting response to a vector

apply_rf0(sp2v, comega, chi0_mv)

This applies the non-interacting response function to a vector (a set of vectors?)

calc_dens_edir_omega(iw, nww, w, edir, vext, tmp_fname=None, inter=False)

Calculate the density change and polarizability for a specific frequency

comp_dens_along_eext(comegas, eext=array([1.0, 0.0, 0.0]), tmp_fname=None, inter=False)

Compute a the average interacting polarizability along the eext direction for the frequencies comegas.

Input Parameters:

comegas (1D array, complex): the real part contains the frequencies at which the polarizability

should be computed. The imaginary part id the width of the polarizability define as self.eps

eext (1D xyz array, real): direction of the external field tmp_fname (string, optional): temporary file to where to store results

Particularly useful for large calculations

inter (bool, optional): Perform interactive calculations (or not)

Other Calculated quantity:

self.p_mat (complex array, dim: [3, 3, comega.size]): store the (3, 3) polarizability matrix

[[Pxx, Pxy, Pxz],

[Pyx, Pyy, Pyz], [Pzx, Pzy, Pzz]] for each frequency.

self.dn (complex array, dim: [3, comegas.size, self.nprod]): store the density change

comp_polariz_nonin_xx_atom_split(comegas)

Compute the non interacting polarizability along the xx direction

initialize_chi0_matvec_GPU()

Initialize GPU variable. Copy xocc, xvrt, ksn2e and ksn2f to the device

write_chi0_mv_timing(fname)

pynao.chi0_matvec.write_occupation(ksn2f, fname='orbitals_occupation.txt')

Write the orbital occupation to the file orbitals_occupation.txt

pynao.coords2sort_order module

pynao.coords2sort_order.coords2sort_order(a2c)

Delivers a list of atom indices which generates a near-diagonal overlap for a given set of atom coordinates

pynao.coords2sort_order.read_xyz(fname)

Reads xyz files

pynao.coords2sort_order.write_xyz(fname, s, ccc)

Writes xyz files

pynao.gw module

class pynao.gw.gw(**kw)

Bases: pynao.scf.scf

G0W0 with integration along imaginary axis

epsilon(ww)

This computes the dynamic dielectric function epsilon(r,r’.omega) = delta(r,r’) - v(r.r”) * chi_0(r”,r’,omega)

etot_gw()

g0w0_eigvals()

This computes the G0W0 corrections to the eigenvalues

get_h0_vh_x_expval()

This calculates the expectation value of Hartree-Fock Hamiltonian, when: self.get_hcore() = -1/2 del^{2}+ V_{ext} self.get_j() = Coloumb operator or Hartree energy vh self.get_k() = Exchange operator/energy mat1 is product of this Hamiltonian and molecular coefficients and it will be diagonalized in expval by einsum

get_snmw2sf(optimize='greedy')

This computes a matrix elements of W_c: <PsiPsi | W_c |PsiPsi>. sf[spin,n,m,w] = X^n V_mu X^m W_mu_nu X^n V_nu X^m, where n runs from s…f, m runs from 0…norbs, w runs from 0…nff_ia, spin=0…1 or 2.

get_wmin_wmax_tmax_ia_def(tol)

This is a default choice of the wmin and wmax parameters for a log grid along imaginary axis. The default choice is based on the eigenvalues.

gw_corr_int(sn2w, eps=None)

This computes an integral part of the GW correction at energies sn2e[spin,len(self.nn)]

rac{1}{2pi}int_{-infty}^{+infty }

sum_m

rac{I^{nm}(iomega{‘})}{E_n + iomega{‘}-E_m^0} domega{‘}

see eq (16) within DOI: 10.1021/acs.jctc.9b00436

gw_corr_res(sn2w)

This computes a residue part of the GW correction at energies sn2w[spin,len(self.nn)]

kernel_gw()

This creates the fields mo_energy_g0w0, and mo_coeff_g0w0

lsofs_inside_contour(ee, w, eps)

Computes number of states the eigenenergies of which are located inside an integration contour. The integration contour depends on w z_n = E^0_n-omega pm ieta

make_mo_g0w0()

This creates the fields mo_energy_g0w0, and mo_coeff_g0w0

report()

Prints the energy levels of meanfield and G0W0

report_ex()

Prints the self energy components in HF levels

report_mf(dm=None)

Prints the energy levels of mean-field hamiltonian

report_t()

Prints spent time in each step of GW class

rf(ww)

Full matrix interacting response from NAO GW class

rf0(ww)

Full matrix response in the basis of atom-centered product functions for parallel spins.

Blas version to speed up matrix matrix multiplication spped up of 7.237 compared to einsum version for C20 system

rf0_cmplx_ref(ww)

Full matrix response in the basis of atom-centered product functions

rf0_cmplx_ref_blk(ww)

Full matrix response in the basis of atom-centered product functions

rf0_cmplx_vertex_ac(ww)

Full matrix response in the basis of atom-centered product functions

rf0_cmplx_vertex_dp(ww)

Full matrix response in the basis of atom-centered product functions

rf_pyscf(ww)

Full matrix interacting response from tdscf class

si_c(ww, use_numba_impl=False)

This computes the correlation part of the screened interaction W_c by solving <self.nprod> linear equations (1-K chi0) W = K chi0 K or v_{ind}sim W_{c} = (1-vchi_{0})^{-1}vchi_{0}v scr_inter[w,p,q], where w in ww, p and q in 0..self.nprod

si_c_via_diagrpa(ww)

This method computes the correlation part of the screened interaction W_c via the interacting response function. The interacting response function, in turn, is computed by diagonalizing RPA Hamiltonian.

spin_square()

write_data(step=None)

writes data in RESTART’hdf5’ format

pynao.gw_iter module

class pynao.gw_iter.gw_iter(**kw)

Bases: pynao.gw.gw

gw_iter object.

Iterative G0W0 with integration along imaginary axis.

References:

1. Koval P. and al., Toward Efficient GW Calculations Using Numerical Atomic Orbitals: Benchmarking and Application to Molecular Dynamics Simulations, J. Chem. Theory Comput. 2019, 15, 8

check_veff(optimize='greedy')

This checks the equality of effective field (scalar potential) given the external scalar potential obtained from lgmres(linearopt, v_ext) and np.solve(dense matrix, vext).

chi0_mv(dvin, comega)

dupl_E(spin, thrs=None)

This returns index of states whose eigenenergy has difference less than thrs with former state

g0w0_eigvals_iter()

This computes the G0W0 corrections to the eigenvalues

get_snmw2sf_iter(nbnd=None, optimize='greedy')

This computes a matrix elements of W_c:

<Psi(r)Psi(r) | W_c(r,r’,omega) |Psi(r’)Psi(r’)>. sf[spin,n,m,w] = X^n V_mu X^m W_mu_nu X^n V_nu X^m,

where n runs from s…f,

m runs from 0…norbs, w runs from 0…nff_ia, spin=0…1 or 2.

1- XVX is calculated using dominant product: gw_xvx(‘dp_coo’) 2- I_nm = W XVX = (1-vchi_0)^{-1}vchi_0v 3- S_nm = XVX W XVX = XVX * I_nm

gw_applykernel(dn)

Apply the kernel to chi0_mv, i.e, matrix vector multiplication. The kernel is stored in pack or pack sparse format.

Input parameters: dn: 1D np.array, complex

the results of chi0_mv product, i.e, chi0*delta_V_ext

Output parameters: vcre: 1D np.array, float

The real part of the resulting matvec product kernel.dot(chi0_mv)

vcim: 1D np.array, float

The imaginary part of the resulting matvec product kernel.dot(chi0_mv)

gw_comp_veff(vext, comega=0j)

This computes an effective field (scalar potential) given the external scalar potential as follows:

(1-vchi_{0})V_{eff} = V_{ext} = X_{a}^{n}V_{mu}^{ab}X_{b}^{m} *

vchi_{0}v * X_{a}^{n}V_{nu}^{ab}X_{b}^{m}

returns V_{eff} as list for all n states(self.nn[s]).

gw_corr_int_iter(sn2w, eps=None)

This computes an integral part of the GW correction at GW class while uses get_snmw2sf_iter

gw_corr_res_iter(sn2w)

This computes a residue part of the GW correction at energies in iterative procedure

gw_vext2veffmatvec(vin)

gw_vext2veffmatvec2(vin)

gw_xvx(algo=None)

calculates basis products Psi(r’)Psi(r’) = XVX[spin,(nn, norbs, nprod)] = X_{a}^{n}V_{

u}^{ab}X_{b}^{m}

using 4-methods:

1- direct multiplication by using np.dot and np.einsum via swapping between axis 2- using atom-centered product basis 3- using atom-centered product basis and BLAS multiplication 4- using sparse atom-centered product basis and BLAS multiplication 5- using dominant product basis 6- using sparse dominant product basis

kernel_gw_iter()

This creates the fields mo_energy_g0w0, and mo_coeff_g0w0

make_mo_g0w0_iter()

This creates the fields mo_energy_g0w0, and mo_coeff_g0w0

precond_lgmres(omega=None)

For a simple V_eff, this calculates linear operator of M = (1-chi_0(omega) = kc_opt)^-1, as a preconditioner for lgmres

si_c_check(tol=1e-05)

This compares np.solve and LinearOpt-lgmres methods for solving linear equation (1-vchi_{0}) * W_c = vchi_{0}v

si_c_iter(ww)

This computes the correlation part of the screened interaction using LinearOpt and lgmres lgmres method is much slower than np.linalg.solve !!

pynao.gw_iter.profile(fnc)

Profiles any function in following class just by adding @profile above function

pynao.log_mesh module

pynao.log_mesh.funct_log_mesh(nr, rmin, rmax, kmax=None)

Initializes log grid in real and reciprocal (momentum) spaces. These grids are used in James Talman’s subroutines.

pynao.log_mesh.get_default_log_mesh_param4gpaw(sp2dic)

Determines the default (optimal) parameters for radial orbitals given on equidistant grid

pynao.log_mesh.get_default_log_mesh_param4gto(gto, tol_in=None)

pynao.log_mesh.get_default_log_mesh_param4ion(sp2ion)

class pynao.log_mesh.log_mesh(**kw)

Bases: object

Constructor of the log grid used with NAOs.

init_log_mesh(**kw)

Taking over the given grid rr and pp

init_log_mesh_gpaw(**kw)

This initializes an optimal logarithmic mesh based on setups from GPAW

init_log_mesh_gto(**kw)

Initialize an optimal logarithmic mesh based on Gaussian orbitals from pySCF

init_log_mesh_ion(**kw)

Initialize an optimal logarithmic mesh based on information from SIESTA ion files

pynao.lsofcsr module

pynao.lsofcsr.lsofcsr(coo3, dtype=<class 'float'>, shape=None, axis=0)

Generate a list of csr matrices out of a 3-dimensional coo format Args:

coo3 : must be a tuple (data, (i1,i2,i3)) in analogy to the tuple (data, (rows,cols)) for a common coo format shape : a tuple of dimensions if they are known or cannot be guessed correctly from the data axis : index (0,1 or 2) along which to construct the list of sparse arrays

Returns:

list of csr matrices

class pynao.lsofcsr.lsofcsr_c(coo3, dtype=<class 'float'>, shape=None, axis=0)

Bases: object

toarray()

pynao.m_ao_eval module

pynao.m_ao_eval.ao_eval(ao, ra, isp, coords)

Compute the values of atomic orbitals on given grid points Args:

ao : instance of ao_log_c class ra : vector where the atomic orbitals from “ao” are centered isp : specie index for which we compute coords: coordinates on which we compute

Returns:

res[norbs,ncoord] : array of atomic orbital values

pynao.m_ao_eval_libnao module

pynao.m_ao_eval_libnao.ao_eval_libnao(ao, ra, isp, coords)

pynao.m_ao_eval_libnao.ao_eval_libnao_(ao, rat, isp, crds, res)

Compute the values of atomic orbitals on given grid points Args:

ao : instance of ao_log_c class rat : vector where the atomic orbitals from “ao” are centered isp : specie index for which we compute crds: coordinates on which we compute

Returns:

res[norbs,ncoord] : array of atomic orbital values

pynao.m_ao_log_hartree module

pynao.m_ao_log_hartree.ao_log_hartree(ao, ao_log_hartree_method='lap')

pynao.m_ao_log_hartree.ao_log_hartree_lap(ao)

Computes radial parts of Hartree potentials generated by the radial orbitals using Laplace transform (r>, r<)

self: class instance of ao_log_c

ao_pot with respective radial parts

pynao.m_ao_log_hartree.ao_log_hartree_lap_libnao(ao)

Computes radial parts of Hartree potentials generated by the radial orbitals using Laplace transform (r>, r<) Uses a call to Fortran version of the calculation.

self: class instance of ao_log_c

ao_pot with respective radial parts

pynao.m_ao_log_hartree.ao_log_hartree_sbt(ao)

Computes radial parts of Hartree potentials generated by the radial orbitals using the spherical Bessel transform

self: class instance of ao_log

ao_pot with respective radial parts

pynao.m_ao_matelem module

class pynao.m_ao_matelem.ao_matelem_c(rr, pp, sv=None, dm=None)

Bases: pynao.m_sbt.sbt_c, pynao.m_c2r.c2r_c, pynao.m_gaunt.gaunt_c

Evaluator of matrix elements given by the numerical atomic orbitals. The class will contain

the Gaunt coefficients, the complex -> real transform (for spherical harmonics) and the spherical Bessel transform.

coulomb_am(sp1, R1, sp2, R2)

coulomb_ni(sp1, R1, sp2, R2, **kvargs)

init_one_set(ao, **kvargs)

Constructor for two-center matrix elements, i.e. one set of radial orbitals per specie is provided

init_two_sets(ao1, ao2, **kvargs)

Constructor for matrix elements between product functions and orbital’s products: two sets of radial orbitals must be provided.

overlap_am(sp1, R1, sp2, R2)

overlap_ni(sp1, R1, sp2, R2, **kvargs)

xc_scalar(sp1, R1, sp2, R2, **kvargs)

pynao.m_ao_matelem.build_3dgrid(me, sp1, R1, sp2, R2, **kw)

pynao.m_ao_matelem.build_3dgrid3c(me, sp1, sp2, R1, R2, sp3, R3, level=3)

pynao.m_aos_libnao module

pynao.m_aos_libnao.aos_libnao(coords, norbs)

pynao.m_c2r module

class pynao.m_c2r.c2r_c(j)

Bases: object

Conversion from complex to real harmonics

c2r_(j1, j2, jm, cmat, rmat, mat)

c2r_moo(j, mab_c, mu2info)

Transform tensor m, orb, orb given in complex spherical harmonics to real spherical harmonic

c2r_vector(clm, l, si, fi)

pynao.m_chi0_noxv module

pynao.m_chi0_noxv.calc_ab2v(mat1, mat2, matin)

calculate ab2v via two matrix-matrix multiplication in an row

pynao.m_chi0_noxv.calc_nm2v(mat1, mat2, matin)

calculate nm2v via two matrix-matrix multiplication in an row

pynao.m_chi0_noxv.calc_sab(mat1, mat2, vec, GPU=False)

Perform two matrix-vector multiplication in an row, i.e

out = M2.dot(M1.dot(in))

pynao.m_chi0_noxv.chi0_mv(self, dvin, comega, dnout=None)

Apply the non-interacting response function Chi_0 to the induced effetive potential delta V_eff. See equation 2.81 and 2.83in Ref 1.

Input Parameters:

self : tddft_iter or tddft_tem class sp2v : vector describing the effective perturbation [spin*product] –> value comega: complex frequency

Reference:

[1]: Plasmon in Nanoparticles: Atomistic Ab initio theory for large Systems, M. Barbry, PhD thesis, 2018

pynao.m_chi0_noxv.chi0_mv_gpu(self, dvin, comega, dnout=None)

Apply the non-interacting response function Chi_0 to the induced effetive potential delta V_eff. See equation 2.81 and 2.83in Ref 1. Performs the matrix-matrix operation on GPU

Input Parameters:

self : tddft_iter or tddft_tem class sp2v : vector describing the effective perturbation [spin*product] –> value comega: complex frequency

Reference:

[1]: Plasmon in Nanoparticles: Atomistic Ab initio theory for large Systems, M. Barbry, PhD thesis, 2018

pynao.m_chi0_noxv.div_eigenenergy(ksn2e, ksn2f, vstart, nfermi, comega, spin, use_numba, nm2v_re, nm2v_im, div_numba=None, GPU=False, blockspergrid=None, threadsperblock=None)

Divide by the energy terms

omega - (E_m - E_n) + i*varepsilon

See equation 2.81 in Ref 1.

Reference:

[1]: Plasmon in Nanoparticles: Atomistic Ab initio theory for large Systems, M. Barbry, PhD thesis, 2018

pynao.m_color module

class pynao.m_color.color

Bases: object

BLUE = '\x1b[94m'

ENDC = '\x1b[0m'

FAIL = '\x1b[91m'

GREEN = '\x1b[92m'

HEADER = '\x1b[95m'

RED = '\x1b[91m'

T = None

WARNING = '\x1b[93m'

disable()

os = <module 'os' from '/home/marc/anaconda3/envs/pynao/lib/python3.7/os.py'>

pynao.m_comp_coulomb_den module

pynao.m_comp_coulomb_den.comp_coulomb_den(sv, ao_log=None, funct=<function coulomb_am>, dtype=<class 'numpy.float64'>, **kvargs)

Computes the matrix elements given by funct, for instance coulomb interaction

Input parameters:

sv: System Variables

this must have arrays of coordinates and species, etc

Output parameters:

res: 2D np.array, dtype

matrix elements (real-space overlap) for the whole system

pynao.m_comp_coulomb_pack module

pynao.m_comp_coulomb_pack.comp_coulomb_pack(sv, ao_log=None, funct=<function coulomb_am>, dtype=<class 'numpy.float64'>, **kw)

Computes the matrix elements given by funct, for instance coulomb interaction

Input parameters: sv: (System Variables), this must have arrays of coordinates and species, etc ao_log : description of functions (either orbitals or product basis functions)

Output parameters:

res: 1D np.array dtype

matrix elements for the whole system in packed form (lower triangular part)

norbs: int

Dimension of the matrix (norbs, norbs)

pynao.m_comp_dm module

pynao.m_comp_dm.comp_dm(ksnac2x, ksn2occ)

Computes the density matrix Args:

ksnac2x : eigenvectors ksn2occ : occupations

Returns:

ksabc2dm :

pynao.m_comp_spatial_distributions module

class pynao.m_comp_spatial_distributions.spatial_distribution(dn, freq, box, excitation='light', **kw)

Bases: pynao.mf.mf

class to calculate spatial distribution of

• density change

• induce potential

• induce electric field

• intensity of the induce Efield

Example:

from __future__ import print_function, division import numpy as np from pynao import tddft_iter from pynao.m_comp_spatial_distributions import spatial_distribution

from ase.units import Ry, eV, Ha, Bohr

# run tddft calculation td = tddft_iter(label=”siesta”, iter_broadening=0.15/Ha, xc_code=’LDA,PZ’)

omegas = np.linspace(0.0, 10.0, 200)/Ha + 1j*td.eps td.comp_dens_inter_along_Eext(omegas, Eext=np.array([1.0, 0.0, 0.0]))

box = np.array([[-10.0, 10.0],

[-10.0, 10.0], [-10.0, 10.0]])/Bohr

dr = np.array([0.5, 0.5, 0.5])/Bohr

# initialize spatial calculations spd = spatial_distribution(td.dn, omegas, box, dr = dr, label=”siesta”)

# compute spatial density change distribution spd.get_spatial_density(3.5/Ha, Eext=np.array([1.0, 0.0, 0.0]))

# compute Efield Efield = spd.comp_induce_field()

# compute potential pot = spd.comp_induce_potential()

# compute intensity intensity = spd.comp_intensity_Efield(Efield)

comp_induce_field()

Compute the induce Electric field corresponding to the density change calculated in get_spatial_density

comp_induce_potential()

Compute the induce potential corresponding to the density change calculated in get_spatial_density

comp_intensity_Efield(Efield)

Compute the intesity of the induce electric field

get_spatial_density(w0, ao_log=None, Eext=array([1.0, 1.0, 1.0]))

Compute the density change fromm the product basis to cartesian bais for the frequency w0 and the external field directed along Eext

pynao.m_comp_spatial_numba module

pynao.m_comp_spatial_numba.get_spatial_density_numba(dn_spatial, mu2dn, meshx, meshy, meshz, atom2sp, atom2coord, atom2s, sp_mu2j, psi_log_rl, sp_mu2s, sp2rcut, rr, res, coeffs, gammin_jt, dg_jt, nr)

pynao.m_comp_spatial_numba.get_spatial_density_numba_parallel(dn_spatial, mu2dn, meshx, meshy, meshz, atom2sp, atom2coord, atom2s, sp_mu2j, psi_log_rl, sp_mu2s, sp2rcut, rr, res, coeffs, gammin_jt, dg_jt, nr)

pynao.m_comp_vext_tem module

pynao.m_comp_vext_tem.comp_vext_tem(nprod, velec, beam_offset, freq, freq_sym, time_range, ao_log=None)

Compute the external potential created by a moving charge using the fortran routine

pynao.m_comp_vext_tem.comp_vext_tem_pyth(self, ao_log=None, numba_parallel=True)

Compute the external potential created by a moving charge Python version

pynao.m_comp_vext_tem_numba module

pynao.m_comp_vext_tem_numba.c2r_lm(conj_c2r, jmx, clm, clmm, m)

clm: sph harmonic l and m clmm: sph harmonic l and -m convert from real to complex spherical harmonic for an unique value of l and m

pynao.m_comp_vext_tem_numba.get_index_lm(l, m)

return the index of an array ordered as [l=0 m=0, l=1 m=-1, l=1 m=0, l=1 m=1, ….]

pynao.m_comp_vext_tem_numba.get_tem_potential_numba(time, R0, vnorm, vdir, center, rcut, inte1, rr, dr, fr_val, conj_c2r, l, m, jmx, ind_lm, ind_lmm, V_time)

Numba version of the computation of the external potential in time for tem calculations

pynao.m_conv_ac_dp module

class pynao.m_conv_ac_dp.conv_ac_dp_c(pb, dtype=<class 'numpy.float64'>)

Bases: object

pynao.m_conv_yzx2xyz module

class pynao.m_conv_yzx2xyz.conv_yzx2xyz_c(mol)

Bases: object

A class to organize a permutation of l=1 matrix elements

conv_yzx2xyz_1d(mat_yzx, m_yzx2m_xyz, sh_slice=None)

Permutation of first index in a possibly multi-dimensional array. This would convert nao->pyscf convention

conv_yzx2xyz_2d(mat_yzx, direction='pynao', ss=None)

conv_yzx2xyz_4d(mat_yzx, direction='pynao', ss=None)

pynao.m_coulomb_am module

pynao.m_coulomb_am.coulomb_am(self, sp1, R1, sp2, R2, **kvargs)

Computes Coulomb overlap for an atom pair. The atom pair is given by a pair of species indices and the coordinates of the atoms. <a|r^-1|b> = iint a(r)|r-r’|b(r’) dr dr’ Args:

self: class instance of ao_matelem_c sp1,sp2 : specie indices, and R1,R2 : respective coordinates

Result:

matrix of Coulomb overlaps

The procedure uses the angular momentum algebra and spherical Bessel transform. It is almost a repetition of bilocal overlaps.

pynao.m_csphar module

pynao.m_csphar.csphar(r, lmax)

Computes (all) complex spherical harmonics up to the angular momentum lmax

Args:

r : Cartesian coordinates defining correct theta and phi angles for spherical harmonic lmax : Integer, maximal angular momentum

Result:

1-d numpy array of complex128 elements with all spherical harmonics stored in order 0,0; 1,-1; 1,0; 1,+1 … lmax,lmax, althogether 0 : (lmax+1)**2 elements.

pynao.m_csphar_talman_libnao module

pynao.m_csphar_talman_libnao.csphar_talman_libnao(rvec, lmax)

pynao.m_csphar_talman_libnao.talman2world(ylm_t)

pynao.m_dens_elec_vec module

pynao.m_dens_elec_vec.dens_elec_vec(sv, crds, dm)

Compute the electronic density using vectorized oprators

pynao.m_dens_libnao module

pynao.m_dens_libnao.dens_libnao(crds, nspin)

Compute the electronic density using library call

pynao.m_density_cart module

pynao.m_density_cart.density_cart(crds)

Compute the values of atomic orbitals on given grid points Args:

sv : instance of system_vars_c class crds : vector where the atomic orbitals from “ao” are centered sab2dm : density matrix

Returns:

res[ncoord] : array of density

pynao.m_dipole_coo module

pynao.m_dipole_coo.dipole_coo(sv, ao_log=None, funct=<function dipole_ni>, **kvargs)

Computes the dipole matrix and returns it in coo format (simplest sparse format to construct) Args:

sv : (System Variables), this must have arrays of coordinates and species, etc

Returns:

overlap (real-space overlap) for the whole system

pynao.m_dipole_ni module

pynao.m_dipole_ni.dipole_ni(me, sp1, R1, sp2, R2, **kvargs)

Computes overlap for an atom pair. The atom pair is given by a pair of species indices and the coordinates of the atoms. Args:

sp1,sp2 : specie indices, and R1,R2 : respective coordinates in Bohr, atomic units

Result:

matrix of orbital overlaps

The procedure uses the numerical integration in coordinate space.

pynao.m_div_eigenenergy_numba module

pynao.m_div_eigenenergy_numba.div_eigenenergy_numba(n2e, n2f, nfermi, vstart, comega, nm2v_re, nm2v_im)

multiply the temporary matrix by (fn - fm) (frac{1.0}{w - (Em-En) -1} -

frac{1.0}{w + (Em - En)})

using numba

pynao.m_div_eigenenergy_numba.mat_mul_numba(a, b)

pynao.m_div_eigenenergy_numba_gpu module

pynao.m_dm module

pynao.m_dm.dm(ksnac2x, ksn2occ)

Computes the density matrix Args:

ksnar2x : eigenvectors ksn2occ : occupations

Returns:

ksabc2dm :

pynao.m_dos_pdos_eigenvalues module

pynao.m_dos_pdos_eigenvalues.eigen_dos(ksn2e, zomegas, nkpoints=1)

Compute the Density of States using the eigenvalues

pynao.m_dos_pdos_eigenvalues.eigen_pdos(ksn2e, zomegas, nkpoints=1)

Compute the Partial Density of States using the eigenvalues

pynao.m_dos_pdos_eigenvalues.plot(d_qp=None)

This plot DOS and PDOS for spin-polarized calculations in both mean-field and GW levels

pynao.m_dos_pdos_eigenvalues.read_qp_molgw(filename)

reading QP energies from MOLGW output for DFT calculations

pynao.m_dos_pdos_ldos module

pynao.m_dos_pdos_ldos.gdos(mf, zomegas, omask=None, mat=None, nkpoints=1)

Compute some masked (over atomic orbitals) or total Density of States or any population analysis

pynao.m_dos_pdos_ldos.lsoa_dos(mf, zomegas, lsoa=None, nkpoints=1)

Compute the Partial Density of States according to a list of atoms

pynao.m_dos_pdos_ldos.omask2wgts_loops(mf, omask, over)

Finds weights

pynao.m_dos_pdos_ldos.pdos(mf, zomegas, nkpoints=1)

Compute the Partial Density of States (resolved in angular momentum of the orbitals) using the eigenvalues and eigenvectors in wfsx

pynao.m_eri2c module

pynao.m_eri2c.eri2c(me, sp1, R1, sp2, R2, **kvargs)

Computes two-center Electron Repulsion Integrals. This is normally done between product basis functions

pynao.m_eri3c module

pynao.m_eri3c.eri3c(me, sp1, sp2, R1, R2, sp3, R3, **kvargs)

Computes three-center Electron Repulsion Integrals. atomic orbitals of second species must contain the Hartree potentials of radial functions (product orbitals)

pynao.m_exc module

pynao.m_exc.exc(sv, dm, xc_code, **kvargs)

Computes the exchange-correlation energy for a given density matrix Args:

sv : (System Variables), this must have arrays of coordinates and species, etc xc_code : is a string must comply with pySCF’s convention PZ

“LDA,PZ” “0.8*LDA+0.2*B88,PZ”

Returns:

exc x+c energy

pynao.m_fact module

pynao.m_fermi_dirac module

pynao.m_fermi_dirac.fermi_dirac(telec, e, fermi_energy)

Fermi-Dirac distribution function

pynao.m_fermi_dirac.fermi_dirac_occupations(telec, ksn2e, fermi_energy)

Occupations according to the Fermi-Dirac distribution function.

pynao.m_fermi_energy module

pynao.m_fermi_energy.fermi_energy(ee, nelec, telec)

Determine Fermi energy

pynao.m_fermi_energy.func1(x, i2e, nelec, telec, norm_const)

pynao.m_fft module

Module to that correct the Discret Fourier Transform in order to get the analytical value.

pynao.m_fft.FT1(t, f, axis=0, norm='asym', d=0)

Calculate the FFT of a ND dimensionnal array over the axes axes corresponding to the FT by using numpy fft. FT can include additional damping term to account for iterative broadening. Two normalizations are possible: symmetric (“sym” - with sqrt(2pi) factor in both FT and iFT and asymmetric (“asym”) with 1/2pi factor in iFT. The Fourier transform is done on the first dimension.

Input:

t, ND numpy array containing the time variable f, ND numpy array containing the data to Fourier transform (f(t, x, ….)) norm, string determining the utilized normalization d, float indicating the amount of iterative broadening

output:

F, ND numpy of the FT along axis

pynao.m_fft.FT2(x, y, f)

2D Fourier Transform

Input Parameters:

x, y (1D array): axis range f (2D numpy array of shape (x.size, y.size)): array to FT

Output Parameters:

F (2D numpy array of shape (x.size, y.size)): FT

Example:

import numpy as np import pynao.m_fft as fft import matplotlib.pyplot as plt

def gauss2D(x, y, a, b):

return np.exp(-(a*x**2 + b*y**2))

x = np.linspace(-5.0, 5.0, 100) y = np.linspace(-7.0, 7.0, 200)

a = 1.0 b = 2.0

xv, yv = np.meshgrid(y, x) f = gauss2D(xv, yv, a, b)

f_FT = fft.FT2(x, y, f)

kx = fft.get_fft_freq(x) ky = fft.get_fft_freq(y) if_FT = fft.iFT2(kx, ky, f_FT)

plt.imshow(f) plt.colorbar() plt.show()

plt.imshow(if_FT.real) plt.colorbar() plt.show()

pynao.m_fft.FT3(x, y, z, f)

3D Fourier Transform

Input Parameters:

x, y, z (1D array): axis range f (3D numpy array of shape (x.size, y.size, z.size)): array to FT

Output Parameters:

F (3D numpy array of shape (x.size, y.size, z.size)): FT

pynao.m_fft.FTconvolve(f, g, *args, domain='time')

Perform a corrected version of the fft convolution from scipy

Input arguments:

f (1, 2 or 3D array): first term of the product g (1, 2 or 3D array): second term of the product (same dim than f) args (list of 1D array): contains the real space array variable,

it must contains at least 1 argument and in maximum 3. the sumber of arguments must match the dimemension of f and g.

domain (string): domain in which the inputs are provided (“time” or

“frequency”)

Output arguments:

conv (1, 2, or 3D array): convolution product

pynao.m_fft.FTderivative(f, *args, domain='time')

Calculate derivative using FFT

Input arguments:

f (1, 2 or 3D array): the array to be differentiated args (list of 1D array): contains the real space array variable,

it must contains at least 1 argument and in maximum 1. the number of arguments must match the dimemension of f.

Output arguments:

deriv (1D array): derivative of f

pynao.m_fft.FTextend(omegas, f)

Computes Fourier transform in the entire domain based on the positive frequency part of FT under assumption that the transformed function is real valued

Input arguments:

omegas (1, 2 or 3D array): positive frequencies f (1, 2 or 3D array): Fourier transform args (list of 1D array): contains the real space array variable,

it must contains at least 1 argument and in maximum 3. the sumber of arguments must match the dimemension of f and g.

Output arguments:

omegatot (1D array): frequencies spectot (1D array): full Fourier transform

pynao.m_fft.calc_postfactor_2D(kx, ky, tmin, F, sign=1.0)

pynao.m_fft.calc_postfactor_3D(kx, ky, kz, tmin, F, sign=1.0)

pynao.m_fft.calc_prefactor_2D(x, y, f, wmin, tmin, dt, sign=1.0)

pynao.m_fft.calc_prefactor_3D(x, y, z, f, wmin, tmin, dt, sign=1.0)

pynao.m_fft.get_fft_freq(t)

return the frequency range of the fft variable, it is a well define array,

dw = 2*pi/(N*dt)

with N the size of the array

pynao.m_fft.iFT1(w, F, axis=0, norm='asym', d=0)

Calculate the FFT of a 1D dimensionnal array corresponding to the FT by using numpy fft iFT can include additional damping term to account for iterative broadening. Two normalizations are possible: symmetric (“sym” - with sqrt(2pi) factor in both FT and iFT and asymmetric (“asym”) with 1/2pi factor in iFT.

w, 1D numpy array containing the frequency variable from fft_freq F, 1D numpy array containing the data to inverse Fourier transform norm, string determining the utilized normalization d, float indicating the amount of iterative broadening

F, 1D numpy array containing the iFT

pynao.m_fft.iFT2(kx, ky, F)

2D Inverse Fourier Transform

Input Parameters:

kx, ky (1D array): axis range F (2D numpy array of shape (x.size, y.size)): array to IFT

Output Parameters:

f (2D numpy array of shape (x.size, y.size)): IFT

pynao.m_fft.iFT3(kx, ky, kz, F)

3D Inverse Fourier Transform

Input Parameters:

kx, ky, kz (1D array): axis range F (3D numpy array of shape (kx.size, ky.size, kz.size)): array to IFT

Output Parameters:

f (3D numpy array of shape (kx.size, ky.size, kz.size)): IFT

pynao.m_fireball_get_HS_dat module

pynao.m_fireball_get_HS_dat.fireball_get_HS_dat(cd, fname='HS.dat')

pynao.m_fireball_get_cdcoeffs_dat module

pynao.m_fireball_get_cdcoeffs_dat.fireball_get_cdcoeffs_dat(cd)

the file does not seems to be complete (either number of k points or number of spins is missing)

pynao.m_fireball_get_cdcoeffs_dat.fireball_get_ucell_cdcoeffs_dat(cd)

the file does not seems to be complete (either number of k points or number of spins is missing)

pynao.m_fireball_get_eigen_dat module

pynao.m_fireball_get_eigen_dat.fireball_get_eigen_dat(cd)

Calls

pynao.m_fireball_hsx module

class pynao.m_fireball_hsx.fireball_hsx(nao, **kw)

Bases: object

deallocate()

pynao.m_fireball_import module

pynao.m_fireball_import.fireball_import(self, **kw)

Calls

pynao.m_g297 module

The following contains a database of 148 small-sized molecules belonging to the G2/97 database Raghavachari, Redfern, and Pople, J. Chem. Phys. Vol. 106, 1063 (1997). See http://www.cse.anl.gov/Catalysis_and_Energy_Conversion/Computational_Thermochemistry.shtml for the original files. Or https://github.com/qsnake/ase/blob/master/ase/data/

pynao.m_g297.mol_cas(casno)

pynao.m_g297.mol_name(in_name)

pynao.m_g297.molgw_input(dirname=None)

writes a standard MOLGW inputs for HF+G0W0 in the given address

pynao.m_g297.pos_latex()

prepares latex file for open-shell systems, number of atoms and their positions

pynao.m_g297.pos_pyscf()

converts ase position lists to a long string which is readable for Pyscf

pynao.m_g297.pyscf_input(dirname=None)

writes the PySCF inputs in the given address

pynao.m_gauleg module

pynao.m_gauleg.gauleg_ab(n=96, a=- 1.0, b=1.0, eps=1e-15, cmx=15)

pynao.m_gauleg.gauss_legendre(n, charge=None, atom_id=None, **kwargs)

Generates the Gauss-Legendre knots and weights for using in 3D grids in pySCF For an example see m_ao_matelem.py

pynao.m_gauleg.leggauss_ab(n=96, a=- 1.0, b=1.0)

pynao.m_gaunt module

class pynao.m_gaunt.gaunt_c(jmax=7)

Bases: object

Computation of Gaunt coefficient (precompute then return)

get_gaunt(j1, m1, j2, m2)

pynao.m_get_atom2bas_s module

pynao.m_get_atom2bas_s.get_atom2bas_s(_bas)

For a given _bas list (see mole and mole_pure from pySCF) constructs a list of atom –> start shell The list is natoms+1 long, i.e. larger than number of atoms The list can be used to get the start,finish indices of pySCF’s multipletts^* This is useful to compose shls_slice arguments for the pySCF integral evaluators .intor(…)

^* pySCF multipletts can be “repeated” number of contractions times

pynao.m_get_sp_mu2s module

pynao.m_get_sp_mu2s.get_sp_mu2s(sp2nmult, sp_mu2j)

Generates list of start indices for atomic orbitals, based on the counting arrays

pynao.m_gpaw_hsx module

class pynao.m_gpaw_hsx.gpaw_hsx_c(sv, calc)

Bases: object

check_overlaps(pyscf_overlaps)

deallocate()

nelec

Comment about the overlap matrix in Gpaw: Marc: > From what I see the overlap matrix is not written in the GPAW output. Do > there is an option to write it or it is not implemented?

Ask answer (asklarsen@gmail.com): It is not implemented. Just write it manually from the calculation script.

If you use ScaLAPACK (or band/orbital parallelization) the array will be distributed, so each core will have only a chunk. One needs to call a function to collect it on master then. But probably you won’t need that.

The array will exist after you call calc.set_positions(atoms), in case you want to generate it without triggering a full calculation.

pynao.m_gpaw_wfsx module

class pynao.m_gpaw_wfsx.gpaw_wfsx_c(calc)

Bases: object

pynao.m_gw_chi0_noxv module

pynao.m_gw_chi0_noxv.gw_chi0_mv(self, dvin, comega)

pynao.m_gw_chi0_noxv.gw_chi0_mv_gpu(self, dvin, comega)

pynao.m_gw_xvx module

pynao.m_gw_xvx.gw_xvx_ac(self)

metod to calculate basis product using atom-centered product basis: V_{mu}^{ab}

pynao.m_gw_xvx.gw_xvx_ac_blas(self)

Metod to calculate basis product using atom-centered product basis: V_{mu}^{ab}

Use Blas to handle matrix-matrix multiplications

pynao.m_gw_xvx.gw_xvx_ac_simple(self)

Simple metod to calculate basis product using atom-centered product basis: V_{mu}^{ab}

This method is used as reference. direct multiplication with np and einsum

pynao.m_gw_xvx.gw_xvx_ac_sparse(self)

Metod to calculate basis product using atom-centered product basis: V_{mu}^{ab}

Use a sparse version of the atom-centered product, allow calculations of larger systems, however, the computational cost to get the sparse version of the atom-centered product is very high (up to 30% of the full simulation time of a GW0 run).

Warning

This method is experimental. For production run on large system, the dp_sparse method will be much more efficient in computational time and memory comsumption.

pynao.m_gw_xvx.gw_xvx_dp(self)

Metod to calculate basis product using dominant product basis V_{mu}^{ab}

pynao.m_gw_xvx.gw_xvx_dp_ndcoo(self)

pynao.m_gw_xvx.gw_xvx_dp_sparse(self)

Method to calculate basis product using dominant product basis V_{mu}^{ab}

Take advantages of the sparsity of the dominant product basis. Numba library must be installed

pynao.m_gw_xvx_dp_sparse module

pynao.m_gw_xvx_dp_sparse.csrmat_denmat_custom2(indptr, indices, data, B, X, N1, N2, N3)

Perform in one shot the following operations (avoiding the use of temporary arrays):

vxdp = v_pd1.dot(xmb.T) vxdp = vxdp.reshape(size,self.norbs, self.norbs) vxdp = np.swapaxes(vxdp,0,1) vxdp = vxdp.reshape(self.norbs,-1) xvx2 = xna.dot(xvx2)

The array vxdp is pretty large: (nfdp*norbs, norbs) and quite dense (0.1%) It is better to avoid its use for large systems, even in sparse format.

Inputs:

indptr: pointer index of v_dab in csr format indices: column indices of v_dab in csr format data: non zeros element of v_dab B: transpose of self.mo_coeff[0, spin, :, :, 0] X: self.mo_coeff[0, spin, self.nn[spin], :, 0] N1: nfdp N2: norbs N3: len(nn[spin])

Outputs:

D: xxv2 store in dense format

pynao.m_gw_xvx_dp_sparse.gw_xvx_sparse_dp(self)

Calculate the basis product using sparse dominant product basis

pynao.m_init_dens_libnao module

pynao.m_init_dens_libnao.init_dens_libnao()

Initilize the auxiliary for computing the density on libnao site

pynao.m_init_dm_libnao module

pynao.m_init_dm_libnao.init_dm_libnao(dm)

pynao.m_init_sab2dm_libnao module

pynao.m_init_sab2dm_libnao.init_dm_libnao(sab2dm)

pynao.m_ion_log module

class pynao.m_ion_log.ion_log_c(ao_log, sp)

Bases: object

holder of radial orbitals on logarithmic grid and bookkeeping info Args:

ion : ion structure (read from ion files from siesta)

Returns:

ao_log:

ion (ion structure from m_siesta_ion or m_siesta_ion_xml). nr (int): number of radial point rr pp mu2ff mu2ff_rl rcut (float): array containing the rcutoff of each specie mu2rcut (array, float) mu2s array containing pointers to the start indices for each radial multiplett norbs number of atomic orbitals interp_rr instance of log_interp_c to interpolate along real-space axis interp_pp instance of log_interp_c to interpolate along momentum-space axis

Examples:

>>> sv = system_vars()
>>> ao = ao_log_c(sv.sp2ion)
>>> print(ao.psi_log.shape)

pynao.m_kernel_utils module

pynao.m_kernel_utils.apply_kernel_nspin1_pack(spmv, nprod, kernel, dn, dtype=<class 'numpy.float64'>)

Perform the matrix vector multiplication for the pack kernel with chi0_mv in the iterative procedure

Input parameters:

spmv: blas function from scipy.linalg.blas

can be blas.sspmv for float or blas.dspmv for double

nprod: int

the dimension of the kernel matrix, i.e., (nprod, nprod)

kernel: 1D np.array, float or double

the kernel in pack format

dn: 1D np.array, complex

the resulting vector from chi0_mv

dtype: numpy data type

the data type of the kernel, np.float32 or np.float64

Output parameters:

vcre: 1D np.array, float

the real part of the matrix vector multiplication

vcim: 1D np.array, float

the imaginary part of the matrix vector multiplication

pynao.m_kernel_utils.apply_kernel_nspin1_sparse(nprod, kernel, kernel_diag, dn, dtype=<class 'numpy.float64'>)

Perform the matrix vector multiplication for the sparse kernel with chi0_mv in the iterative procedure

Input parameters:

nprod: int

the dimension of the kernel matrix, i.e., (nprod, nprod)

kernel: scipy csr_matrix, flat or double

the kernel in sparse format

dn: 1D np.array, complex

the resulting vector from chi0_mv

dtype: numpy data type

the data type of the kernel, np.float32 or np.float64

Output parameters:

vcre: 1D np.array, float or double

the real part of the matrix vector multiplication

vcim: 1D np.array, float or double

the imaginary part of the matrix vector multiplication

pynao.m_kernel_utils.apply_kernel_nspin2_pack(spmv, nprod, nspin, ss2kernel, dn, dtype=<class 'numpy.float64'>)

Perform the matrix vector multiplication for the pack kernel with chi0_mv in the iterative procedure for polarized spin calculation

Input parameters:

spmv: blas function from scipy.linalg.blas

can be blas.sspmv for float or blas.dspmv for double

nprod: int

the dimension of the kernel matrix, i.e., (nprod, nprod)

nspin: int

the number of spin

ss2kernel: list

the spin dependent kernel in pack format ss2kernel = [[kkk[0], kkk[1]],

[kkk[1], kkk[2]]

dn: 1D np.array, complex

the resulting vector from chi0_mv

dtype: numpy data type

the data type of the kernel, np.float32 or np.float64

Output parameters:

vcre: 2D np.array, float or double

the real part of the matrix vector multiplication per spin

vcim: 2D np.array, float or double

the imaginary part of the matrix vector multiplication per spin

pynao.m_kernel_utils.apply_kernel_nspin2_sparse(nprod, nspin, ss2kernel, ss2kernel_diag, dn, dtype=<class 'numpy.float64'>)

Perform the matrix vector multiplication for the sparse kernel with chi0_mv in the iterative procedure for polarized spin calculation

Input parameters:

nprod: int

the dimension of the kernel matrix, i.e., (nprod, nprod)

nspin: int

the number of spin

ss2kernel: list

the spin dependent kernel in sparse csr_matrix format ss2kernel = [[kkk[0], kkk[1]],

[kkk[1], kkk[2]]

dn: 1D np.array, complex

the resulting vector from chi0_mv

dtype: numpy data type

the data type of the kernel, np.float32 or np.float64

Output parameters:

vcre: 2D np.array, float or double

the real part of the matrix vector multiplication per spin

vcim: 2D np.array, float or double

the imaginary part of the matrix vector multiplication per spin

pynao.m_kernel_utils.get_VXC_kernel(self, Nspin, kernel, **kw)

Get the exchange-correlation part of the kernel and add it to Hartree kernel

Input parameters:

self: scf (scf.py) or tddft_iter (tddft_iter.py) class

the contructor that needs kernel initialization. Should contains all the necessary field to initialize the kernel

Nspin: int

Spin number

kernel: 1D np.array or sparse matrix

The Hartree kernel

kernel_format: string

the storage format of the kernel, can be * pack: store the kernel as a 1D array (lower part of the matrix) * sparse: store the lower part of the kernel in a sparse format.

kernel_threshold: float

Threshold used for sparse kernel in order to drop small matrix element and increase the sparsity of the matrix

pynao.m_kernel_utils.get_diagonal_triangle_matrix(arr, n)

pynao.m_kernel_utils.kernel_initialization(self, **kw)

Load the kernel from file or compute it

Inputs:

self: scf (scf.py) or tddft_iter (tddft_iter.py) class

the contructor that needs kernel initialization. Should contains all the necessary field to initialize the kernel

kernel_format: string

the storage format of the kernel, can be * pack: store the kernel as a 1D array (lower part of the matrix) * sparse: store the lower part of the kernel in a sparse format.

load_kernel: bool

Load the kernel from file (pack format only)

pynao.m_kernel_utils.load_kernel_pack(nprod, dtype, kernel_fname='', kernel_file_format='npy', kernel_path_hdf5='')

Loads the kernel from file in a pack format and initializes field

Inputs parameters:

nprod: int

kernel shape is (nprod, nprod)

dtype: numpy data type

data type of the saved kernel, should be np.float32 or np.float64

kernel_fname: string

file name of the kernel

kernel_file_format: string

format of the file: txt, npy or hdf5

kernel_path_hdf5: string

where to find the kernel in the hdf5 file

Outputs parameters:

kernel: 1D np.array, float

the loaded kernel in pack format, size nprod*(nprod + 1)//2

kernel_dim: tuple, int

the dimension (shape) of the unpacked kernel, (nprod, nprod)

pynao.m_kmat_den module

pynao.m_kmat_den.kmat_den(mf, dm=None, algo=None, **kw)

Computes the matrix elements of Fock exchange operator Args: mf : (System Variables), this must have arrays of coordinates and species, etc Returns: matrix elements

pynao.m_kmat_den.sm0_sum(pb, mf, hk, dm)

This algorithm is using two sparse representations of the product vertex V^ab_mu. The algorithm was not realized before and it seems to be superior to the algorithm sm0_prd (see above).

pynao.m_laplace_am module

pynao.m_laplace_am.laplace_am(self, sp1, R1, sp2, R2)

Computes brakets of Laplace operator for an atom pair. The atom pair is given by a pair of species indices and the coordinates of the atoms. Args:

self: class instance of ao_matelem_c sp1,sp2 : specie indices, and R1,R2 : respective coordinates

Result:

matrix of Laplace operator brakets

pynao.m_libnao module

Module to import the C and Fortran libraries

pynao.m_local_vertex module

class pynao.m_local_vertex.local_vertex_c(ao_log)

Bases: pynao.m_ao_matelem.ao_matelem_c

Constructor of the local product functions and the product vertex coefficients.

get_local_vertex(sp)

pynao.m_log_interp module

pynao.m_log_interp.comp_coeffs(self, r)

pynao.m_log_interp.comp_coeffs_(self, r, i2coeff)

Interpolation of a function given on the logarithmic mesh (see m_log_mesh how this is defined) 6-point interpolation on the exponential mesh (James Talman) Args:

r : radial coordinate for which we want the intepolated value

Result:

Array of weights to sum with the functions values to obtain the interpolated value coeff and the index k where summation starts sum(ff[k:k+6]*coeffs)

pynao.m_log_interp.log_interp(ff, r, rho_min_jt, dr_jt)

Interpolation of a function given on the logarithmic mesh (see m_log_mesh how this is defined) 6-point interpolation on the exponential mesh (James Talman) Args:

ff : function values to be interpolated r : radial coordinate for which we want intepolated value rho_min_jt : log(rr[0]), i.e. logarithm of minimal coordinate in the logarithmic mesh dr_jt : log(rr[1]/rr[0]) logarithmic step of the grid

Result:

Interpolated value

Example:

nr = 1024 rr,pp = log_mesh(nr, rmin, rmax, kmax) rho_min, dr = log(rr[0]), log(rr[1]/rr[0]) y = interp_log(ff, 0.2, rho, dr)

class pynao.m_log_interp.log_interp_c(gg)

Bases: object

Interpolation of radial orbitals given on a log grid (m_log_mesh)

coeffs(r)

Interpolation pointers and coefficients

coeffs_csr(rrs, rcut)

Compute a sparse array of interpolation coefficients (nr, rrs.shape) The subroutine returns also list of indices of non-zero rrs

coeffs_rcut(rr, rcut)

Compute an array of interpolation coefficients (6, rr.shape)

coeffs_vv(rr)

Compute an array of interpolation coefficients (6, rr.shape)

diff(za)

Return array with differential za : input array to be differentiated

interp_csr(ff, rrs, rcut=None)

Interpolation of vector data ff[…,:] and vector arguments rrs[:]. The function can accept also a scalar argument rrs

interp_rcut(ff, rr, rcut=None)

Interpolation of vector data ff[…,:] and vector arguments rr[:]

interp_v(ff, r)

Interpolation right away

interp_vv(ff, rr)

Interpolation of vector data ff[…,:] and vector arguments rr[:]

pynao.m_lorentzian module

pynao.m_lorentzian.limag_limag(x, w1, w2, e1, e2)

Product of two Lorentzians’ imaginary parts

pynao.m_lorentzian.llc_imag(x, w1, w2, e1, e2)

Product of two Lorentzians one of which is conjugated

pynao.m_lorentzian.llc_real(x, w1, w2, e1, e2)

Product of two Lorentzians one of which is conjugated

pynao.m_lorentzian.lorentzian(x, w, e)

An elementary Lorentzian

pynao.m_lorentzian.lreal_lreal(x, w1, w2, e1, e2)

Product of two Lorentzians’ real parts

pynao.m_lorentzian.overlap(ww, eps)

Overlap matrix between a set of Lorentzians using numerical integration

pynao.m_lorentzian.overlap_imag(ww, eps, wmax=inf)

Overlap matrix between a set of imaginary parts of Lorentzians using numerical integration

pynao.m_lorentzian.overlap_real(ww, eps, wmax=inf)

Overlap matrix between a set of real parts of Lorentzians using numerical integration

pynao.m_ls_contributing module

pynao.m_ls_contributing.ls_contributing(pb, sp12, ra12)

List of contributing centers prod_basis_c : instance of the prod_basis_c containing parameters .ac_rcut_ratio and .ac_npc_max and .sv providing instance of system_vars_c which provides the coordinates, unit cell vectors, species etc. and .prod_log prividing information on the cutoffs for each specie. sp12 : a couple of species ra12 : a couple of coordinates, correspondingly

pynao.m_ls_part_centers module

pynao.m_ls_part_centers.is_overlapping(rv1, rcut1, rv2, rcut2, rv3, rcut3, ac_rcut_ratio=1.0)

For a given atom pair (1,2) and a center 3 tell whether function at center overlaps with the product.

pynao.m_ls_part_centers.ls_part_centers(sv, ia1, ia2, ac_rcut_ratio=1.0)

For a given atom pair, defined with it’s atom indices, list of the atoms within a radii: list of participating centers. The subroutine goes over radial orbitals of atoms ia1 and ia2 computes the radii of products, the positions of centers of these products and then goes through atoms choosing these which overlap with the product.

pynao.m_next235 module

pynao.m_next235.next235(base)

pynao.m_numba_utils module

pynao.m_numba_utils.comp_coeffs_numba(gammin_jt, dg_jt, nr, r, i2coeff)

Interpolation of a function given on the logarithmic mesh (see m_log_mesh how this is defined)

6-point interpolation on the exponential mesh (James Talman)

Args:

r : radial coordinate for which we want the intepolated value

Result: Array of weights to sum with the functions values to obtain the interpolated value coeff and the index k where summation starts sum(ff[k:k+6]*coeffs)

pynao.m_numba_utils.csphar_numba(r, lmax)

Computes (all) complex spherical harmonics up to the angular momentum lmax

Args:

r : Cartesian coordinates defining correct theta and phi angles for spherical harmonic lmax : Integer, maximal angular momentum

Result:

1-d numpy array of complex128 elements with all spherical harmonics stored in order 0,0; 1,-1; 1,0; 1,+1 … lmax,lmax, althogether 0 : (lmax+1)**2 elements.

pynao.m_numba_utils.ls_contributing_numba(atom2coord, ia2dist)

pynao.m_openmx_import_scfout module

pynao.m_openmx_import_scfout.openmx_import_scfout(self, **kw)

Calls libopenmx to get the data and interpret it then

pynao.m_openmx_mat module

class pynao.m_openmx_mat.openmx_mat_c(natoms, Total_NumOrbs, FNAN, natn)

Bases: object

fromfile(f, out=None, dtype=<class 'float'>)

Read from an open file f

get_dims()

pynao.m_overlap_am module

pynao.m_overlap_am.overlap_am(self, sp1, R1, sp2, R2)

Computes overlap for an atom pair. The atom pair is given by a pair of species indices and the coordinates of the atoms. Args:

self: class instance of ao_matelem_c sp1,sp2 : specie indices, and R1,R2 : respective coordinates

Result:

matrix of orbital overlaps

The procedure uses the angular momentum algebra and spherical Bessel transform to compute the bilocal overlaps.

pynao.m_overlap_coo module

pynao.m_overlap_coo.overlap_coo(sv, ao_log=None, funct=<function overlap_ni>, ao_log2=None, **kw)

Computes the overlap matrix and returns it in coo format (simplest sparse format to construct)

Input parameters:

sv: System Variables, like nao class (pynao.nao)

this must have arrays of coordinates and species, etc

Output parameters:

overlap: coo_matrix (scipy.sparse.coo_matrix)

overlap (real-space overlap) for the whole system in coo sparse format

pynao.m_overlap_lil module

pynao.m_overlap_lil.overlap_lil(sv, ao_log=None, funct=<function overlap_ni>, **kvargs)

Computes the overlap matrix and returns it in List of Lists format (easy to index) Args: sv : (System Variables), this must have arrays of coordinates and species, etc Returns: overlap (real-space overlap) for the whole system

pynao.m_overlap_ni module

pynao.m_overlap_ni.overlap_ni(me, sp1, R1, sp2, R2, **kw)

Computes overlap for an atom pair. The atom pair is given by a pair of species indices and the coordinates of the atoms. Args:

sp1,sp2 : specie indices, and R1,R2 : respective coordinates in Bohr, atomic units

Result:

matrix of orbital overlaps

The procedure uses the numerical integration in coordinate space.

pynao.m_pack2den module

well described here http://www.netlib.org/lapack/lug/node123.html

pynao.m_pack2den.cp_block_pack_u(bmat, s1, f1, s2, f2, gmat, add=False)

pynao.m_pack2den.ij2pack_l(i, j, dim)

pynao.m_pack2den.ij2pack_u(i, j)

pynao.m_pack2den.pack2den_l(pack, dtype=<class 'numpy.float64'>)

Unpacks a packed format to dense format

pynao.m_pack2den.pack2den_u(pack, dtype=<class 'numpy.float64'>)

Unpacks a packed format to dense format

pynao.m_pack2den.size2dim(pack_size)

Computes dimension of the matrix by it’s packed size. Packed size is n*(n+1)//2

pynao.m_pb_ae module

pynao.m_pb_ae.pb_ae(self, sv, tol_loc=1e-05, tol_biloc=1e-06, ac_rcut_ratio=1.0)

It should work with GTOs as well.

pynao.m_phonons module

pynao.m_phonons.normal2cartesian(ma2xyz_nm, a2z)

Converts from normal coordinates (multiplied with sqrt of atomic masses) to Cartesian coordinates

pynao.m_phonons.read_vibra_vectors(fname=None)

Reads siesta.vectors file — output of VIBRA utility

pynao.m_phonons.read_xv(fname)

Reads siesta.XV file. All quantities are in Hartree atomic units

pynao.m_phonons.read_xyz(fname)

Reads xyz files

pynao.m_phonons.write_xyz(fname, s, ccc, line=None)

Writes xyz file

pynao.m_polariz_inter_ave module

pynao.m_polariz_inter_ave.polariz_freq_osc_strength(t2w, t2osc, comega)

Calculates polarizability for gives oscillator energy and strength. Useful with the TDDFT objects from pyscf. After executing tddft.kernel(), one can do

t2osc = tddft.oscillator_strength() t2e = tddft.e

and feed to this subroutine. The heavy lifting will be in PySCF.

pynao.m_polariz_inter_ave.polariz_inter_ave(mf, gto, tddft, comega)

pynao.m_polariz_inter_ave.polariz_inter_xx(mf, gto, tddft, comega)

pynao.m_polariz_inter_ave.polariz_nonin_ave(mf, gto, comega)

pynao.m_prod_basis_obsolete module

pynao.m_prod_biloc module

class pynao.m_prod_biloc.prod_biloc_c(atoms, vrtx, cc2a, cc2s, cc)

Bases: object

Holder of bilocal product vertices and conversion coefficients. Args:

atoms : atom pair (atom indices) vrtx : dominant product vertex coefficients: product,orb1,orb0 cc2a : contributing center -> atom index cc2s : contributing center -> start of the local product’s counting cc : conversion coefficients: product, atom-centered product

Returns:

structure with these fields

pynao.m_prod_biloc.set_prod_biloc(atoms, vrtx, cc2a, cc2s, cc)

pynao.m_prod_talman module

class pynao.m_prod_talman.prod_talman_c(lm=None, jmx=7, ngl=96, lbdmx=14)

Bases: pynao.log_mesh.log_mesh

prdred(phia, la, ra, phib, lb, rb, rcen)

Reduce two atomic orbitals given by their radial functions phia, phib, angular momentum quantum numbers la, lb and their centers ra,rb. The expansion is done around a center rcen.

prdred_further(ja, ma, jb, mb, rcen, jtb, clbdtb, lbdtb, rhotb)

Evaluate the Talman’s expansion at given Cartesian coordinates

prdred_further_scalar(ja, ma, jb, mb, rcen, jtb, clbdtb, lbdtb, rhotb)

Evaluate the Talman’s expansion at given Cartesian coordinates

prdred_libnao(phia, la, ra, phib, lb, rb, rcen)

By calling a subroutine

prdred_terms(la, lb)

Compute term-> Lambda,Lambda’,j correspondence

pynao.m_report module

pynao.m_report.get_energy_levels_FHL(self, egwev)

Return the values of the Fermi, HOMO and LUMO energies in eV

pynao.m_report.report_gw(self)

Prints the energy levels of mean-field and G0W0

pynao.m_report.report_gw_t(self, ret=None)

Lists spent time within main GW calculation

pynao.m_report.report_mfx(self, dm1=None)

This collects the h_core, Hartree (K) and Exchange(J) and Fock expectation value with given density matrix.

pynao.m_report.report_mocc(self)

Lists indeces of occupied and virtual orbitals in a limited range

pynao.m_report.sigma_xc(self)

This calculates the Exchange expectation value and correlation part of self energy, when: self.get_k() = Exchange operator/energy mat1 is product of this operator and molecular coefficients and it will be diagonalized in expval by einsum Sigma_c = E_GW - E_HF

pynao.m_restart module

pynao.m_restart.check_res(filename)

pynao.m_restart.find(pattern, path)

pynao.m_restart.read_rst_h5py(value, filename=None, arr=None)

pynao.m_restart.read_rst_yaml(filename=None)

pynao.m_restart.write_rst_h5py(data, value, filename=None)

pynao.m_restart.write_rst_yaml(data, filename=None)

pynao.m_rf0_den module

pynao.m_rf0_den.add_outer(a, b)

pynao.m_rf0_den.calc_XVX(X, VX, dtype=<class 'numpy.float64'>)

pynao.m_rf0_den.calc_XVX_numba(X, VX, dtype=<class 'numpy.float64'>)

pynao.m_rf0_den.calc_part_rf0(xvx, zxvx_re, zxvx_im)

pynao.m_rf0_den.calc_part_rf0_numba(xvx, zxvx_re, zxvx_im)

pynao.m_rf0_den.rf0_cmplx_ref(self, ww)

Full matrix response in the basis of atom-centered product functions

pynao.m_rf0_den.rf0_cmplx_ref_blk(self, ww)

Full matrix response in the basis of atom-centered product functions

pynao.m_rf0_den.rf0_cmplx_vertex_ac(self, ww)

Full matrix response in the basis of atom-centered product functions

pynao.m_rf0_den.rf0_cmplx_vertex_dp(self, ww)

Full matrix response in the basis of atom-centered product functions

pynao.m_rf0_den.rf0_den(self, ww)

Full matrix response in the basis of atom-centered product functions for parallel spins.

Blas version to speed up matrix matrix multiplication spped up of 7.237 compared to einsum version for C20 system

pynao.m_rf0_den.rf0_den_einsum(self, ww)

Full matrix response in the basis of atom-centered product functions for parallel spins

einsum version, slow

pynao.m_rf0_den.rf0_den_numba(rf0, comega, X, ksn2f, ksn2e, pab2v_den, nprod, norbs, bsize, nspin, nfermi, vstart, dtype=<class 'numpy.complex128'>)

Full matrix response in the basis of atom-centered product functions for parallel spins.

Blas version to speed up matrix matrix multiplication spped up of 7.237 compared to einsum version for C20 system

pynao.m_rf0_den.run_blasDGEMM(alpha, A, B, **kwargs)

pynao.m_rf0_den.run_blasZGEMM(alpha, A, B, **kwargs)

pynao.m_rf0_den.si_correlation(rf0, si0, ww, kernel_sq, nprod)

This computes the correlation part of the screened interaction W_c by solving <self.nprod> linear equations (1-K chi0) W = K chi0 K or v_{ind}sim W_{c} = (1-vchi_{0})^{-1}vchi_{0}v scr_inter[w,p,q], where w in ww, p and q in 0..self.nprod

pynao.m_rf0_den.si_correlation_numba(si0, ww, X, kernel_sq, ksn2f, ksn2e, pab2v_den, nprod, norbs, bsize, nspin, nfermi, vstart)

This computes the correlation part of the screened interaction W_c by solving <self.nprod> linear equations (1-K chi0) W = K chi0 K or v_{ind}sim W_{c} = (1-vchi_{0})^{-1}vchi_{0}v scr_inter[w,p,q], where w in ww, p and q in 0..self.nprod

pynao.m_rf_den module

pynao.m_rf_den.rf_den(self, ww)

Full matrix interacting response from NAO GW class

pynao.m_rf_den.rf_den_via_rf0(self, rf0, v)

Whole matrix of the interacting response via non-interacting response and interaction

pynao.m_rf_den_pyscf module

pynao.m_rf_den_pyscf.rf_den_pyscf(self, ww)

Full matrix interacting response from tdscf class

pynao.m_rsphar module

pynao.m_rsphar.rsphar(r, lmax, res)

Computes (all) real spherical harmonics up to the angular momentum lmax Args:

r : Cartesian coordinates defining correct theta and phi angles for spherical harmonic lmax : Integer, maximal angular momentum

Result:

1-d numpy array of float64 elements with all spherical harmonics stored in order 0,0; 1,-1; 1,0; 1,+1 … lmax,lmax, althogether 0 : (lmax+1)**2 elements.

pynao.m_rsphar_libnao module

pynao.m_rsphar_libnao.rsphar(r, lmax, res)

Computes (all) real spherical harmonics up to the angular momentum lmax Args:

r : Cartesian coordinates defining correct theta and phi angles for spherical harmonic lmax : Integer, maximal angular momentum

Result:

1-d numpy array of float64 elements with all spherical harmonics stored in order 0,0; 1,-1; 1,0; 1,+1 … lmax,lmax, althogether 0 : (lmax+1)**2 elements.

pynao.m_rsphar_libnao.rsphar_exp_vec(rvs, lmax)

Computes (all) real spherical harmonics up to the angular momentum lmax Args:

rvs : Cartesian coordinates defining correct theta and phi angles for spherical harmonic lmax : Integer, maximal angular momentum

Result:

1-d numpy array of float64 elements with all spherical harmonics stored in order 0,0; 1,-1; 1,0; 1,+1 … lmax,lmax, althogether 0 : (lmax+1)**2 elements.

pynao.m_rsphar_libnao.rsphar_vec(rvs, lmax)

Computes (all) real spherical harmonics up to the angular momentum lmax Args:

rvs : Cartesian coordinates defining correct theta and phi angles for spherical harmonic lmax : Integer, maximal angular momentum

Result:

1-d numpy array of float64 elements with all spherical harmonics stored in order 0,0; 1,-1; 1,0; 1,+1 … lmax,lmax, althogether 0 : (lmax+1)**2 elements.

pynao.m_rsphar_vec module

pynao.m_rsphar_vec.rsphar_vec(rvecs, lmax)

Computes (all) real spherical harmonics up to the angular momentum lmax Args:

rvecs : A list of Cartesian coordinates defining the theta and phi angles for spherical harmonic lmax : Integer, maximal angular momentum

Result:

2-d numpy array of float64 elements with all spherical harmonics stored in order 0,0; 1,-1; 1,0; 1,+1 … lmax,lmax, althogether 0 : (lmax+1)**2 elements.

pynao.m_sbt module

class pynao.m_sbt.sbt_c(rr, kk, lmax=12, with_sqrt_pi_2=True, fft_flags=None)

Bases: object

Spherical Bessel Transform by James Talman. Functions are given on logarithmic mesh See m_log_mesh Args:

nr : integer, number of points on radial mesh rr : array of points in coordinate space kk : array of points in momentum space lmax : integer, maximal angular momentum necessary with_sqrt_pi_2 : if one, then transforms will be multiplied by sqrt(pi/2) fft_flags : ??

Returns:

a class preinitialized to perform the spherical Bessel Transform

Examples:

label = ‘siesta’ sv = system_vars_c(label) sbt = sbt_c(sv.ao_log.rr, sv.ao_log.pp) print(sbt.exe(sv.ao_log.psi_log[0,0,:], 0))

sbt(ff, am, direction=1, npow=0)

Args:

ff : numpy array containing radial orbital (values of radial orbital on logarithmic grid) to be transformed. The data must be on self.rr grid or self.kk grid provided during initialization. am : angular momentum of the radial orbital ff[:] direction : 1 – real-space –> momentum space transform; -1 – momentum space –> real-space transform. npow : additional power for the shape of the orbital

f(xyz) = rr[i]**npow * ff[i] * Y_lm( xyz )

Result:

gg : numpy array containing the result of the Spherical Bessel Transform gg(k) = int_0^infty ff(r) j_{am}(k*r) r**2 dr ( direction == 1 ) gg(r) = int_0^infty ff(k) j_{am}(k*r) k**2 dk ( direction == -1 )

pynao.m_sf2f_rf module

pynao.m_sf2f_rf.sf2f_rf(ff, eps, w2f, wab2sf)

pynao.m_siesta2blanko_csr module

pynao.m_siesta2blanko_denvec module

pynao.m_siesta_eig module

pynao.m_siesta_hsx module

class pynao.m_siesta_hsx.siesta_hsx_c(**kw)

Bases: object

Interface to the HSX data from Siesta calculations

deallocate()

Deallocate hsx variables

pynao.m_siesta_hsx.siesta_hsx_read(fname, force_gamma=None)

Read the siesta.HSX file

pynao.m_siesta_hsx_bloch_mat module

pynao.m_siesta_hsx_bloch_mat.siesta_hsx_bloch_mat(csr, hsx, kvec=array([0.0, 0.0, 0.0]))

pynao.m_siesta_hsx_bloch_mat.siesta_hsx_bloch_mat_py(csr, hsx, kvec=[0.0, 0.0, 0.0], prec=64)

pynao.m_siesta_ion module

class pynao.m_siesta_ion.siesta_ion_c(splabel)

Bases: object

pynao.m_siesta_ion_add_sp2 module

pynao.m_siesta_ion_xml module

pynao.m_siesta_ion_xml.extract_field_elements(doc, field=None)

extract the different pao element of the xml file Input Parameters: —————–

field: field name of the node doc (minidom.parse)

pao (dict): the following keys are added to the ion dict:

npts delta cutoff data orbital

pynao.m_siesta_ion_xml.extract_orbital(orb_xml)

extract the orbital

pynao.m_siesta_ion_xml.extract_projector(pro_xml)

extract the projector

pynao.m_siesta_ion_xml.getNodeText(node)

pynao.m_siesta_ion_xml.get_data_elements(name, dtype)

return the right type of the element value

pynao.m_siesta_ion_xml.siesta_ion_xml(fname)

Read the ion.xml file of a specie Input parameters: —————– fname (str): name of the ion file Output Parameters: —————— ion (dict): The ion dictionnary contains all the data

from the ion file. Each field of the xml file give one key. The different keys are:

‘lmax_basis’: int ‘self_energy’: float ‘z’: int ‘symbol’: str ‘label’: str ‘mass’: flaot ‘lmax_projs’: int ‘basis_specs’: str ‘norbs_nl’: int ‘valence’: float ‘nprojs_nl: int

The following keys give the pao field, ‘npts’: list of int ‘delta’:list of float ‘cutoff’: list of float ‘data’:list of np.arrayof shape (npts[i], 2) ‘orbital’: list of dictionary ‘projector’: list of dictionary

pynao.m_siesta_ion_xml.str2float(string)

pynao.m_siesta_ion_xml.str2int(string)

pynao.m_siesta_units module

Units conversion for Siesta outputs

Used the definition from Scipy https://docs.scipy.org/doc/scipy/reference/constants.html

pynao.m_siesta_units.get_conversion_factors(definition_year=2014)

Possibility to used different definition of the physicale constants as defined by the CODATA

By default, PyNAO used the definitions from the year 2014 (even if they are not the latest values). This is because the tests for PyNAO were set up with this definition. Changing the year may make some tests fails.

References

Theoretical and experimental publications relevant to the fundamental constants and closely related precision measurements published since the mid 1980s, but also including many older papers of particular interest, some of which date back to the 1800s. To search bibliography visit

https://physics.nist.gov/cuu/Constants/

pynao.m_siesta_utils module

pynao.m_siesta_utils.get_pseudo(sp, suffix='')

return the path to the pseudopotential of a particular specie

pynao.m_siesta_utils.get_siesta_command(label, directory='')

pynao.m_siesta_utils.install_siesta()

pynao.m_siesta_utils.runbash(cmd)

Run a bash command with subprocess, fails if ret != 0

pynao.m_siesta_wfsx module

pynao.m_siesta_wfsx.siesta_wfsx_book_read_py(fname, nreim)

Creates buffer for integer data from .WFSX files

class pynao.m_siesta_wfsx.siesta_wfsx_c(**kw)

Bases: object

pynao.m_siesta_wfsx.siesta_wfsx_dread(w, nreim)

pynao.m_siesta_wfsx.siesta_wfsx_sread(w, sdata, nreim)

pynao.m_siesta_xml module

pynao.m_siesta_xml.extract_energy_parameter(parent, field)

Extract an energy parameter from the siesta.xml file and convert it to Hartree

pynao.m_siesta_xml.siesta_xml(fname='siesta.xml')

Read the siesta.xml file

pynao.m_siesta_xml_print module

pynao.m_siesta_xml_print.siesta_xml_print(label='siesta')

pynao.m_simulation module

class pynao.m_simulation.simulation_c(sv, **kvargs)

Bases: object

do_overlap_check_of_pb(**kvargs)

pynao.m_sparsetools module

pynao.m_sparsetools.count_nnz_spmat_denmat(csr, ncolB)

pynao.m_sparsetools.csc_matvec(csc, x)

Matrix vector multiplication using csc format

pynao.m_sparsetools.csc_matvecs(csc, B, transB=False, order='C')

Matrix matrix multiplication using csc format

pynao.m_sparsetools.csr_matvec(csr, x, y=None)

pynao.m_sparsetools.spmat_denmat(csr, B)

sparse matrix dense matrix multiplication directly giving a sparse matrix in CSR format.

pynao.m_spline_diff2 module

pynao.m_spline_diff2.spline_diff2(h, yin, yp1, ypn)

subroutine spline(delt,y,n,yp1,ypn,y2) !! Cubic Spline Interpolation. !! Adapted from Numerical Recipes routines for a uniform grid !! D. Sanchez-Portal, Oct. 1996. !! Alberto Garcia, June 2000 !! Peter Koval, Dec 2009

implicit none !! external integer, intent(in) :: n real(8), intent(in) :: delt, yp1, ypn, y(:) real(8), intent(out) :: y2(:)

!! internal integer i, k real(8) sig, p, qn, un

real(8), allocatable :: u(:) allocate(u(n));

if (yp1.eq. huge(1D0)) then

y2(1)=0 u(1)=0

else

y2(1)=-0.5D0 u(1)=(3.0D0/delt)*((y(2)-y(1))/delt-yp1)

endif

do i=2,n-1

sig=0.5D0 p=sig*y2(i-1)+2 y2(i)=(sig-1)/p u(i)=(3*( y(i+1)+y(i-1)-2*y(i) )/(delt*delt)-sig*u(i-1))/p

enddo

if (ypn.eq.huge(1D0)) then

qn=0; un=0

else

qn=0.5D0; un=(3/delt)*(ypn-(y(n)-y(n-1))/delt)

endif

y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1) do k=n-1,1,-1

y2(k)=y2(k)*y2(k+1)+u(k)

enddo

end subroutine !spline

pynao.m_spline_interp module

pynao.m_spline_interp.spline_interp(h, yy, yy_diff2, x)

!subroutine splint(delt,ya,y2a,n,x,y,dydx) !! Cubic Spline Interpolation. !! Adapted from Numerical Recipes for a uniform grid.

implicit none !! external integer, intent(in) :: n real(8), intent(in) :: delt, ya(n), y2a(n), x real(8), intent(out) :: y

! real(dp), intent(out) :: y, dydx

!! internal integer :: nlo, nhi real(8) :: a, b

nlo=max(int(x/delt)+1,1) if(nlo>n-1) then; y=0; return; endif !if(nlo>n-1) then; y=0; dydx=0; return; endif nhi=min(nlo+1,n) a=nhi-x/delt-1 b=1.0D0-a y=a*ya(nlo)+b*ya(nhi)+((a**3-a)*y2a(nlo)+(b**3-b)*y2a(nhi))*(delt**2)/6D0

! dydx=(ya(nhi)-ya(nlo))/delt + (-((3*(a**2)-1._dp)*y2a(nlo))+ (3*(b**2)-1._dp)*y2a(nhi))*delt/6._dp end subroutine ! splint

pynao.m_sv_chain_data module

pynao.m_sv_chain_data.sv_chain_data(sv)

This subroutine creates a buffer of information to communicate the system variables to libnao.

pynao.m_sv_diag module

pynao.m_sv_diag.sv_diag(sv, kvec=[0.0, 0.0, 0.0], spin=0, prec=64)

pynao.m_sv_get_denmat module

pynao.m_sv_get_denmat.sv_get_denmat(sv, mattype='hamiltonian', prec=64, kvec=[0.0, 0.0, 0.0], spin=0)

pynao.m_system_vars_deprecated module

pynao.m_system_vars_dos module

pynao.m_system_vars_dos.eigenvalues2dos(ksn2e, zomegas, nkpoints=1)

Compute the Density of States using the eigenvalues

pynao.m_system_vars_dos.system_vars_ados(sv, zomegas, ls_atom_groups, nkpoints=1)

Compute a Partial Density of States (resolved in atomic indices) using the eigenvalues and eigenvectors in wfsx

pynao.m_system_vars_dos.system_vars_dos(sv, zomegas, nkpoints=1)

Compute the Density of States using the eigenvalues

pynao.m_system_vars_dos.system_vars_pdos(sv, zomegas, nkpoints=1)

Compute the Partial Density of States (resolved in angular momentum of the orbitals) using the eigenvalues and eigenvectors in wfsx

pynao.m_system_vars_gpaw module

pynao.m_system_vars_gpaw.system_vars_gpaw(self, **kw)

Initialise system variables with gpaw inputs: Input parameters: —————–

calc: GPAW object kvargs: optional arguments, need a list of precise arguments somewhere

pynao.m_tddft_iter_gpu module

class pynao.m_tddft_iter_gpu.tddft_iter_gpu_c(GPU, X4, ksn2f, ksn2e, norbs, nfermi, nprod, vstart)

Bases: object

calc_nb2v_from_nm2v_imag()

calc_nb2v_from_nm2v_real()

calc_nb2v_from_sab(reim)

calc_nm2v_imag()

calc_nm2v_real()

calc_sab(reim)

clean_gpu()

countGPUs()

Return the number of devices available for the calculations

cpy_sab_to_device(sab, Async=- 1)

Async can take the following values:

• 0 default stream

• 1 real stream

• 2 imag stream

• -1 or any other value: Not using Async, just blocking memcpy

cpy_sab_to_host(sab, Async=- 1)

Async can take the following values:

• 0 default stream

• 1 real stream

• 2 imag stream

• -1 or any other value: Not using Async, just blocking memcpy

div_eigenenergy_gpu(comega)

getDevice()

setDevice(gpu_id)

pynao.m_thrj module

pynao.m_thrj.comp_number_of3j(lmax)

Computes number of 3j coefficients irreducible by symmetry (as implemented in thrj(l1,l2,l3,m1,m2,m3) )

pynao.m_thrj.thrj(l1i, l2i, l3i, m1i, m2i, m3i)

Wigner3j symbol. Written by James Talman.

pynao.m_thrj.thrj_nobuf(l1, l2, l3, m1, m2, m3)

Wigner3j symbol without buffer. Written by James Talman.

pynao.m_tools module

modules containing tools and utility functions

pynao.m_tools.is_power2(n)

Check if n is a power of 2

pynao.m_tools.read_xyz(fname)

Reads xyz files

pynao.m_tools.rtp2xyz(r, t, p)

pynao.m_tools.transformData2newCoordinate(oldCoordinates, newCoordinates, data, transform=<function rtp2xyz>)

transform a 3D array from a coodinate system to another. For example, transforming from cartesian to spherical coordinates:

from __future__ import division import numpy as np from pynao.m_tools import transformData2newCoordinate

dims = (10, 5, 6) x = np.linspace(-5, 5, dims[0]) y = np.linspace(-2, 2, dims[1]) z = np.linspace(-3, 3, dims[2])

dn = np.random.randn(dims[0], dims[1], dims[2])

r = np.arange(0.0, 2.0, 0.1) phi = np.arange(0.0, 2*np.pi, 0.01) theta = np.arange(0.0, np.pi, 0.01)

dn_new = transformData2newCoordinate((x, y, z), (r, phi, theta), dn)

pynao.m_tools.write_xyz(fname, s, ccc)

Writes xyz files

pynao.m_tools.xyz2rtp(x, y, z)

pynao.m_valence module

pynao.m_valence.get_str_fin(self, frozen_core, **kw)

pynao.m_valence.n_core(self)

pynao.m_valence.n_core_vale(self)

pynao.m_valence.n_valence(self)

pynao.m_vertex_loop module

class pynao.m_vertex_loop.vertex_loop_c(pb)

Bases: object

Holder of dominant product vertices in the form which must be easy to deploy in iterative tddft or in Fock operator Args:

instance of prod_basis_c

Returns:

%.vdata : all product vertex blocks packed into an one dimensional array %.i2inf : integer array iteration –> info (start in vdata, atom1,atom2,product center,start functions,finish functions)

Examples:

pynao.m_vhartree_coo module

pynao.m_vhartree_coo.dm2j_fullmat(n, v_dab, da2cc, hk, dm)

pynao.m_vhartree_coo.vhartree_coo(mf, dm=None, **kw)

Computes the matrix elements of Hartree potential

Input parameters:

mf: mf class

this must have arrays of coordinates and species, etc

Output parameters:

vh_coo: coo_matrix (scipy.sparse.coo_matrix)

matrix elements stored in coordinate sparse format

pynao.m_vhartree_pbc module

pynao.m_vhartree_pbc.apply_inv_G2(vh, gg0, gg1, gg2)

Only critics to this solution is the memory concerns

pynao.m_vhartree_pbc.apply_inv_G2_ref(vh, gg0, gg1, gg2)

pynao.m_vhartree_pbc.vhartree_pbc(self, dens, **kw)

Compute Hartree potential for the density given in an equidistant grid

pynao.m_vnucele_coo module

pynao.m_vnucele_coo.vnucele_coo(sv, algo=None, **kvargs)

Computes the nucleus-electron attraction matrix elements Args:

sv : (System Variables), this must have arrays of coordinates and species, etc

Returns:

matrix elements

These are tricky to define. In case of all-electron calculations it is well known, but in case of pseudo-potential calculations we need some a specification of pseudo-potential and a specification of (Kleinman-Bulander) projectors to compute explicitly. Practically, we will subtract the computed matrix elements from the total Hamiltonian to find out the nuclear-electron interaction in case of SIESTA import. This means that Vne is defined by

Vne = H_KS - T - V_H - V_xc

pynao.m_vnucele_coo_subtract module

pynao.m_vnucele_coo_subtract.vnucele_coo_subtract(sv, **kvargs)

Computes the matrix elements defined by

Vne = H_KS - T - V_H - V_xc

which serve as nuclear-electron attraction matrix elements for pseudo-potential DFT calculations Args:

sv : (System Variables), this must have arrays of coordinates and species, etc

Returns:

matrix elements

pynao.m_vxc_lil module

pynao.m_vxc_pack module

pynao.m_vxc_pack.vxc_pack(self, **kw)

Computes the exchange-correlation matrix elements packed version (upper triangular) Args:

sv : (System Variables), this must have arrays of coordinates and species, etc

Returns:

vxc,exc

pynao.m_x_zip module

pynao.m_x_zip.detect_maxima(ww, dos)

Detects maxima of a function given on a grid and lists the arguments at maxima

pynao.m_x_zip.ee2dos(ee, eps)

convert a set of energies to a grid and a density of states

pynao.m_x_zip.ee_xx_oo2dos(m2e, ma2x, ab2o, eps)

convert a set of energies, eigenvectors and overlap to a grid and a density of states

pynao.m_x_zip.x_zip(n2e, na2x, eps, emax)

Construct a ‘presummed’ set of eigenvectors, the effect must be a smaller number of eigenvectors

pynao.m_xc_scalar_ni module

pynao.m_xc_scalar_ni.xc_scalar_ni(me, sp1, R1, sp2, R2, xc_code, deriv, **kw)

pynao.m_xjl module

pynao.m_xjl.xjl(xx, lc)

Spherical bessel functions Computes a table of j_l(x) for fixed xx, Eq. (39) Args:

xx : float lc : integer angular momentum

Result:

xj[0:lc] : float

pynao.m_xjl_numba module

pynao.m_xjl_numba.c2r_numba(j, tr_c2r, conj_c2r, j1, j2, jm, cmat, rmat, mat)

pynao.m_xjl_numba.calc_oo2co(bessel_pp, dg_jt, ao1_sp2info_sp1, ao1_sp2info_sp2, ao1_psi_log_mom_sp1, ao1_psi_log_mom_sp2, njm, gaunt_iptr, gaunt_data, ylm, j, jmx, tr_c2r, conj_c2r, l2S, cS, rS, cmat, oo2co)

pynao.m_xjl_numba.get_bessel_xjl_numba(kk, dist, j, nr)

Calculate spherical bessel functions j_l(k*r) for a given r and on a grid of momentas k given by the array kk Args: kk : 1D array (float): k grid r : (float) radial coordinate l : (integer) angular momentum nr: (integer) k grid dimension Result: xj[1:2*j+1, 1:nr] : 2D array (float)

pynao.mesh_affine_equ module

class pynao.mesh_affine_equ.mesh_affine_equ(**kw)

Bases: object

get_3dgrid()

Generate 3d Grid a la PySCF with .coords and .weights fields

get_all_coords()

get_mesh_center()

write(fname, **kw)

pynao.mf module

class pynao.mf.mf(**kw)

Bases: pynao.nao.nao

dens_elec(coords, dm)

diag_check(atol=1e-05, rtol=0.0001)

dos(comegas, **kw)

Ordinary Density of States (from the current mean-field eigenvalues)

exc(dm=None, xc_code=None, **kw)

gdos(comegas, **kw)

Some molecular orbital population analysis

get_hamiltonian()

get_init_guess(mol=None, key=None)

Compute an initial guess for the density matrix.

get_occupations(telec=None, ksn2e=None, fermi_energy=None)

Compute occupations of electron levels according to Fermi-Dirac distribution

get_vertex_pov()

Computes the occupied-virtual-product basis vertex

init_libnao(wfsx=None)

Initialization of data on libnao site

init_mo_coeff_fireball(**kw)

Constructor a mean-field class from the preceeding FIREBALL calculation

init_mo_from_pyscf(**kw)

Initializing from a previous pySCF mean-field calc.

lsoa_dos(comegas, **kw)

Partial Density of States (contributions from a given list of atoms)

nonin_osc_strength()

pdos(comegas, **kw)

Partial Density of States (resolved in angular momentum of atomic orbitals)

plot_contour(w=0.0)

Plot contour with poles of Green’s function in the self-energy SelfEnergy(w) = G(w+w’)W(w’) with respect to w’ = Re(w’)+Im(w’) Poles of G(w+w’) are located: w+w’-(E_n-Fermi)+i*eps sign(E_n-Fermi)==0 ==> w’= (E_n-Fermi) - w -i eps sign(E_n-Fermi)

polariz_nonin_ave_matelem(comega)

read_wfsx(fname, **kw)

An occasional reading of the SIESTA’s .WFSX file

spin_square(mo_coeff=None, mo_occ=None)

vhartree_pbc(dens, **kw)

Compute Hartree potential for the density given in an equidistant grid

vhartree_pbc_coo(density_factors=[1, 0], **kw)

Compute matrix elements of Hartree potential for the density given in an equidistant grid

vxc_lil(**kw)

Compute matrix elements of exchange-correlation potential

pynao.nao module

exception pynao.nao.SystemNotCentered

Bases: Exception

pynao.nao.check_occupation(self)

Check that the calculated number of electron obtained from siesta wfsx correspond to the expected number of electrons

pynao.nao.check_system_centered(coordinates, raise_error=True, eps=1.0, verbose=2)

Check if the system is centered. Non centered system can lead at wrong induced density, particularly for large system. See issue for further details.

Input parameters:

coordinates ([Natoms, 3] float array): Coordinate of the atoms in a.u.

raise_error (bool): Set to False to skip the check

eps (float): Threshold for the check, by default 1.0 a.u.

pynao.nao.get_nelec_occ(self, correction=1.0)

calculate the number of electrons from WFSX file

pynao.nao.get_orb2j(sv)

pynao.nao.get_orb2m(sv)

class pynao.nao.nao(**kw)

Bases: object

ao_loc_nr()

atom_charge(ia)

atom_charges()

atom_coord(ia)

atom_coords()

atom_nelec_core(ia)

atom_symbol(ia)

atom_symbols()

build_3dgrid_ae(level=3)

Build a global grid and weights for a molecular integration (integration in 3-dimensional coordinate space)

build_3dgrid_pp(level=3)

Build a global grid and weights for a molecular integration (integration in 3-dimensional coordinate space)

comp_aos_csr(coords, tol=1e-08, ram=160000000.0)

Compute the atomic orbitals for a given set of (Cartesian) coordinates. The sparse format CSR is used for output and the computation is organized block-wise. Thence, larger molecules can be tackled right away

Input parameters:

coords: 2D np.array, float

set of Cartesian coordinates

tol: float

tolerance for dropping the values

ram: float

size of the allowed block (in bytes)

Output parameters:

co2v: sparse csr matrix (scipy.sparse.csr_matrix), float

CSR matrix of shape (coordinate, atomic orbital)

comp_aos_den(coords)

Compute the atomic orbitals for a given set of (Cartesian) coordinates.

comp_aos_py(coords)

Compute the atomic orbitals for a given set of (Cartesian) coordinates.

comp_vnuc_coulomb(coords)

dipole_coo(**kw)

Compute dipole matrix elements for the given system

energy_nuc(charges=None, coords=None)

Potential energy of electrostatic repulsion of point nuclei

get_init_guess(key=None)

Compute an initial guess for the density matrix. ????

get_orb2j()

get_orb2m()

get_symbols()

get_vkb()

Compose the vector of Kleinman-Bylander energies v^p = v^KB_ln, where p is a global projector index

init_fireball(**kw)

init_gpaw(**kw)

Use the data from a GPAW LCAO calculations as input to initialize system variables.

init_gto(**kw)

Interpret previous PySCF calculation

init_guess_by_minao(key=None)

Compute an initial guess for the density matrix. ????

init_label(**kw)

init_libnao_orbs()

Initialization of data on libnao site

init_mo_coeff_label(**kw)

Constructor a mean-field class from the preceeding SIESTA calculation

init_openmx(**kw)

init_wfsx_fname(**kw)

Initialise system var starting with a given WFSX file

init_xyz_list(**kw)

This is simple constructor which only initializes geometry info

intor_symmetric(intor, comp=None)

laplace_coo()

Compute matrix of Laplace brakets for the whole molecule

make_rdm1(mo_coeff=None, mo_occ=None)

matelem_int3d_coo(g, v)

Compute matrix elements of a potential v given on the 3d grid g using blocks along the grid

matelem_int3d_coo_ref(g, v)

Compute matrix elements of a potential v given on the 3d grid g

nao_nr()

property nelectron

overlap_check(**kw)

Works for pyscf and after init_siesta_xml()

overlap_coo(**kw)

Compute overlap matrix for the molecule

overlap_lil(**kw)

Compute overlap matrix in list of lists format

ucell_mom()

vna(coords, **kw)

Compute the neutral-atom potential V_NA(coords) for a set of Cartesian coordinates coords. The subroutine could be also used for computing the non-linear core corrections or some other atom-centered fields.

vna_coo(**kw)

Compute matrix elements of a potential which is given as superposition of central fields from each nuclei

vnl_coo()

Non-local part of the Hamiltonian due to Kleinman-Bylander projectors

vnucele_coo(**kw)

vnucele_coo_coulomb(**kw)

vnucele_coo_pseudo(**kw)

Compute matrix elements of attraction by forces from pseudo atom

pynao.nao.overlap_check(sv, tol=1e-05, **kvargs)

pynao.ndcoo module

class pynao.ndcoo.ndcoo(inp, **kw)

Bases: object

tocoo_a_pb(descr)

converts shape of sparse matrix () into (a, p*b)

tocoo_b_pa(descr)

converts shape of sparse matrix () into (b, p*a)

tocoo_p_ab(descr)

converts shape of sparse matrix () into (p, a*b)

tocoo_pa_b(descr)

converts shape of sparse matrix () into (p*a, b)

pynao.prod_basis module

class pynao.prod_basis.prod_basis(nao=None, **kw)

Bases: object

prod_basis object.

Holder of local and bilocal product functions and vertices.

Input parameters:

nao: nao class (pynao.nao)

holder of the geometry, and orbitals description

tol_loc: float

tolerance for local basis

tol_biloc: float

tolerance for bilocal basis

tol_elim: float

tolerance for ??

ac_rcut_ratio: float

ac rcut ratio??

ac_npc_max: int

maximal number of participating centers

pb_algorithm: string

algorithm to use for product basis

jcutoff: int

number of angular orbital?? included

metric_type: int

type of metric ??

optimize_centers: int

optimze the center or not

ngl: int

???

Output parameters:

For each specie returns a set of radial functions defining a product basis These functions are sufficient to represent the products of original atomic orbitals via a product vertex coefficients and conversion coefficients.

prod_log: class prod_log (pynao.prod_log)

Holder of (local) product functions and vertices

hkernel_csr: csr sparse matrix (scipy.sparse.csr)

hartree kernel: local part of Coulomb interaction

c2s: 1D np.array, int, size natoms+1

global product Center (atom) -> start in case of atom-centered basis

bp2info: list

some information including indices of atoms, list of contributing centres, conversion coefficients

dpc2s: list

product Center -> list of the size of the basis set in this center, of center’s types,of product species

dpc2t: list

product Center -> list of the size of the basis set in this center, of center’s types,of product species

dpc2sp: list

product Center -> list of the size of the basis set in this center, of center’s types,of product species

chain_data()

This subroutine creates a buffer of information to communicate the system variables and the local product vertex to libnao. Later, one will be able to generate the bilocal vertex and conversion coefficient for a given pair of atom species and their coordinates.

comp_apair_pp_libint(a1, a2)

Get’s the vertex coefficient and conversion coefficients for a pair of atoms given by their atom indices

comp_coulomb_den(**kw)

Computes the dense (square) version of the Hartree kernel

comp_coulomb_pack(**kw)

Computes the packed version of the Hartree kernel

comp_coulomb_sparse(**kw)

Computes the dense (square) version of the Hartree kernel

comp_fci_den(hk, dtype=<class 'numpy.float64'>)

Compute the four-center integrals and return it in a dense storage

comp_moments(dtype=<class 'numpy.float64'>)

Computes the scalar and dipole moments for the all functions in the product basis

dipole_check(**kw)

Our standard minimal check

generate_png_chess_dp_vertex()

Produces pictures of the dominant product vertex a chessboard convention

generate_png_spy_dp_vertex()

Produces pictures of the dominant product vertex in a common black-and-white way

get_ac_vertex_array(dtype=<class 'numpy.float64'>)

Returns the product vertex coefficients as 3d array (dense table)

get_da2cc_den(dtype=<class 'numpy.float64'>)

Returns Conversion Coefficients as dense matrix

get_da2cc_nnz()

Computes the number of non-zero matrix elements in the conversion matrix ac <=> dp

get_da2cc_sparse(dtype=<class 'numpy.float64'>, sparseformat=<class 'scipy.sparse.coo.coo_matrix'>)

Returns Conversion Coefficients as sparse COO matrix

get_dp_vertex_array(dtype=<class 'numpy.float64'>)

Returns the product vertex coefficients as 3d array for dominant products

get_dp_vertex_doubly_sparse(dtype=<class 'numpy.float64'>, sparseformat=<class 'pynao.lsofcsr.lsofcsr_c'>, axis=0)

Returns the product vertex coefficients for dominant products as a one-dimensional array of sparse matrices

get_dp_vertex_doubly_sparse_loops(dtype=<class 'numpy.float64'>, sparseformat=<class 'pynao.lsofcsr.lsofcsr_c'>, axis=0)

Returns the product vertex coefficients for dominant products as a one-dimensional array of sparse matrices, slow version

get_dp_vertex_nnz()

Number of non-zero elements in the dominant product vertex: can be speedup, but…

get_dp_vertex_sparse(dtype=<class 'numpy.float64'>, sparseformat=<class 'scipy.sparse.coo.coo_matrix'>)

Returns the product vertex coefficients as 3d array for dominant products, in a sparse format coo(p,ab) by default

get_dp_vertex_sparse2(dtype=<class 'numpy.float64'>, sparseformat=<class 'scipy.sparse.coo.coo_matrix'>)

Returns the product vertex coefficients as 3d array for dominant products, in a sparse format coo(pa,b) by default

get_dp_vertex_sparse_loops(dtype=<class 'numpy.float64'>, sparseformat=<class 'scipy.sparse.coo.coo_matrix'>)

Returns the product vertex coefficients as 3d array for dominant products, in a sparse format coo(p,ab) slow version

get_nprod()

Number of (atom-centered) products

init_c2s_domiprod()

Compute the array of start indices for dominant product basis set

init_inp_param_prod_log_dp(sv, tol_loc=1e-05, tol_biloc=1e-06, ac_rcut_ratio=1.0, ac_npc_max=8, jcutoff=14, metric_type=2, optimize_centers=0, ngl=96, **kw)

Talman’s procedure should be working well with a pseudo-potential hamiltonians. This subroutine prepares the class for a later atom pair by atom pair generation of the dominant product vertices and the conversion coefficients by calling subroutines from the library libnao.

init_prod_basis_pp(sv, **kvargs)

Talman’s procedure should be working well with Pseudo-Potential starting point.

init_prod_basis_pp_batch(nao, **kw)

Talman’s procedure should be working well with Pseudo-Potential starting point.

ls_contributing(a1, a2)

Get the list of contributing centers

overlap_check(**kw)

Our standard minimal check comparing with overlaps

reshape_COO(shape)

Reshape the sparse matrix (a) and returns a coo_matrix with favourable shape.

pynao.prod_log module

pynao.prod_log.dipole_check(sv, prod_log, dipole_funct=<function dipole_ni>, **kvargs)

Computes the allclose(), mean absolute error and maximal error of the dipoles reproduced by the (local) vertex.

pynao.prod_log.overlap_check(prod_log, overlap_funct=<function overlap_ni>, **kvargs)

Computes the allclose(), mean absolute error and maximal error of the overlap reproduced by the (local) vertex.

class pynao.prod_log.prod_log(**kw)

Bases: pynao.ao_log.ao_log

prod_log object.

Holder of (local) product functions and vertices.

Input parameters:

ao_log: ao_log class (pynao.ao_log) or pyscf Mole object (pyscf.gto.mole)

holder of the numerical orbitals

tol: float

tolerance to keep the linear combinations

Output parameters:

For each specie returns a set of radial functions defining a product basis. These functions are sufficient to represent the products of original atomic orbitals via a product vertex coefficients.

hartree_pot(**kvargs)

Compute Hartree potential of the radial orbitals and return another ao_log_c storage with these potentials.

init_prod_log_df(**kw)

Initializes the radial functions from pyscf

init_prod_log_dp(**kw)

Builds linear combinations of the original orbital products

lambda_check_coulomb()

Check the equality (p|q)<q,cd> = [p,ab] <ab|q>(q|r)<r|cd>

lambda_check_overlap(overlap_funct=<function overlap_am>, **kvargs)

Check the equality (p) = [p,ab] S^ab, i.e. scalar moments are recomputed with inversed vertex from the ao’s overlap

load_prod_log_dp(ao_log, sp2charge, tol_loc=1e-05, fname='local2functs_vertex.hdf5')

Builds linear combinations of the original orbital products Testing the calculations using same kernel than from fortran

overlap_check(overlap_funct=<function overlap_ni>, **kvargs)

Recompute the overlap between orbitals using the product vertex and scalar moments of product functions

pynao.qchem_inter_rf module

class pynao.qchem_inter_rf.qchem_inter_rf(**kw)

Bases: pynao.scf.scf

Quantum-chemical interacting response function

inter_rf(ww)

This delivers the interacting response function in the product basis

kernel_qchem_inter_rf(**kw)

kernel_qchem_inter_rf_pos_neg(**kw)

This is constructing the E_m-E_n and E_n-E_m matrices

pynao.scf module

class pynao.scf.scf(**kw)

Bases: pynao.chi0_matvec.chi0_matvec

add_pb_hk(**kw)

get_fock(h1e=None, dm=None, **kw)

Fock matrix for the given density matrix (matrix elements of the Fockian).

get_hcore(mol=None, **kw)

get_j(dm=None, **kw)

Compute J matrix for the given density matrix (matrix elements of the Hartree potential).

get_jk(mol=None, dm=None, hermi=1, **kw)

get_k(dm=None, **kw)

Compute K matrix for the given density matrix.

get_ovlp(sv=None)

get_veff(dm=None, **kw)

Hartree-Fock potential matrix for the given density matrix

kernel_scf(dump_chk=False, **kw)

This does the actual SCF loop so far only HF

vhartree_coo(**kw)

vhartree_den(**kw)

Compute matrix elements of the Hartree potential and return dense matrix compatible with RHF or UHF

vnucele_coo(**kw)

Compute matrix elements of nuclear-electron interaction (attraction) it subtracts the computed matrix elements from the total Hamiltonian to find out the nuclear-electron interaction in case of SIESTA import. This means that Vne is defined by Vne = H_KS - T - V_H - V_xc

pynao.scf_dos module

pynao.scf_dos.eigenvalues2dos(ksn2e, zomegas, nkpoints=1)

Compute the Density of States using the eigenvalues

pynao.scf_dos.scf_dos(scf, zomegas, nkpoints=1)

Compute the Density of States using the eigenvalues

pynao.scf_dos.system_vars_ados(sv, zomegas, ls_atom_groups, nkpoints=1)

Compute a Partial Density of States (resolved in atomic indices) using the eigenvalues and eigenvectors in wfsx

pynao.tddft_iter module

class pynao.tddft_iter.tddft_iter(**kw)

Bases: pynao.chi0_matvec.chi0_matvec

tddft_iter object.

Iterative TDDFT a la PK, DF, OC JCTC

References:

1. [Barbry2018thesis]

2. [koval2018pynao]

3. [koval2010-mbptlcao]

apply_kernel(dn)

Apply the kernel to chi0_mv, i.e, matrix vector multiplication. The kernel is stored in pack or pack sparse format.

Input parameters:

dn: 1D np.array, complex

the results of chi0_mv product, i.e, chi0*delta_V_ext

Output parameters:

vcre: 1D np.array, float

The real part of the resulting matvec product kernel.dot(chi0_mv)

vcim: 1D np.array, float

The imaginary part of the resulting matvec product kernel.dot(chi0_mv)

comp_polariz_inter_ave(comegas, tmp_fname=None)

Compute average interacting polarizability

comp_polariz_inter_edir(comegas, eext=array([1.0, 1.0, 1.0]), tmp_fname=None)

Compute average interacting polarizability

comp_polariz_inter_xx(comegas, tmp_fname=None)

Compute the interacting polarizability along the xx direction

comp_polariz_nonin_ave(comegas, tmp_fname=None)

Compute average interacting polarizability

comp_polariz_nonin_edir(comegas, eext=array([1.0, 1.0, 1.0]), tmp_fname=None)

Compute average interacting polarizability

comp_polariz_nonin_xx(comegas, tmp_fname=None)

Compute the non-interacting polarizability along the xx direction

comp_veff(vext, comega=0j, prev_sol=None)

This computes an effective field (scalar potential) given the external scalar potential.

Solves

Ax = b

with

• b the external potential delta V_ext

• x the effective potential delta V_eff

• A = (1 - K_{Hxc}Chi_0)

get_eels_spectrum(freq, velec=array([1.0, 0.0, 0.0]), beam_offset=array([0.0, 0.0, 0.0]), dr=array([0.3, 0.3, 0.3]), tmp_fname=None, Vext=None, inter=True)

Calculate the interacting TEM spectra for a given electron trajectory. Typically, can be used to calculate valence Electron Energy Loss Spectroscopy (eels). The perturbation is created by a moving charge. See chapter 6 of [Barbry2018thesis] for more details.

Input Parameters:

freq: 1D np.array, float

Frequency range for which the eels spectra is calculated in atomic units.

velec: 1D np.array, float

xyz component of the electron velocity in atomic unit

beam_offset: 1D np.array, float

xyz components of the beam offset, must be orthogonal to velec in atomic unit

dr: 1D np.array, float

spatial resolution for the electron trajectory in atomic unit. Warning: This parameter influence the accuracy of the calculations. If taken too large the results will be highly inacurate.

tmp_fname: string

filename to store temporary results of the eels spectra. if None, the results will not be saved during the run

Vext: 1D np.array, Complex

Optional, if not given, the external potential created by the charge will be calculated. Can be used to avoid to recalculate Vext if doing calculation where only the norm of the velocity is changing.

inter: bool

Perform interactive or non-interactive calculations.

Output Parameters:

Vext: 1D np.array, complex

The external potential created by the charge in frequency domain. Can be used again for further calculations.

eels: 1D np.array, complex

The calculated eels spectrum

The velocity is given in atomic units, it can be converted from energy (in Hartree) using the relativistic energy kinetic formulation. Below is a function doing such conversion

Convert electron energy (in eV) to velocity:

>>> def E2vel(E):
>>>     # Fine-structure
>>>     alpha = 0.0072973525693
>>>     num = np.sqrt((alpha**2)*E**2 + 2*E)
>>>     den = E*alpha**2 + 1
>>>     return num/den
>>> # Convert electron energy in eV to velocity in atomic unit
>>> Ha = 27.211
>>> E = 1e5/Ha      # E = 100 keV
>>> v = E2vel(E)
>>> print(v)

polariz_inter_ave(comegas, tmp_fname=None)

Compute average interacting polarizability

preconditioning(comega, vin, x0)

vext2veff_matvec(vin)

vext2veff_matvec2(vin)

pynao.tddft_iter_2ord module

class pynao.tddft_iter_2ord.tddft_iter_2ord(**kw)

Bases: pynao.tddft_iter.tddft_iter

tddft_iter_2ord object.

Iterative TDDFT with a high-energy part of the KS eigenvectors compressed

kchi0_mv(v, comega=None)

Operator O = K chi0

polariz_dckcd(comegas)

Compute a term <d chi0 K chi0 d>

polariz_upkc(comegas)

Compute interacting polarizability along the xx direction using an alternative algorighm with chi = chi0 (1+K chi0) [1-K chi0 K chi0]^(-1)

solve_umkckc(vext, comega=0j, x0=None)

This solves a system of linear equations

(1 - K chi0 K chi0 ) X = vext

or computes

X = (1 - K chi0 K chi0 )^{-1} vext

umkckc_mv(v)

Operator O = [1- K chi0 K chi0]

upkc_mv(v)

Operator O = [1 + K chi0]

pynao.tddft_iter_x_zip module

class pynao.tddft_iter_x_zip.tddft_iter_x_zip(**kw)

Bases: pynao.tddft_iter.tddft_iter

Iterative TDDFT with a high-energy part of the KS eigenvectors compressed

build_x_zip()

define compressed eigenvectors

pynao.tddft_tem module

pynao.tddft_tem.calc_external_potential(nprod, velec, beam_offset, freq, freq_sym, time_range, ao_log=None, verbosity=0)

Calculate the external potential created by a moving charge

pynao.tddft_tem.check_collision(velec, beam_offset, atom2coord, threshold=1e-06, verbosity=1)

Check if the electron collide with an atom if the trajectory of the electron is passing too close of an atom, i.e., if at any moment in time the distance between the electron and an atom is below the value given by threshold.

Also check that the electron velocity and the beam ofset are perpendicular.

Input Parameters:

atom2coord: 2D np.array, float

Atoms coordinates

threshold: float

Threshold use to determine if the beam is colliding an atom.

pynao.tddft_tem.comp_tem_spectrum(self, comegas, V_freq, tmp_fname=None, inter=True)

Compute the interacting tem spectrum

Input Parameters:

comegas: 1D np.array, complex

frequency range (in Hartree) for which the polarizability is computed. The imaginary part control the width of the signal. For example, >>> td = tddft_iter_c(…) >>> comegas = np.arange(0.0, 10.05, 0.05) + 1j*td.eps

V_freq: 2D np.array, complex

The perturbation to the system

inter: boolean

Perform TDDFT calculation without or with interaction

tmp_fname: string

temporary file to save polarizability at each frequency. Can be a life saver for large systems. The format of the file is the following, # energy (Hartree) Re(gamma) Im(gamma)

Output Parameters:

dn: 2D np.array, complex

computed density change in prod basis

gamma: 1D np.array, complex

computed eels spectrum

pynao.tddft_tem.get_time_range(velec, beam_offset, dr, freq, verbosity=0)

Get the time and symmetric frequency range for the electron passing close to the particle. The time range is a symmetric array around 0.0. At t = 0, the electron is at its closest position from the molecule. This array will depend on the frequency range and the spatial precision dr. To respect the Fourier transform convention, the following relationshiip must be fulfill,

N = 2*pi/(dw*dt)

with N the number of element of t. N must be an odd number in order that t is symmetric

Input Parameters:

dr: 1D np.array, float, size 3

spatial resolution for the electron trajectory in atomic unit. Warning: This parameter influence the accuracy of the calculations. If taken too large the results will highly inacurate.

freq: 1D np.array, float

Frequency range (in atomic unit), freq[0] must be 0.0!!

Module contents

Numerical Atomic Orbitals