atoMEC package

Subpackages

Submodules

atoMEC.check_inputs module

The check_inputs module checks the validity of all user-defined inputs.

If inputs are invalid, InputError exceptions are raised. It also assigns appropriate default inputs where none are supplied.

Classes

exception atoMEC.check_inputs.InputError[source]

Bases: Exception

Exit atoMEC and print relevant input error message.

ELF_error()[source]

Raise exception if error in ELF inputs.

Parameters:

err_msg (str) – the error message printed

Return type:

None

SCF_error()[source]

Raise exception if error in SCF inputs.

Parameters:

err_msg (str) – the error message printed

Return type:

None

bands_error()[source]

Raise exception if nkpts not positive int or dE_spc not positive number.

Parameters:

err_msg (str) – the error message printed

Return type:

None

bc_error()[source]

Raise exception if unbound not str or in permitted values.

Parameters:

err_msg (str) – the error message printed

Return type:

None

charge_error()[source]

Raise an exception if charge is not an integer.

Parameters:

None –

Return type:

None

conv_error()[source]

Raise exception if error in convergence inputs.

Parameters:

err_msg (str) – the error message printed

Return type:

None

density_error()[source]

Raise an exception if density is not a float or negative.

Parameters:

err_msg (str) – error message printed

Return type:

None

grid_error()[source]

Raise exception if error in grid inputs.

Parameters:

err_msg (str) – the error message printed

Return type:

None

species_error()[source]

Raise an exception if there is an invalid species.

Parameters:

err_msg (str) – error message printed

Return type:

None

spinmag_error()[source]

Raise an exception if density is not a float or negative.

Parameters:

err_msg (str) – error message printed

Return type:

None

spinpol_error()[source]

Raise exception if spinpol not a boolean.

Parameters:

err_msg (str) – the error message printed

Return type:

None

temp_error()[source]

Raise an exception if temperature is not a float.

Parameters:

err_msg (str) – error message printed

Return type:

None

unbound_error()[source]

Raise exception if unbound not str or in permitted values.

Parameters:

err_msg (str) – the error message printed

Return type:

None

v_shift_error()[source]

Raise exception if error in SCF inputs.

Parameters:

err_msg (str) – the error message printed

Return type:

None

xc_error()[source]

Raise an exception if density is not a float or negative.

Parameters:

err_msg (str) – error message printed

Return type:

None

class atoMEC.check_inputs.Atom[source]

Bases: object

Check the inputs from the Atom class.

check_charge(charge)[source]

Check the net charge is an integer.

Parameters:

charge (int) – the net charge

Returns:

charge – the net charge (if input valid)

Return type:

int

Raises:

InputError.charge_error – if charge is not an integer

check_density(density)[source]

Check the mass density is valid and reasonable.

Parameters:

density (float or int) – mass density (in \(\mathrm{g\ cm}^{-3}\))

Returns:

density – mass density (in \(\mathrm{g\ cm}^{-3}\)) if input accepted

Return type:

float

Raises:

InputError.density_error – if the density is not a positive number <= 100

check_rad_dens_init(atom, radius, density, units_radius, units_density)[source]

Check that at least one of radius or density is specified and reasonable.

In case both are specified, check they are compatible.

Parameters:

  • Atom (atoMEC.Atom) – the main Atom object

  • radius (float or int) – Wigner-Seitz radius

  • density (float or int) – mass density

  • units_radius (str) – units of radius

  • units_density (str) – units of density

Returns:

radius, density – the Wigner-Seitz radius and mass density if inputs are valid

Return type:

tuple of floats

Raises:

InputError.density_error – if neither density nor radius is not given, or if one is invalid, or if both are given and they are incompatible

check_radius(radius, units_radius)[source]

Check the Wigner-Seitz radius is valid and reasonable.

Parameters:

  • radius (float or int) – Wigner-Seitz radius (in input units)

  • units_radius (str) – input units of radius

Returns:

radius – Wigner-Seitz radius in Hartree units (Bohr)

Return type:

float

Raises:

InputError.density_error – if the radius is not a positive number > 0.1

check_species(species)[source]

Check the species is a string and corresponds to an actual element.

Parameters:

species (str) – chemical symbol for atomic species

Return type:

None

Raises:

InputError.species_error – Chemical symbol is not valid

check_temp(temp, units_temp)[source]

Check the temperature is a float within a sensible range.

Parameters:

  • temp (float) – temperature (in any accepted units)

  • units_temp (str) – units of temperature

Returns:

temp – temperature in units of Hartree

Return type:

float

Raises:

  • InputError.temp_error – input temperature is not a positive number

  • InputWarning.temp_warning – input temperature is not inside a well-tested range

check_units_density(units_density)[source]

Check the units of density are accepted.

Parameters:

units_density (str) – units of density

Returns:

units_density – units of density (if accepted) converted to lowercase

Return type:

str

Raises:

InputError.density_error – if units of density are not one of “g/cm3” or “gcm3”

check_units_radius(units_radius)[source]

Check the units of radius are accepted.

Parameters:

units_radius (str) – units of radius

Returns:

units_radius – units of radius (if accepted) converted to lowercase

Return type:

str

Raises:

InputError.density_error – if units of radius are not one of “bohr”, “angstrom” or “ang”

check_units_temp(units_temp)[source]

Check the units of temperature are accepted.

Parameters:

units_temp (str) – units of temperature

Returns:

units_temp – units of temperature (if valid input) converted to lowercase

Return type:

str

Raises:

InputError.temp_error – unit of temperature is not accepted, i.e. not one of “ha”, “ev” or “k”

dens_to_radius(atom, density)[source]

Convert the mass density to a Wigner-Seitz radius.

Parameters:

  • atom (atoMEC.Atom) – the main Atom object

  • density (float) – the mass density

Returns:

radius – the Wigner-Seitz radius

Return type:

float

radius_to_dens(atom, radius)[source]

Convert the Voronoi sphere radius to a mass density.

Parameters:

  • atom (atoMEC.Atom) – the main Atom object

  • radius (float) – the Wigner-Seitz radius

Returns:

density – the mass density

Return type:

float

class atoMEC.check_inputs.EnergyCalcs[source]

Bases: object

Check inputs for CalcEnergy calculations.

static check_band_params(input_params)[source]

Check if band parameters are reasonable, or assign if empty.

Parameters:

input_params (dict) – can contain the keys maxscf and mixfrac for max scf cycle and potential mixing fraction

Returns:

band_params – dictionary for band parameters as follows: { nkpts (int) : number of levels per band, de_min (float) : minimum energy gap to make a band }

Return type:

dict

Raises:

InputError.bands_error – if band parameters are of invalid type or range

static check_conv_params(input_params)[source]

Check convergence parameters are reasonable, or assigns if empty.

Parameters:

input_params (dict of floats) – Can contain the keys econv, nconv and vconv, for energy, density and potential convergence parameters

Returns:

conv_params – dictionary of convergence parameters as follows: { econv (float) : convergence for total energy, nconv (float) : convergence for density, vconv (float) : convergence for electron number, eigtol (float) : tolerance for eigenvalues, bandtol (float) : tolerance for n(l)max warning }

Return type:

dict of floats

Raises:

InputError.conv_error – if a convergence parameter is invalid (not float or negative).

static check_grid_params(grid_params)[source]

Check grid parameters are reasonable, or assigns if empty.

Parameters:

grid_params (dict) – Can contain the keys ngrid (int, number of grid points), x0 (float, LHS grid point for log grid), and s0 (float, LHS grid point for sqrt grid)

Returns:

grid_params – dictionary of grid parameters as follows: { ngrid (int) : number of grid points, x0 (float) : LHS grid point takes form \(r_0=\exp(x_0)\); \(x_0\) can be specified s0 (float) : LHS grid point takes form \(r_0=s0^2\); \(s_0\) can be specified ngrid_coarse (int) : (smaller) number of grid points for estimation of eigenvalues with full diagonalization }

Return type:

dict

Raises:

  • InputError.grid_error – if grid inputs are invalid or outside a reasonable range

  • InputError.ngrid_warning – if ngrid is outside a reasonable convergence range

static check_grid_type(grid_type)[source]

Check grid type.

Parameters:

grid_type (str) – the grid type

Returns:

grid_type – the grid type

Return type:

str

Raises:

InputError.grid_type_error – if grid type not one of “log” or “sqrt”

static check_scf_params(input_params)[source]

Check scf parameters are reasonable, or assigns if empty.

Parameters:

input_params (dict) – can contain the keys maxscf and mixfrac for max scf cycle and potential mixing fraction

Returns:

scf_params – dictionary with the following keys: { maxscf (int) : max number scf cycles, mixfrac (float) : mixing fraction }

Return type:

dict

Raises:

InputError.SCF_error – if the SCF parameters are not of correct type or in valid range

class atoMEC.check_inputs.ISModel[source]

Bases: object

Check the inputs for the IS model class.

calc_nele(nele, spinpol)[source]

Calculate the electron number in each spin channel (if spin polarized).

Parameters:

  • spinmag (int) – the spin magnetization

  • nele (int) – total electron number

  • spinpol (bool) – spin polarization

Returns:

nele – number of electrons in each spin channel if spin-polarized, else just total electron number

Return type:

array of ints

check_bc()[source]

Check the boundary condition is accepted.

Parameters:

bc (str) – the boundary condition used to solve KS eqns (can be either “dirichlet” or “neumann”)

Returns:

bc – the boundary condition used to solve KS eqns (lowercase)

Return type:

str

Raises:

InputError.bc_error – if the boundary condition is not recognised

check_spinmag(nele)[source]

Check the spin magnetization is compatible with the total electron number.

Also compute a default value if none is specified.

Parameters:

  • spinmag (int) – the spin magnetization (e.g. 1 for a doublet state)

  • nele (int) – the total number of electrons

Returns:

spinmag – the spin magnetization if input valid

Return type:

int

Raises:

InputError.spinmag_error – if spinmag input is not an integer or incompatible with electron number

check_spinpol()[source]

Check the spin polarization is a boolean.

Parameters:

spinpol (bool) – whether spin polarized calculation is done

Returns:

spinpol – same as input unless error raised

Return type:

bool

Raises:

InputError.spinpol_error – if the spin polarization is not a bool

check_unbound(bc)[source]

Check the unbound electron input is accepted.

Parameters:

  • unbound (str) – defines the treatment of the unbound electrons

  • bc (str) – the boundary condition

Returns:

unbound – treatment of unbound electrons (if input valid)

Return type:

str

Raises:

InputError.unbound_error – if the treatment of unbound electrons is not a valid input

check_v_shift()[source]

Check the potential shift is a boolean.

Parameters:

v_shift (bool) – whether potential is shifted or not

Returns:

  • v_shift (bool)

  • same as input unless error raised

Raises:

InputError.v_shift_error – if the potential shift is not a bool

check_xc(xc_type)[source]

Check the exchange and correlation functionals are accepted.

Parameters:

  • xc_func (str or int) – the libxc name or id of the x/c functional

  • xc_type (str) – type i.e. “exchange” or “correlation”

Returns:

xc_func – the libxc name of the x/c functional (if valid input)

Return type:

str

Raises:

InputError.xc_error – if xc functional is not a valid libxc input or is not supported by the current version of atoMEC

class atoMEC.check_inputs.InputWarning[source]

Bases: object

Warn if inputs are considered outside of typical ranges, but proceed anyway.

ngrid_warning(err2)[source]

Warn if grid params outside of sensible range.

Parameters:

  • err1 (str) – first custom part of error message

  • err2 (str) – second custom part of error message

Returns:

warning – full error message

Return type:

str

norbs_warning()[source]

Warn if either nmax or lmax is too low for convergence.

Parameters:

quantum_num (str) – lmax or nmax

Returns:

warning – full error message

Return type:

str

temp_warning()[source]

Warn if temperature outside of sensible range.

Parameters:

err (str) – custom part of error message

Returns:

warning – full errror message

Return type:

str

atoMEC.config module

Configuration file to store global parameters.

atoMEC.convergence module

Contains classes and functions used to compute and store aspects related to convergence.

So far, the only procedure requring convergence is the static SCF cycle. More will be added in future.

Classes

  • SCF : holds the SCF convergence attributes and calculates them for the given cycle

class atoMEC.convergence.SCF(xgrid, grid_type)[source]

Bases: object

Convergence attributes and functions related to SCF energy procedures.

Notes

Contains the private attributes _energy, _potential and _density

check_conv(E_free, pot, dens, iscf)[source]

Compute and check the changes in energy, integrated density and integrated potential.

If the convergence tolerances are all simultaneously satisfied, the complete variable returns True as the SCF cycle is complete.

Parameters:

  • E_free (float) – the total free energy

  • v_s (ndarray) – the KS potential

  • dens (ndarray) – the electronic density

  • iscf (int) – the iteration number

Returns:

conv_vals – Dictionary of convergence parameters as follows: { conv_energy : float, conv_rho : ndarray, conv_pot : ndarray, complete : bool }

Return type:

dict

atoMEC.mathtools module

Low-level module containing miscalleneous mathematical functions.

Functions

atoMEC.mathtools.chem_pot(orbs)[source]

Determine the chemical potential by enforcing charge neutrality (see notes).

Uses scipy.optimize.root_scalar with brentq implementation.

Parameters:

orbs (staticKS.Orbitals) – the orbitals object

Returns:

mu – chemical potential (spin-dependent)

Return type:

array_like

Notes

Finds the roots of equation

\[\sum_{nl} (2l+1) f_{fd}(\epsilon_{nl},\beta,\mu) + N_{ub}(\beta,\mu) - N_e = 0.\]

The number of unbound electrons \(N_{ub}\) depends on the implementation choice.

atoMEC.mathtools.f_root_id(mu, eigvals, lbound, nele)[source]

Functional input for the chemical potential root finding function (ideal approx).

See notes for function returned, the ideal approximation is used for free electrons.

Parameters:

  • mu (array_like) – chemical potential

  • eigvals (ndarray) – the energy eigenvalues

  • lbound (ndarray) – the lbound matrix \((2l+1)\Theta(\epsilon_{nl}^\sigma)\)

  • nele (union(int, float)) – the number of electrons for given spin

Returns:

f_root – the difference of the predicted electron number with given mu and the actual electron number

Return type:

float

Notes

The returned function is

\[f = \sum_{nl} (2l+1) f_{fd}(\epsilon_{nl},\beta,\mu) + N_{ub}(\beta,\mu) - N_e\]
atoMEC.mathtools.f_root_qu(mu, eigvals, occ_weight, nele)[source]

Functional input for the chemical potential root finding function (ideal approx).

See notes for function returned, the ideal approximation is used for free electrons.

Parameters:

  • mu (array_like) – chemical potential

  • eigvals (ndarray) – the energy eigenvalues

  • lbound (ndarray) – the lbound matrix \((2l+1)\Theta(\epsilon_{nl}^\sigma)\)

  • nele (union(int, float)) – the number of electrons for given spin

Returns:

f_root – the difference of the predicted electron number with given mu and the actual electron number

Return type:

float

Notes

The returned function is

\[f = \sum_{nl} (2l+1) f_{fd}(\epsilon_{nl},\beta,\mu) + N_{ub}(\beta,\mu) - N_e\]
atoMEC.mathtools.fd_int_complete(mu, beta, n)[source]

Compute complete Fermi-Dirac integral for given order (see notes for function form).

Parameters:

  • mu (float) – chemical potential

  • beta (float) – inverse temperature

  • n (int) – order of Fermi-Dirac integral (see notes)

Returns:

I_n – the complete fermi-dirac integral

Return type:

float

Notes

Complete Fermi-Dirac integrals are of the form

\[I_{n}(\mu,\beta)=\int_0^\infty\mathrm{d}\epsilon\ \epsilon^{n/2}f_{fd} (\mu,\epsilon,\beta)\]

where n is the order of the integral

atoMEC.mathtools.fermi_dirac(eps, mu, beta, n=0)[source]

Compute the Fermi-Dirac function, see notes for functional form.

Parameters:

  • mu (array_like) – the chemical potential

  • beta (float) – the inverse potential

  • eps (array_like) – the energies

  • n (int) – energy is raised to power n/2 in the numerator (see notes)

Returns:

f_fd – the fermi dirac occupation(s)

Return type:

array_like

Notes

The FD function is defined as:

\[f^{(n)}_{fd}(\epsilon, \mu, \beta) = \frac{\epsilon^{n/2}}{1+\exp(1+ \beta(\epsilon - \mu))}\]
atoMEC.mathtools.grad_func(den, xgrid, grid_type)[source]

Compute the gradient of a function on the log / sqrt grid.

Parameters:

  • den (ndarray) – density array or any other function that is integrated

  • xgrid (ndarray) – log / sqrt grid

  • grid_type (str) – grid type, log or sqrt

Returns:

grad – The gradient of the density w.r.t. the radial grid.

Return type:

ndarray

atoMEC.mathtools.ideal_entropy(eps, mu, beta, n=0)[source]

Define the integrand to be used in ideal_entropy_int() (see notes).

Parameters:

  • eps (array_like) – the energies

  • mu (array_like) – the chemical potential

  • beta (float) – the inverse potential

  • n (int) – energy is raised to power n/2 in the numerator (see notes)

Returns:

f_ent – the entropy integrand function

Return type:

array_like

Notes

The ideal entropy integrand is defined as

\[f_n(\epsilon,\mu,\beta) = \epsilon^{n/2} (f_\mathrm{fd}\log{f_\mathrm{fd}} + (1-f_\mathrm{fd}) \log(1-f_\mathrm{fd}) ),\]

where \(f_\mathrm{fd}=f_\mathrm{fd}(\epsilon,\mu,\beta)\) is the Fermi-Dirac distribution.

atoMEC.mathtools.ideal_entropy_int(mu, beta, n)[source]

Compute the entropy for the ideal electron gas (without prefactor) - see notes.

Parameters:

  • mu (float) – chemical potential

  • beta (float) – inverse temperature

  • n (int) – order of Fermi-Dirac integral (see notes)

Returns:

I_n – the complete fermi-dirac integral

Return type:

float

Notes

The entropy of an ideal electron gas is defined as

\[I_n(\mu,\beta) = \int_0^\infty \mathrm{d}\epsilon\ \epsilon^{n/2} (f_\mathrm{fd}\log{f_\mathrm{fd}} + (1-f_\mathrm{fd}) \log(1-f_\mathrm{fd}) ),\]

where \(f_\mathrm{fd}=f_\mathrm{fd}(\epsilon,\mu,\beta)\) is the Fermi-Dirac distribution.

atoMEC.mathtools.int_sphere(fx, xgrid, grid_type)[source]

Compute integral over sphere defined by input grid.

Parameters:

  • fx (array_like) – The function (array) to be integrated

  • xgrid (ndarray) – The log / sqrt radial grid

Returns:

I_sph – The value of the integrand

Return type:

float

Notes

The integral formula on the log grid is given by

\[I = 4 \pi \int \mathrm{d}x\ e^{3x} f(x)\]

and on the sqrt grid is given by

\[I = 4 \pi \int \mathrm{d}s\ s^5 f(s)\]
atoMEC.mathtools.laplace(y, x, axis=-1)[source]

Compute the second-order derivative \(d^2 y(x) / dx^2\).

Derivative can be computed over any given axis.

Parameters:

  • y (ndarray) – array y(x) on which laplacian is computed

  • x (ndarray) – x array

  • axis (int, optional) – axis over which derivatives are taken default : -1

Returns:

grad2_y – the laplacian of y

Return type:

ndarray

atoMEC.mathtools.lorentzian(x, x0, gamma)[source]

Compute the Lorentzian function.

Parameters:

  • x (array_like) – the “x” variable

  • x0 (array_like) – the position of the peak

  • gamma (float) – half the peak width / smoothing factor

Returns:

lorentzian_ – the lorentzian function

Return type:

array_like

Notes

The Lorentzian function is defined as

\[\mathcal{L}(x, x_0, \gamma)=\frac{\gamma}{\pi}\frac{1}{\gamma^2+(x-x_0)^2}\]
atoMEC.mathtools.normalize_orbs(eigfuncs_x, xgrid, grid_type)[source]

Normalize the KS orbitals within the chosen sphere.

Parameters:

  • eigfuncs_x (ndarray) – The radial KS eigenfunctions \(X_{nl}^{\sigma}(x)\)

  • xgrid (ndarray) – The log / sqrt grid over which normalization is performed

Returns:

eigfuncs_x_norm – The radial KS eigenfunctions normalized over the chosen sphere

Return type:

ndarray

atoMEC.models module

Contains models used to compute properties of interest from the Atom object.

So far, the only model implemented is the ISModel. More models will be added in future releases.

Classes

  • ISModel : Ion-sphere type model, static properties such as KS orbitals, density and energy are directly computed

class atoMEC.models.ISModel(atom, xfunc_id='lda_x', cfunc_id='lda_c_pw', bc='bands', spinpol=False, spinmag=-1, unbound='quantum', v_shift=True, write_info=True)[source]

Bases: object

The ISModel represents a particular family of AA models known as ion-sphere models.

The implementation in atoMEC is based on the model described in [1].

Parameter inputs for this model are related to particular choices of approximation, e.g. boundary conditions or exchange-correlation functional, rather than fundamental physical properties.

Parameters:

  • atom (atoMEC.Atom) – The main atom object

  • xfunc_id (str or int, optional) – The exchange functional, can be the libxc code or string, or special internal value Default : “lda_x”

  • cfunc_id (str or int, optional) – The correlation functional, can be the libxc code or string, or special internal value Default : “lda_c_pw”

  • bc (str, optional) – The boundary condition, can be “dirichlet” or “neumann” Default : “dirichlet”

  • spinpol (bool, optional) – Whether to run a spin-polarized calculation Default : False

  • spinmag (int, optional) – The spin-magentization Default: 0 for nele even, 1 for nele odd

  • unbound (str, optional) – The way in which the unbound electron density is computed Default : “ideal”

  • v_shift (bool, optional) – Shifts the potential vertically by its value at the boundary, v(R_s) Default : True

  • write_info (bool, optional) – Writes information about the model parameters Default : True

nele_tot

total number of electrons

Type:

int

nele

number of electrons per spin channel (or total if spin unpolarized)

Type:

array_like

References

CalcEnergy(nmax, lmax, grid_params={}, conv_params={}, scf_params={}, band_params={}, force_bound=[], grid_type='sqrt', verbosity=0, write_info=True, write_density=True, write_potential=True, write_eigs_occs=True, write_dos=True, density_file='density.csv', potential_file='potential.csv', eigs_occs_file='eigs_occs.csv', dos_file='dos.csv', guess=False, guess_pot=0)[source]

Run a self-consistent calculation to minimize the Kohn-Sham free energy.

Parameters:

  • nmax (int) – maximum no. eigenvalues to compute for each value of angular momentum

  • lmax (int) – maximum no. angular momentum eigenfunctions to consider

  • grid_params (dict, optional) – dictionary of grid parameters as follows: { ngrid (int) : number of grid points, x0 (float) : LHS grid point takes form \(r_0=\exp(x_0)\); \(x_0\) can be specified, ngrid_coarse (int) : (smaller) number of grid points for estimation of eigenvalues with full diagonalization }

  • conv_params (dict, optional) – dictionary of convergence parameters as follows: { econv (float) : convergence for total energy, nconv (float) : convergence for density, vconv (float) : convergence for electron number, eigtol (float) : convergence for eigenvalues, bandtol (float) : tolerance for n(l)max warning }

  • scf_params (dict, optional) – dictionary for scf cycle parameters as follows: { maxscf (int) : maximum number of scf cycles, mixfrac (float) : density mixing fraction }

  • band_params (dict, optional) – dictionary for band parameters as follows: { nkpts (int) : number of levels per band, de_min (float) : minimum energy gap to make a band }

  • force_bound (list of list of ints, optional) – force certain levels to be bound, for example: force_bound = [0, 1, 0] forces the orbital with quantum numbers \(\sigma=0,\ l=1,\ n=0\) to be always bound even if it has positive energy. This prevents convergence issues.

  • grid_type (str, optional) – the transformed radial grid used for the KS equations. can be ‘sqrt’ (default) or ‘log’

  • verbosity (int, optional) – how much information is printed at each SCF cycle. verbosity=0 prints the total energy and convergence values (default). verbosity=1 prints the above and the KS eigenvalues and occupations.

  • write_info (bool, optional) – prints the scf cycle and final parameters default: True

  • write_density (bool, optional) – writes the density to a text file default: True

  • write_potential (bool, optional) – writes the potential to a text file default: True

  • density_file (str, optional) – name of the file to write the density to default: density.csv

  • potential_file (str, optional) – name of the file to write the potential to default: potential.csv

  • eigs_occs_file (str, optional) – filename for the orbital energies and occupations default: eigs_occs

  • dos_file (str, optional) – filename for the density-of-states and fd distribution default: dos

  • guess (bool, optional) – use coulomb pot (guess=False) or given pot (guess=True) as initial guess

  • guess_pot (numpy array, optional) – initial guess for the potential

Returns:

output_dict – dictionary containing final KS quantities as follows: { energy (staticKS.Energy) : total energy object, density (staticKS.Density) : density object, potential (staticKS.Potential) : potential object, orbitals (staticKS.Orbitals) : orbitals object }

Return type:

dict

CalcPressure(atom, energy_output, nmax=None, lmax=None, conv_params={}, scf_params={}, band_params={}, force_bound=[], write_info=False, verbosity=0, dR=0.01, method='A')[source]

Calculate the electronic pressure using the finite differences method.

N.B.: This is just a wrapper for the postprocess.pressure.finite_diff function. It is maintained to support backwards compatibility until next major release.

Parameters:

  • atom (atoMEC.Atom) – The main atom object

  • energy_output (dict) – output parameters of the function CalcEnergy

  • conv_params (dict, optional) – dictionary of convergence parameters as follows: { econv (float) : convergence for total energy, nconv (float) : convergence for density, vconv (float) : convergence for electron number, eigtol (float) : tolerance for eigenvalues }

  • scf_params (dict, optional) – dictionary for scf cycle parameters as follows: { maxscf (int) : maximum number of scf cycles, mixfrac (float) : density mixing fraction }

  • force_bound (list of list of ints, optional) – force certain levels to be bound, for example: force_bound = [0, 1, 0] forces the orbital quith quantum numbers \(\sigma=0,\ l=1,\ n=0\) to be always bound even if it has positive energy. This prevents convergence issues.

  • verbosity (int, optional) – how much information is printed at each SCF cycle. verbosity=0 prints the total energy and convergence values (default) verbosity=1 prints the above and the KS eigenvalues and occupations.

  • write_info (bool, optional) – prints the scf cycle and final parameters defaults to False

  • dR (float, optional) – radius difference for finite difference calculation defaults to 0.01

Returns:

P_e – electronic pressure in Ha

Return type:

float

property bc

boundary condition for solving the KS equations in a finite sphere.

Type:

str

property cfunc_id

correlation functional shorthand id.

Type:

str

property info

formatted description of the ISModel attributes.

Type:

str

property spinmag

the spin magentization (difference in no. up/down spin electrons).

Type:

int

property spinpol

Whether calculation will be spin-polarized.

Type:

bool

property unbound

the treatment of unbound/free electrons.

Type:

str

property v_shift

shift the potential vertically by a constant.

Type:

bool

property xfunc_id

exchange functional shorthand id.

Type:

str

atoMEC.numerov module

The numerov module handles the routines required to solve the KS equations.

So far, there is only a single implementation which is based on a matrix diagonalization. There is an option for parallelizing over the angular quantum number l.

Classes

  • Solver : Solve the radial KS equation via matrix diagonalization of Numerov’s method. Can use either logarathmic or sqrt grid.

class atoMEC.numerov.Solver(grid_type)[source]

Bases: object

Class holding the numerov solver for KS equations.

Parameters:

grid_type (str) – The grid type (sqrt or log)

KS_matsolve_parallel(T_sparse, B_sparse, v, xgrid, bc, solve_type, eigs_min_guess)[source]

Solve the KS matrix diagonalization by parallelizing over config.numcores.

Parameters:

  • T (ndarray) – kinetic energy array

  • B (ndarray) – off-diagonal array (for RHS of eigenvalue problem)

  • v (ndarray) – KS potential array

  • xgrid (ndarray) – the numerical grid

  • solve_type (str) – whether to do a “full” or “guess” calculation: “guess” estimates the lower bounds of the eigenvalues

  • eigs_min_guess (ndarray) – input guess for the lowest eigenvalues for given l and spin, should be dimension config.spindims * config.lmax

Returns:

  • eigfuncs (ndarray) – radial KS wfns

  • eigvals (ndarray) – KS eigenvalues

Notes

The parallelization is done via the joblib.Parallel class of the joblib library, see here for more information.

For “best” performance (i.e. exactly one core per call of the diagonalization routine plus one extra for the “master” node), the number of cores should be chosen as config.numcores = 1 + config.spindims * config.lmax. However, if this number is larger than the total number of cores available, performance is hindered.

Therefore for “good” performance, we can suggest: config.numcores = max(1 + config.spindimgs * config.lmax, n_avail), where n_avail is the number of cores available.

The above is just a guide for how to choose config.numcores, there may well be better choices. One example where it might not work is for particularly large numbers of grid points, when the memory required might be too large for a single core.

N.B. if config.numcores=-N then joblib detects the number of available cores n_avail and parallelizes into n_avail + 1 - N separate jobs.

KS_matsolve_serial(T_sparse, B_sparse, v, xgrid, bc, solve_type, eigs_min_guess)[source]

Solve the KS equations via matrix diagonalization in serial.

Parameters:

  • T (ndarray) – kinetic energy array

  • B (ndarray) – off-diagonal array (for RHS of eigenvalue problem)

  • v (ndarray) – KS potential array

  • xgrid (ndarray) – the numerical grid

  • solve_type (str) – whether to do a “full” or “guess” calculation: “guess” estimates the lower bounds of the eigenvalues

  • eigs_min_guess (ndarray) – input guess for the lowest eigenvalues for given l and spin, should be dimension config.spindims * config.lmax

Returns:

  • eigfuncs (ndarray) – radial KS wfns

  • eigvals (ndarray) – KS eigenvalues

calc_eigs_min(v, xgrid, bc)[source]

Compute an estimate for the minimum values of the KS eigenvalues.

This estimate uses full diagonalization of the Hamiltonian on a coarse grid. The eigenvalue estimates are used in Scipy’s sparse eigenvalue solver (for the full grid) which optimizes the performance of the solver.

Parameters:

  • v (ndarray) – the KS potential

  • xgrid (ndarray) – the numerical grid (full)

Returns:

eigs_min – array containing estimations of the lowest eigenvalue for each value of l angular quantum number and spin quantum number

Return type:

ndarray

calc_wfns_e_grid(xgrid, v, e_arr, eigfuncs_l, eigfuncs_u)[source]

Compute all KS orbitals defined on the energy grid.

This routine is used to propagate a set of orbitals defined with a fixed set of energies. It is used for the bands boundary condition in the models.ISModel class.

Parameters:

  • xgrid (ndarray) – the spatial (numerical) grid

  • v (ndarray) – the KS potential

  • e_arr (ndarray) – the energy grid

Returns:

eigfuncs_e – the KS orbitals with defined energies

Return type:

ndarray

diag_H(p, T_sparse, B_sparse, v, xgrid, nmax, bc, eigs_guess, solve_type)[source]

Diagonilize the Hamiltonian for the input potential v[p].

Uses Scipy’s sparse matrix solver scipy.sparse.linalg.eigs. This searches for the lowest magnitude nmax eigenvalues, so care must be taken to converge calculations wrt nmax.

Parameters:

  • p (int) – the desired index of the input array v to solve for

  • T (ndarray) – the kinetic energy matrix

  • B (ndarray) – the off diagonal matrix multiplying V and RHS

  • v (ndarray) – KS potential array

  • xgrid (ndarray) – the numerical grid

  • nmax (int) – number of eigenvalues returned by the sparse matrix diagonalization

  • bc (str) – the boundary condition

Returns:

  • evecs (ndarray) – the KS radial eigenfunctions

  • evals (ndarray) – the KS eigenvalues

static mat_convert_dirichlet(diag_mat)[source]

Truncate a scipy diags matrix from (N,N) to (N-1, N-1).

Parameters:

diag_mat (scipy.sparse.diags object) – diagonal matrix to truncate

Returns:

diag_mat – truncated diagonal matrix

Return type:

scipy.sparse.diags object

matrix_solve(v, xgrid, bc, solve_type='full', eigs_min_guess=None)[source]

Solve the radial KS equation via matrix diagonalization of Numerov’s method.

See notes for details of the implementation.

Parameters:

  • v (ndarray) – the KS potential on the log/sqrt grid

  • xgrid (ndarray) – the numerical grid

  • solve_type (str, optional) – whether to do a “full” or “guess” calculation: “guess” estimates the lower bounds of the eigenvalues

  • eigs_min_guess (ndarray, optional) – input guess for the lowest eigenvalues for given l and spin, should be dimension config.spindims * config.lmax

Returns:

  • eigfuncs (ndarray) – the radial KS eigenfunctions on the log/sqrt grid

  • eigvals (ndarray) – the KS eigenvalues

Notes

The implementation is based on [2].

The matrix diagonalization is of the form:

\[\begin{split}\hat{H} \lvert X \rangle &= \lambda \hat{B} \lvert X \rangle\ , \\ \hat{H} &= \hat{T} + \hat{B}\times\hat{V}\ , \\ \hat{T} &= -0.5\times\hat{p}\times\hat{A}\ .\end{split}\]

where \(\hat{p}=\exp(-2x)\). See [2] for definitions of the matrices \(\hat{A}\) and \(\hat{B}\).

References

static num_propagate(xgrid, W, e_arr, eigfuncs_init, grid_type)[source]

Propagate the wfn manually for fixed energy with numerov scheme.

Parameters:

  • xgrid (ndarray) – the numerical grid

  • W (ndarray) – flattened potential array (with angular momentum term)

  • e_arr (ndarray) – flattened energy array

  • eigfuncs_init (ndarray) – initial values for eigenfunctions

Returns:

Psi_norm – normalized wavefunction

Return type:

ndarray

static update_orbs(l_eigfuncs, l_eigvals, xgrid, bc, K, grid_type)[source]

Sort the eigenvalues and functions by ascending energies and normalize orbs.

Parameters:

  • l_eigfuncs (ndarray) – input (unsorted and unnormalized) eigenfunctions (for given l and spin)

  • l_eigvals (ndarray) – input (unsorted) eigenvalues (for given l and spin)

  • xgrid (ndarray) – the numerical grid

  • bc (str) – the boundary condition

Returns:

  • eigfuncs (ndarray) – sorted and normalized eigenfunctions

  • eigvals (ndarray) – sorted eigenvalues in ascending energy

atoMEC.staticKS module

The staticKS module handles time-indepenent KS quantities such as the KS orbitals.

This module is a ‘workhouse’ for much of the average-atom functionality. Computation of static KS quantities is designed to be quite general. The main assumption used everywhere is that we have a logarithmically spaced grid.

Classes

  • Orbitals : Holds the KS orbitals (transformed to the log grid) and related quantities such as eigenvalues and occupation numbers

  • Density : Holds the KS density (bound/unbound components) and routines to compute it.

  • Potential : Holds the KS potential (split into individual components) and the routines to compute them.

  • Energy : Holds the free energy and all internal components (KS quantities and entropy) and the routines required to compute them.

  • EnergyAlt : Holds the free energy and all internal components (KS quantities and entropy) and the routines required to compute them. N.B. the EnergyAlt class constructs the energy functional in an alternative manner to the main Energy class.

  • GramSchmidt : Holds the Gram-Schmidt orthoganalization procedure

Functions

  • log_grid() : Sets up the logarithmic (and associated real-space) grid on which all computations rely.

class atoMEC.staticKS.Density(orbs)[source]

Bases: object

Class holding the static KS density and routines required to compute its components.

Parameters:

orbs (Orbitals) – The KS orbitals object

static construct_rho_orbs(eigfuncs, occnums, xgrid)[source]

Construct a density from a set of discrete KS orbitals.

Parameters:

  • eigfuncs (ndarray) – the radial eigenfunctions on the xgrid

  • occnums (ndarray) – the orbital occupations

  • xgrid (ndarray) – the log / sqrt grid

Returns:

dens – contains the keys rho and N denoting the density and number of electrons respectively

Return type:

dict of ndarrays

static construct_rho_unbound()[source]

Construct the unbound part of the density.

Parameters:

  • orbs (ndarray) – the radial eigenfunctions on the xgrid

  • xgrid (ndarray) – the log / sqrt grid

Returns:

rho_unbound – contains the keys rho and N denoting the unbound density and number of unbound electrons respectively

Return type:

dict of ndarrays

property MIS

the mean ionization state.

Type:

ndarray

property bound

Bound part of KS density.

Contains the keys rho and N denoting the bound density and number of bound electrons respectively.

Type:

dict of ndarrays

property total

Total KS density \(n(r)\) or \(n(x)\).

Type:

ndarray

property unbound

Unbound part of KS density.

Contains the keys rho and N denoting the unbound density and number of unbound electrons respectively

Type:

dict of ndarrays

class atoMEC.staticKS.Energy(orbs, dens)[source]

Bases: object

Class holding information about the KS total energy and relevant routines.

static calc_E_en(density, xgrid, grid_type)[source]

Compute the electron-nuclear energy.

Parameters:

  • density (ndarray) – the (spin) KS density

  • xgrid (ndarray) – the log / sqrt grid

Returns:

E_en – the electron-nuclear energy

Return type:

float

Notes

Electron-nuclear energy is given by

\[E_{en} = 4\pi\int_0^{r_s} \mathrm{d}{r} r^2 n(r) v_{en}(r),\]

where \(r_s\) denotes the radius of the sphere of interest

static calc_E_ha(density, xgrid, grid_type)[source]

Compute the Hartree energy.

Parameters:

  • density (ndarray) – the (spin) KS density

  • xgrid (ndarray) – the log / sqrt grid

Returns:

E_ha – the Hartree energy

Return type:

float

Notes

The Hartree energy is given by

\[E_\mathrm{ha} = 2\pi\int_0^{r_s}\mathrm{d}r r^2 n(r) v_\mathrm{ha}(r),\]

where \(v_\mathrm{ha}(r)\) is the Hartree potential

calc_E_kin(orbs, xgrid)[source]

Compute the kinetic energy.

Kinetic energy is in general different for bound and unbound components. This routine is a wrapper calling the respective functions.

Parameters:

  • orbs (Orbitals) – the KS orbitals object

  • xgrid (ndarray) – the log / sqrt grid

Returns:

E_kin – Contains tot, bound and unbound keys

Return type:

dict of ndarrays

static calc_E_kin_dens(eigfuncs, occnums, xgrid, grid_type, method='A')[source]

Calculate the local kinetic energy density (KED).

There are multiple definitions in the literature of the local KED, see notes.

Parameters:

  • eigfuncs (ndarray) – the radial KS orbitals on the log grid

  • occnums (ndarray) – the orbital occupations

  • xgrid (ndarray) – the log / sqrt grid

  • method (str, optional) – the definition used for KED, can be ‘A’ or ‘B’ (see notes).

Returns:

e_kin_dens – the local kinetic energy density

Return type:

ndarray

Notes

The methods ‘A’ and ‘B’ in this function are given according to the definitions in [3] and [4].

Both methods should integrate to the same kinetic energy, in the case of Neumann or Diriclet boundary conditions; for the bands condition they will be different. The definition ‘B’ is the one used in the usual definition of the electron localization function [5]. It is given by formula (B.8) in [11].

References

static calc_E_kin_orbs(eigfuncs, occnums, xgrid, grid_type)[source]

Compute the kinetic energy contribution from discrete KS orbitals.

Parameters:

  • eigfuncs (ndarray) – the radial KS orbitals on the log grid

  • occnums (ndarray) – the orbital occupations

  • xgrid (ndarray) – the log / sqrt grid

Returns:

E_kin – the kinetic energy

Return type:

float

static calc_E_kin_unbound()[source]

Compute the contribution from unbound (continuum) electrons to kinetic energy.

Parameters:

  • orbs (Orbitals) – the KS orbitals object

  • xgrid (ndarray) – the log / sqrt grid

Returns:

E_kin_unbound – the unbound kinetic energy

Return type:

float

Notes

Currently only “ideal” (uniform) approximation for unbound electrons supported.

\[T_\mathrm{ub} = \sum_\sigma \frac{N^\sigma\times V}{\sqrt{2}\pi^2}\ I_{3/2}(\mu,\beta),\]

where \(I_{3/2}(\mu,\beta)\) denotes the complete Fermi-Diract integral

static calc_S_orbs(occnums, degen)[source]

Compute the KS (non-interacting) entropy for specified orbitals (see notes).

Parameters:

  • occnums (ndarray) – orbital occupation numbers

  • degen (ndarray) – product of dos and degeneracy array (\((2l+1)\))

Returns:

S – the bound contribution to the entropy

Return type:

float

Notes

The entropy of non-interacting (KS) electrons is given by:

\[S_\mathrm{b} = -\sum_{s,l,n} g_{ln} (2l+1) [ f_{nls} \log(f_{nls}) \ + (1-f_{nls}) (\log(1-f_{nls}) ]\]
static calc_S_unbound()[source]

Compute the unbound contribution to the entropy.

Parameters:

orbs (Orbitals) – the KS orbitals object

Returns:

S_unbound – the unbound entropy term

Return type:

float

Notes

Currently only “ideal” (uniform) treatment of unbound electrons is supported.

\[S_\mathrm{ub} = \sum_\sigma \frac{V}{\sqrt{2}\pi^2} I_{1/2}(\mu,\beta)\]

where \(I_{1/2}(\mu,\beta)\) is the complete Fermi-Dirac integral of order \(1/2\)

calc_entropy(orbs)[source]

Compute the KS entropy.

Entropy is in general different for bound and unbound orbitals. This function is a wrapper which calls the respective functions.

Parameters:

orbs (Orbitals) – the KS orbitals object

Returns:

S – contains tot, bound and unbound keys

Return type:

dict of floats

property E_en

The electron-nuclear energy.

Type:

float

property E_ha

The Hartree energy.

Type:

float

property E_kin

KS kinetic energy.

Contains bound and unbound keys.

Type:

dict of floats

property E_tot

The total KS internal energy.

Given by \(E=T_\mathrm{s}+E_\mathrm{en}+E_\mathrm{ha}+F_\mathrm{xc}\).

Type:

float

property E_xc

The exchange-correlation energy.

Contains the keys x, c and xc for exchange, correlation and exchange + correlation respectively

Type:

dict of floats

property F_tot

The free energy.

Contains the keys F, E and S for free energy, internal energy and total entropy. \(F = E - TS\)

Type:

dict of floats

property entropy

The total entropy.

Contains bound and unbound keys.

Type:

dict of floats

class atoMEC.staticKS.EnergyAlt(orbs, dens, pot)[source]

Bases: object

Class holding information about the KS energy and associated routines.

N.B. this computes the total energy functional in an alternative way to the main Energy object. In this class, the total energy is first calculated from the sum of the eigenvalues and then \(\int \mathrm{d}r v_\mathrm{Hxc}(r) n(r)\) is subtracted to obtain the sum over kinetic and electron-nuclear energies.

In general, this gives a different result compared to the Energy, due to the way the unbound electrons are treated. Which result is ‘correct’ depends on how we interpret the unbound energy (i.e., is it purely a kinetic term or not).

static calc_E_v_hxc(dens, pot, xgrid, grid_type)[source]

Compute the compensating integral term over the Hxc-potential (see notes).

Parameters:

  • dens (ndarray) – the (spin) KS density

  • pot (Potential) – the KS potential object

  • xgrid (ndarray) – the log / sqrt grid

Returns:

E_v_hxc – the compensating Hxc-potential integral

Return type:

float

Notes

In this construction of the KS energy functional, the sum over the KS eigenvalues must be compensated by an integral of the Hxc-potential. This integral is given by:

\[E = 4\pi\sum_\sigma\int \mathrm{d}r r^2 n^\sigma(r) v_\mathrm{Hxc}^\sigma(r)\]
calc_entropy(orbs)[source]

Compute the KS entropy.

Entropy is in general different for bound and unbound orbitals. This function is a wrapper which calls the respective functions.

Parameters:

orbs (Orbitals) – the KS orbitals object

Returns:

S – contains tot, bound and unbound keys

Return type:

dict of floats

property E_en

Electron-nuclear attraction energy.

Type:

float

property E_eps

The sum of the (weighted) eigenvalues.

Given by \(E=\sum_{nl\sigma} (2l+1) f_{nl}^\sigma \epsilon_{nl}^\sigma\)

Type:

float

property E_ha

The Hartree energy.

Type:

float

property E_kin

Kinetic energy components.

Type:

Dict of floats

property E_tot

The total KS internal energy.

Given by \(E=T_\mathrm{s}+E_\mathrm{en}+E_\mathrm{ha}+F_\mathrm{xc}\).

Type:

float

property E_unbound

The energy of the unbound part of the electron density.

Type:

float

property E_v_hxc

The integral \(\int \mathrm{d}r v_\mathrm{Hxc}(r) n(r)\).

Type:

float

property E_xc

The exchange-correlation energy.

Contains the keys x, c and xc for exchange, correlation and exchange + correlation respectively

Type:

dict of floats

property F_tot

The free energy.

Contains the keys F, E and S for free energy, internal energy and total entropy. \(F = E - TS\)

Type:

dict of floats

property entropy

The total entropy.

Contains bound and unbound keys.

Type:

dict of floats

class atoMEC.staticKS.GramSchmidt(eigfuncs, xgrid, grid_type)[source]

Bases: object

Class holding Gram-Schmidt orthoganalization process.

make_ortho()[source]

Make eigenfunctions orthonormal using Gram-Schmidt orthoganalization.

Parameters:

None –

Returns:

eigfuncs_ortho – orthonormal eigenfunctions

Return type:

ndarray

static prod_eigfuncs(phi0, phi1, xgrid)[source]

Compute the product of two KS eigenfunctions and integrate.

Parameters:

  • phi_0 (ndarray) – first orbital

  • phi_1 (ndarray) – second orbital

  • xgrid (ndarray) – log grid

Returns:

prod_eigfuncs_ – integrated product of phi_0 and phi_1

Return type:

ndarray

static proj_eigfuncs(phi0, phi1, xgrid)[source]

Compute the projection of one eigenfunction onto another.

Parameters:

  • phi_0 (ndarray) – first orbital

  • phi_1 (ndarray) – second orbital

  • xgrid (ndarray) – log grid

Returns:

proj_eigfuncs_ – projection of phi_0 onto phi_1

Return type:

ndarray

class atoMEC.staticKS.Orbitals(xgrid, grid_type)[source]

Bases: object

Class holding the KS orbitals, associated quantities and relevant routines.

static calc_DOS_sum(eigs_min, eigs_max, ldegen)[source]

Compute the summed density-of-states using the method of Massacrier et al.

This function sums the DOS over l and n quantum numbers and is used for writing the output of the DOS only. Refer to make_DOS_bands for details on the functional form of the DOS.

Parameters:

  • eigs_min (ndarray) – the lower bound of the energy band

  • eigs_max (ndarray) – the upper bound of the energy band

  • ldegen (ndarray) – degeneracy matrix

Returns:

  • e_arr (ndarray) – sorted array of energies

  • fd_dist (ndarray) – Fermi-Dirac occupation numbers

  • DOS_sum (ndarray) – Density-of-states multiplied and summed over (l,n) quanutum numbers.

calc_bands(v, eigfuncs_l, eigfuncs_u, solver)[source]

Compute the eigenfunctions which fill the energy bands.

Parameters:

  • v (ndarray) – the KS potential

  • eigfuncs_l (ndarray) – the lower bound (neumann) eigenfunctions

Returns:

  • eigvals (ndarray) – the orbital energies across all the bands

  • eigfuncs (ndarray) – the KS eigenfunctions for all energies covered

static calc_occnums(eigvals, mu)[source]

Compute the Fermi-Dirac occupations for the eigenvalues.

Parameters:

  • eigvals (ndarray) – the KS eigenvalues

  • lbound (ndarray) – the array of bound eigenvalues

  • mu (aray_like) – the chemical potential

Returns:

occnums – the orbital occupations multiplied with (2l+1) degeneracy

Return type:

ndarray

static check_orbs(occnums_w, threshold)[source]

Check the values of nmax and lmax are sufficient.

Finds the values of the occupations of the final orbitals in lmax and nmax directions. If either is above the threshold, returns False (which triggers a warning elsewhere).

Parameters:

  • occnums_w (np.ndarray) – weighted orbital occupations

  • threshold (float) – the threshold occupation number at which to trigger a warning

Returns:

norbs_ok, lorbs_ok – whether the values of nmax and lmax are sufficient

Return type:

tuple of bools

compute(potential, bc, init=False, eig_guess=False)[source]

Compute the orbitals and their eigenvalues with the given potential.

Parameters:

  • potential (ndarray) – the KS (or alternatively chosen) potential

  • init (bool, optional) – whether orbitals are being computed for first time

Return type:

None

static make_DOS_bands(eigs_min, eigs_max, eigvals)[source]

Compute the density-of-states using the method of Massacrier (see notes, [6]).

Parameters:

  • eigs_min (ndarray) – the lower bound of the energy band

  • eigs_max (ndarray) – the upper bound of the energy band

  • eigvals (ndarray) – the KS energy eigenvalues

Returns:

dos – the density-of-states

Return type:

ndarray

Notes

The density-of-states is defined in this model as:

\[g(\epsilon) = \frac{2}{\pi \delta^2} \sqrt{(\epsilon^+ - \epsilon)\ (\epsilon - \epsilon^-)},\ \delta = \frac{1}{2} (\epsilon^+ - \epsilon^-),\]

where \(\epsilon^\pm\) are the upper and lower band limits respectively.

References

static make_e_arr(eigvals_min, eigvals_max, sp)[source]

Make the energy array for the bands boundary condition.

Energies are populated from the lowest band up to the maximum specified energy, with band gaps (forbidden energies) accounted for. The array is non-linear in order to optimize computation time.

Parameters:

  • eigvals_min (ndarray) – the lower bound of the energy bands

  • eigvals_max (ndarray) – the upper bound of the energy bands

Returns:

e_tot_arr – the banded energy array

Return type:

ndarray

static make_ldegen(eigvals)[source]

Construct the lbound matrix denoting the bound states and their degeneracies.

For each spin channel, \(L^\mathrm{B}_{ln}=(2l+1)\times\Theta(\epsilon_{nl})\)

Parameters:

eigvals (ndarray) – the KS orbital eigenvalues

Returns:

lbound_mat – the lbound matrix

Return type:

ndarray

occupy()[source]

Occupy the KS orbitals according to Fermi-Dirac statistics.

Parameters:

None –

Return type:

None

property DOS

Density of states (DOS) matrix.

Type:

ndarray

property eigfuncs

The radial KS eigenfunctions on the log / sqrt grid.

Related to real-grid KS radial orbitals for log grid by: \(X_{nl}^{\sigma}(x)=e^{x/2} R_{nl}^{\sigma}(r).\)

And for the sqrt grid by (no transformation): \(X_{nl}^{\sigma}(s) = R_{nl}^{\sigma}(r).\)

Type:

ndarray

property eigvals

The KS eigenvalues \(\epsilon_{nl}^\sigma\).

Type:

ndarray

property eigvals_max

Upper bound (for bands bc) of KS eigenvalues.

Type:

ndarray

property eigvals_min

Lower bound (for bands bc) of KS eigenvalues.

Type:

ndarray

property kpt_int_weight

The integration weight for the Fermi-Dirac / DOS integral.

Type:

ndarray

property ldegen

Angular momentum degeneracy matrix.

Type:

ndarray

property occ_weight

KS occupation number weightings.

The occupation weighting is the product of the DOS, degeneracy of the angular momentum, and the integration weight.

Type:

ndarray

property occnums

Bare KS occupation numbers (Fermi-Dirac).

Type:

ndarray

property occnums_w

Weighted KS occupation numbers.

Type:

ndarray

class atoMEC.staticKS.Potential(density)[source]

Bases: object

Class holding the KS potential and the routines required to compute it.

static calc_v_en(xgrid, grid_type)[source]

Construct the electron-nuclear potential.

The electron-nuclear potential is given by \(v_\mathrm{en} (x) = -Z * \exp(-x)\) on the log / sqrt grid

static calc_v_ha(density, xgrid, grid_type)[source]

Construct the Hartree potential (see notes).

Parameters:

  • density (ndarray) – the total KS density

  • xgrid (ndarray) – the log / sqrt grid

Notes

\(v_\mathrm{ha}\) is defined on the r-grid as:

\[v_\mathrm{ha}(r) = 4\pi\int_0^r \mathrm{d}r' r'^2 \frac{n(r')}{\max(r,r')}\]

On the log-grid:

\[v_\mathrm{ha}(x) = 4\pi\Big\{\exp(-x)\int_{x0}^x \mathrm{d}x' n(x') \ \exp(3x') -\int_x^{\log(r_s)} \mathrm{d}x' n(x') \exp(2x') \Big\}\]
property v_en

The electron-nuclear potential.

Type:

ndarray

property v_ha

The Hartree potential.

Type:

ndarray

property v_s

The full KS potential.

Given by \(v_\mathrm{s} = v_\mathrm{en} + v_\mathrm{ha} + v_\mathrm{xc}\).

Type:

ndarray

property v_xc

The exchange-correlation potential.

Contains the keys x, c and xc denoting exchange, correlation, and exchange + correlation respectively.

Type:

dict of ndarrays

atoMEC.staticKS.log_grid(x_r)[source]

Set up the logarithmic (and real) grid, defined up to x_r.

The leftmost grid point is given by the config.grid_params[“x0”] parameter.

Parameters:

x_r (float) – The RHS grid point (in logarithmic space)

Returns:

xgrid, rgrid – The grids in logarithmic (x) and real (r) space

Return type:

tuple of ndarrays

atoMEC.staticKS.sqrt_grid(s_r)[source]

Set up the sqrt (and real) grid, defined up to d_r.

The leftmost grid point is given by the config.grid_params[“s0”] parameter.

Parameters:

s_r (float) – The RHS grid point (in sqrt space)

Returns:

xgrid, rgrid – The grids in sqrt (x) and real (r) space

Return type:

tuple of ndarrays

atoMEC.unitconv module

Contains unit conversions.

atoMEC.writeoutput module

Contains the routines which generate printed output (to screen and separate files).

Classes

  • SCF : Write information about the self-consistent field (SCF) cycle.

Functions

class atoMEC.writeoutput.SCF[source]

Bases: object

Write information about the self-consistent field (SCF) cycle.

static write_cycle(iscf, E_free, conv_vals)[source]

Write output string for each SCF iteration.

Parameters:

  • iscf (int) – the iteration number

  • E_free (float) – the total free energy

  • conv_vals (dict) – dictionary of convergence values

Returns:

output_str – free energy and convergence values for the i-th iteration

Return type:

str

write_final(energy, orbitals, density, conv_vals)[source]

Write final post-SCF information about energy, orbitals and convergence.

Parameters:
Returns:

output_str – information about whether the calculation converged, the final energies (broken down into components), and the orbital energies and occupations

Return type:

str

static write_final_energies(energy)[source]

Write formatted KS energy information (by component).

Parameters:

energy (staticKS.Energy) – the KS energy object

Returns:

output_str – formatted output string of energies by component

Return type:

str

static write_init()[source]

Write the initial spiel for an SCF calculation.

Parameters:

None –

Returns:

output_str – header string for SCF cycle

Return type:

str

static write_orb_info(orbitals)[source]

Write formatted KS orbital information.

Information up to to the highest occupied level +1 (for both n and l) is included; remaining levels are truncated.

Parameters:

orbitals (staticKS.Orbitals) – the KS orbitals object

Returns:

eigs_occs – a tuple containing formatted eigenvalue table eigval_tbl, and formatted occupation numbers table occnum_tbl

Return type:

tuple[str,str]

atoMEC.writeoutput.density_to_csv(rgrid, density, filename)[source]

Write the density (on the r-grid) to file.

Parameters:

  • rgrid (ndarray) – real-space grid

  • density (staticKS.Density) – the KS density object

  • filename (str) – name of the file to write to

atoMEC.writeoutput.dos_to_csv(orbitals, filename)[source]

Write the energy eigenvalues, Fermi-Dirac occupations and DOS to file.

Parameters:

  • orbitals (staticKS.Orbitals) – the orbitals object

  • filename (str) – name of the file to write to

Return type:

None

atoMEC.writeoutput.eigs_occs_to_csv(orbitals, filename)[source]

Write all the orbital energies and their occupations to file.

Parameters:

  • orbitals (staticKS.Orbitals) – the orbitals object

  • filename (str) – name of the file to write to

Return type:

None

atoMEC.writeoutput.potential_to_csv(rgrid, potential, filename)[source]

Write the potential (on the r-grid) to file.

Parameters:

  • rgrid (ndarray) – real-space grid

  • density (staticKS.Potential) – the KS potential object

  • filename (str) – name of the file to write to

atoMEC.writeoutput.timing(f)[source]

Generate timing information for given input function.

Parameters:

f (func) – the function to be timed

atoMEC.writeoutput.write_ISModel_data(ISModel)[source]

Write information about the approximations used for the IS model.

Parameters:

ISModel (models.ISModel) – The ISModel object

Returns:

output_str – formatted description of the ISModel attributes

Return type:

str

atoMEC.writeoutput.write_atomic_data(atom)[source]

Write information about the main Atom object.

Parameters:

atom (atoMEC.Atom) – the main Atom object

Returns:

output_str – formatted description of the Atom attributes

Return type:

str

atoMEC.xc module

The xc module calculates exchange-correlation energies and potentials.

Mostly it sets up inputs and makes appropriate calls to the libxc package. It also includes a class for customized xc-functionals and checks the validity of inputs.

Classes

  • XCFunc : handles atoMEC-defined functionals (not part of libxc package)

Functions

  • check_xc_func() : check the xc func input is recognised and valid for atoMEC

  • set_xc_func() : initialize the xc functional object

  • v_xc() : return the xc potential on the grid

  • E_xc() : return the xc energy

  • calc_xc() : compute the xc energy or potential depending on arguments

class atoMEC.xc.XCFunc(xc_code)[source]

Bases: object

Defines the XCFunc object for ‘special’ (non-libxc) funcs.

This is desgined to align with some proprerties of libxc func objects in order to make certain general calls more straightforward

Parameters:

xc_code (str) – the string identifier for the x/c func

Notes

The XCFunc object contains no public attributes, but it does contain the following private attributes (named so to match libxc):

_xc_codestr

the string identifier for the x/c functional (not a libxc attribute)

_xc_func_namestr

name of the functional

_numberint

functional id

_familystr

the family to which the xc functional belongs

atoMEC.xc.E_xc(density, xgrid, xfunc, cfunc, grid_type)[source]

Retrieve the xc energy.

Parameters:

  • density (ndarray) – the KS density on the log/sqrt grid

  • xgrid (ndarray) – the log/sqrt grid

  • xfunc (XCFunc or pylibxc.LibXCFunctional) – the exchange functional object

  • cfunc (XCFunc or pylibxc.LibXCFunctional) – the correlation functional object

Returns:

_E_xc – dictionary containing terms x, c and xc for exchange, correlation and exchange + correlation respectively

Return type:

dict of floats

atoMEC.xc.calc_xc(density, xgrid, xcfunc, xctype)[source]

Compute the x/c energy density or potential depending on xctype input.

Parameters:

  • density (ndarray) – the KS density on the log/sqrt grid

  • xgrid (ndarray) – the log/sqrt grid

  • xcfunc (XCFunc or pylibxc.LibXCFunctional) – the exchange or correlation functional object

  • xctype (str) – the quantity to return, e_xc for energy density or v_xc for potential

Returns:

xc_arr – xc energy density or potential depending on xctype

Return type:

ndarray

atoMEC.xc.check_xc_func(xc_code, id_supp)[source]

Check the xc input string or id is valid.

Parameters:

  • xc_code (str or int) – the name or libxc id of the xc functional

  • id_supp (list of ints) – supported families of xc functional

Returns:

xc_func_name – the xc functional name

Return type:

str or None

atoMEC.xc.gga_pot_chainrule(libxc_output, grad_0, grad_1, xgrid, spindims)[source]

Calculate a component for the spin-polarized gga xc potential.

Parameters:

  • libxc_output (dict of libxc objects) – The output of the libxc calculation

  • grad_0 (ndarray) – The gradient of the density of the 0 spin channel.

  • grad_1 (ndarray) – The gradient of the density of the 1 spin channel.

  • xgrid (ndarray) – the log / sqrt grid.

  • spindim (int) – the number of spin dimentions in the calculation.

Returns:

gga_addition2 – The addition to obtain the spin polarized gga xc potential

Return type:

nd array

Notes

Calculates the second component of the chain rule expansion for the gga xc pot calculation for a spin polarized calculation:

\[v_{xc}^{\sigma}(r)=\frac{\partial e_{xc}[\rho_0,\rho_1]} {\partial\rho_{\sigma}}- 2\nabla[\nabla \rho_{\sigma}(r)v_{2}^{\sigma,\sigma}(r)+ \nabla \rho_{\sigma'}(r)v_{2}^{\sigma,\sigma'}(r)]\]

where

\[v_{2}^{\sigma,\sigma'}= \frac{\partial e_{xc}[\rho_0,\rho_1]} {\partial (\nabla \rho_{\sigma} \nabla \rho_{\sigma'})}\]

The output of the function is the square brackets after being acted upon with the divergence.

atoMEC.xc.set_xc_func(xc_code)[source]

Initialize the xc functional object.

Parameters:

xc_code (str or int) – the name or id of the libxc functional

Returns:

xc_func

Return type:

XCFunc or pylibxc.LibXCFunctional

atoMEC.xc.v_xc(density, xgrid, xfunc, cfunc, grid_type)[source]

Retrive the xc potential.

Parameters:

  • density (ndarray) – the KS density on the log/sqrt grid

  • xgrid (ndarray) – the log/sqrt grid

  • xfunc (XCFunc or pylibxc.LibXCFunctional) – the exchange functional object

  • cfunc (XCFunc or pylibxc.LibXCFunctional) – the correlation functional object

Returns:

_v_xc – dictionary containing terms x, c and xc for exchange, correlation and exchange + correlation respectively

Return type:

dict of ndarrays

Module contents

atoMEC: Average-atom code for matter under extreme conditions.

Copyright (c) 2021 (in alphabetical order), Tim Callow, Attila Cangi, Eli Kraisler. All rights reserved.

atoMEC is a python-based average-atom code for simulations of high energy density phenomena such as in warm dense matter. Please see the README or the project wiki (https://atomec-project.github.io/atoMEC/) for more information.

Classes

  • Atom : the main object for atoMEC calculations, containing information about physical material properties

class atoMEC.Atom(species, temp, radius=-1, density=-1, charge=0, units_temp='ha', units_radius='bohr', units_density='g/cm3', write_info=True)[source]

Bases: object

The principal object in atoMEC calculations which defines the material properties.

The Atom contains key information about the physical properties of the material such as temperature, density, and charge. It does not contain any information about approximations or choices of model.

Parameters:

  • species (str) – The chemical symbol for the atomic species, eg “He”

  • temp (float) – The electronic temperature in hartree, eV or Kelvin

  • radius (float, optional) – The radius of the Wigner-Seitz sphere, defined as 0.5*a_i, where a_i is the average inter-atomic distance

  • density (float, optional) – The mass density of the material in \(\mathrm{g\ cm}^{-3}\)

  • charge (int, optional) – The overall net charge

  • units_temp (str, optional) – The units of temperature, must be one of “ha”, “ev” or “k”

  • units_radius (str, optional) – The units of radius, must be one of “ang” or “bohr”

  • units_density (str, optional) – The units of density, currently only “g/cm3” is supported

  • write_output (bool, optional) – Whether to print atomic information, defaults True

property E_Fermi

the Fermi energy.

Type:

float

property WS_radius

the Wigner-Seitz radius, or the electron coupling parameter.

The Wigner-Seitz radius differs from the Voronoi radius because it depends on the free electron density, defined in atoMEC as the valence electron density.

Type:

float

property at_chrg

the atomic charge Z.

Type:

int

property at_mass

the atomic mass (in a.u.).

Type:

float

property charge

the net charge of the atom.

Type:

int

property density

the mass density of the material in \(\mathrm{g\ cm}^{-3}\).

Type:

float

property gamma_ion

the ionic coupling parameter.

Type:

float

property info

formatted information about the Atom’s attributes.

Type:

str

property nele

the number of electrons in the atom.

The total electron number is given by the sum of at_chrg and charge.

Type:

int

property nvalence

the number of valence electrons.

Type:

int

property radius

radius of the Voronoi sphere.

The radius is defined as \(a_i /2\), where \(a_i\) is the average inter-atomic distance.

Type:

float

property species

the chemical symbol for the atomic species.

Type:

str

property temp

the electronic temperature in Hartree.

Type:

float

property theta_e

the electron degeneracy parameter.

Type:

float

property units_density

the units of the atomic density.

Type:

str

property units_radius

the units of the atomic radius.

Type:

str

property units_temp

the units of temperature.

Type:

str