# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at https://mozilla.org/MPL/2.0/.
from __future__ import annotations
from typing import Literal, Optional
from sisl.messages import deprecate_argument
r"""Electron related functions and classes
==========================================
In sisl electronic structure calculations are relying on routines
specific for electrons. For instance density of states calculations from
electronic eigenvalues and other quantities.
This module implements the necessary tools required for calculating
DOS, PDOS, band-velocities and spin moments of non-colinear calculations.
One may also plot real-space wavefunctions.
DOS
PDOS
COP
berry_phase
berry_curvature
conductivity
wavefunction
spin_moment
spin_contamination
Supporting classes
------------------
Certain classes aid in the usage of the above methods by implementing them
using automatic arguments.
For instance, the PDOS method requires the overlap matrix in non-orthogonal
basis sets at the :math:`\mathbf k`-point corresponding to the eigenstates. Hence, the
argument ``S`` must be :math:`\mathbf S(\mathbf k)`. The `EigenstateElectron` class
automatically passes the correct ``S`` because it knows the states :math:`\mathbf k`-point.
CoefficientElectron
StateElectron
StateCElectron
EigenvalueElectron
EigenvectorElectron
EigenstateElectron
"""
from functools import reduce
import numpy as np
from numpy import (
add,
ceil,
conj,
cos,
dot,
einsum,
empty,
exp,
floor,
int32,
log,
ogrid,
pi,
sin,
sort,
zeros,
)
from scipy.sparse import csr_matrix, hstack, identity, issparse
import sisl._array as _a
from sisl import BoundaryCondition as BC
from sisl import Geometry, Grid, Lattice, constant, units
from sisl._core.oplist import oplist
from sisl._help import dtype_complex_to_real, dtype_real_to_complex
from sisl._indices import indices_le
from sisl._internal import set_module
from sisl._math_small import xyz_to_spherical_cos_phi
from sisl.linalg import det
from sisl.linalg import eigvals as la_eigvals
from sisl.linalg import sqrth, svd_destroy
from sisl.messages import SislError, info, progressbar, warn
from .distribution import get_distribution
from .sparse import SparseOrbitalBZSpin
from .spin import Spin
from .state import Coefficient, State, StateC, _FakeMatrix, degenerate_decouple
__all__ = ["DOS", "PDOS", "COP"]
__all__ += ["spin_moment", "spin_contamination"]
__all__ += ["berry_phase", "berry_curvature"]
__all__ += ["conductivity"]
__all__ += ["wavefunction"]
__all__ += ["CoefficientElectron", "StateElectron", "StateCElectron"]
__all__ += ["EigenvalueElectron", "EigenvectorElectron", "EigenstateElectron"]
[docs]
@set_module("sisl.physics.electron")
def DOS(E, eig, distribution="gaussian"):
r"""Calculate the density of states (DOS) for a set of energies, `E`, with a distribution function
The :math:`\mathrm{DOS}(E)` is calculated as:
.. math::
\mathrm{DOS}(E) = \sum_i D(E-\epsilon_i) \approx\delta(E-\epsilon_i)
where :math:`D(\Delta E)` is the distribution function used. Note that the distribution function
used may be a user-defined function. Alternatively a distribution function may
be retrieved from `~sisl.physics.distribution`.
Parameters
----------
E : array_like
energies to calculate the DOS at
eig : array_like
electronic eigenvalues
distribution : func or str, optional
a function that accepts :math:`\Delta E` as argument and calculates the
distribution function.
See Also
--------
sisl.physics.distribution : a selected set of implemented distribution functions
COP : calculate COOP or COHP curves
PDOS : projected DOS (same as this, but projected onto each orbital)
spin_moment : spin moment
Returns
-------
numpy.ndarray
DOS calculated at energies, has same length as `E`
"""
if isinstance(distribution, str):
distribution = get_distribution(distribution)
return reduce(lambda DOS, eig: DOS + distribution(E - eig), eig, 0.0)
[docs]
@set_module("sisl.physics.electron")
def PDOS(E, eig, state, S=None, distribution="gaussian", spin=None):
r""" Calculate the projected density of states (PDOS) for a set of energies, `E`, with a distribution function
The :math:`\mathrm{PDOS}(E)` is calculated as:
.. math::
\mathrm{PDOS}_i(E) = \sum_\alpha \psi^*_{\alpha,i} [\mathbf S | \psi_{\alpha}\rangle]_i D(E-\epsilon_\alpha)
where :math:`D(\Delta E)` is the distribution function used. Note that the distribution function
used may be a user-defined function. Alternatively a distribution function may
be aquired from `~sisl.physics.distribution`.
In case of an orthogonal basis set :math:`\mathbf S` is equal to the identity matrix.
Note that `DOS` is the sum of the orbital projected DOS:
.. math::
\mathrm{DOS}(E) = \sum_i\mathrm{PDOS}_i(E)
For non-colinear calculations (this includes spin-orbit calculations) the PDOS is additionally
separated into 4 components (in this order):
- Total projected DOS
- Projected spin magnetic moment along :math:`x` direction
- Projected spin magnetic moment along :math:`y` direction
- Projected spin magnetic moment along :math:`z` direction
These are calculated using the Pauli matrices :math:`\boldsymbol\sigma_x`, :math:`\boldsymbol\sigma_y` and :math:`\boldsymbol\sigma_z`:
.. math::
\mathrm{PDOS}_i^\sigma(E) &= \sum_\alpha \psi^*_{\alpha,i} \boldsymbol\sigma_z \boldsymbol\sigma_z [\mathbf S | \psi_\alpha\rangle]_i D(E-\epsilon_\alpha)
\\
\mathrm{PDOS}_i^x(E) &= \sum_\alpha \psi^*_{\alpha,i} \boldsymbol\sigma_x [\mathbf S | \psi_\alpha\rangle]_i D(E-\epsilon_\alpha)
\\
\mathrm{PDOS}_i^y(E) &= \sum_\alpha \psi^*_{\alpha,i} \boldsymbol\sigma_y [\mathbf S | \psi_\alpha\rangle]_i D(E-\epsilon_\alpha)
\\
\mathrm{PDOS}_i^z(E) &= \sum_\alpha\psi^*_{\alpha,i} \boldsymbol\sigma_z [\mathbf S | \psi_\alpha\rangle]_i D(E-\epsilon_\alpha)
Note that the total PDOS may be calculated using :math:`\boldsymbol\sigma_\gamma\boldsymbol\sigma_\gamma` where :math:`\gamma` may be either of :math:`x`,
:math:`y` or :math:`z`.
Parameters
----------
E : array_like
energies to calculate the projected-DOS from
eig : array_like
eigenvalues
state : array_like
eigenvectors
S : array_like, optional
overlap matrix used in the :math:`\langle\psi|\mathbf S|\psi\rangle` calculation. If `None` the identity
matrix is assumed. For non-colinear calculations this matrix may be halve the size of ``len(state[0, :])`` to
trigger the non-colinear calculation of PDOS.
distribution : func or str, optional
a function that accepts :math:`E-\epsilon` as argument and calculates the
distribution function.
spin : str or Spin, optional
the spin configuration. This is generally only needed when the eigenvectors correspond to a non-colinear
calculation.
See Also
--------
sisl.physics.distribution : a selected set of implemented distribution functions
DOS : total DOS (same as summing over orbitals)
COP : calculate COOP or COHP curves
spin_moment : spin moment
Returns
-------
numpy.ndarray
projected DOS calculated at energies, has dimension ``(1, state.shape[1], len(E))``.
For non-colinear calculations it will be ``(4, state.shape[1] // 2, len(E))``, ordered as
indicated in the above list.
"""
if isinstance(distribution, str):
distribution = get_distribution(distribution)
# Figure out whether we are dealing with a non-colinear calculation
if S is None:
class S:
__slots__ = []
shape = (state.shape[1], state.shape[1])
@staticmethod
def dot(v):
return v
if spin is None:
if S.shape[1] == state.shape[1] // 2:
spin = Spin("nc")
S = S[::2, ::2]
else:
spin = Spin()
# check for non-colinear (or SO)
if spin.kind > Spin.POLARIZED:
# Non colinear eigenvectors
if S.shape[1] == state.shape[1]:
# Since we are going to reshape the eigen-vectors
# to more easily get the mixed states, we can reduce the overlap matrix
S = S[::2, ::2]
# Initialize data
PDOS = empty(
[4, state.shape[1] // 2, len(E)], dtype=dtype_complex_to_real(state.dtype)
)
# Do spin-box calculations:
# PDOS[0] = total DOS (diagonal)
# PDOS[1] = x == < psi | \sigma_x S | psi >
# PDOS[2] = y == < psi | \sigma_y S | psi >
# PDOS[3] = z == < psi | \sigma_z S | psi >
d = distribution(E - eig[0]).reshape(1, -1)
cs = conj(state[0]).reshape(-1, 2)
v = S.dot(state[0].reshape(-1, 2))
D1 = (cs * v).real # uu,dd PDOS
PDOS[0, :, :] = D1.sum(1).reshape(-1, 1) * d # total DOS
PDOS[3, :, :] = (D1[:, 0] - D1[:, 1]).reshape(-1, 1) * d # z-dos
D1 = (cs[:, 1] * v[:, 0]).reshape(-1, 1) # d,u
D2 = (cs[:, 0] * v[:, 1]).reshape(-1, 1) # u,d
PDOS[1, :, :] = (D1.real + D2.real) * d # x-dos
PDOS[2, :, :] = (D2.imag - D1.imag) * d # y-dos
for i in range(1, len(eig)):
d = distribution(E - eig[i]).reshape(1, -1)
cs = conj(state[i]).reshape(-1, 2)
v = S.dot(state[i].reshape(-1, 2))
D1 = (cs * v).real
PDOS[0, :, :] += D1.sum(1).reshape(-1, 1) * d
PDOS[3, :, :] += (D1[:, 0] - D1[:, 1]).reshape(-1, 1) * d
D1 = (cs[:, 1] * v[:, 0]).reshape(-1, 1)
D2 = (cs[:, 0] * v[:, 1]).reshape(-1, 1)
PDOS[1, :, :] += (D1.real + D2.real) * d
PDOS[2, :, :] += (D2.imag - D1.imag) * d
else:
PDOS = (conj(state[0]) * S.dot(state[0])).real.reshape(-1, 1) * distribution(
E - eig[0]
).reshape(1, -1)
for i in range(1, len(eig)):
PDOS += (conj(state[i]) * S.dot(state[i])).real.reshape(
-1, 1
) * distribution(E - eig[i]).reshape(1, -1)
PDOS.shape = (1, *PDOS.shape)
return PDOS
[docs]
@set_module("sisl.physics.electron")
def COP(E, eig, state, M, distribution="gaussian", tol=1e-10):
r"""Calculate the Crystal Orbital Population for a set of energies, `E`, with a distribution function
The :math:`\mathrm{COP}(E)` is calculated as:
.. math::
\mathrm{COP}_{i,j}(E) = \sum_\alpha \psi^*_{\alpha,i}\psi_{\alpha,j} \mathbf M e^{i\mathbf k\cdot \mathbf R} D(E-\epsilon_\alpha)
where :math:`D(\Delta E)` is the distribution function used. Note that the distribution function
used may be a user-defined function. Alternatively a distribution function may
be aquired from `~sisl.physics.distribution`.
The COP curves generally refers to COOP or COHP curves.
COOP is the Crystal Orbital Overlap Population with `M` being the overlap matrix.
COHP is the Crystal Orbital Hamiltonian Population with `M` being the Hamiltonian.
Parameters
----------
E : array_like
energies to calculate the COP from
eig : array_like
eigenvalues
state : array_like
eigenvectors
M : array_like
matrix used in the COP curve.
distribution : func or str, optional
a function that accepts :math:`E-\epsilon` as argument and calculates the
distribution function.
tol : float, optional
tolerance for considering an eigenstate to contribute to an energy point,
a higher value means that more energy points are discarded and so the calculation
is faster.
Notes
-----
This is not tested for non-collinear states.
This requires substantial amounts of memory for big systems with lots of energy points.
This method is considered experimental and implementation may change in the future.
See Also
--------
sisl.physics.distribution : a selected set of implemented distribution functions
DOS : total DOS
PDOS : projected DOS over all orbitals
spin_moment : spin moment
Returns
-------
~sisl.oplist.oplist
COP calculated at energies, has dimension ``(len(E), *M.shape)``.
"""
if isinstance(distribution, str):
distribution = get_distribution(distribution)
assert len(eig) == len(
state
), "COP: number of eigenvalues and states are not consistent"
# get default dtype
dtype = dtype_complex_to_real(state.dtype)
# initialize the COP values
no = M.shape[0]
n_s = M.shape[1] // M.shape[0]
if isinstance(M, _FakeMatrix):
# A fake matrix equals the identity matrix.
# Hence we can do all calculations only on the diagonal,
# then finally we recreate the full matrix dimensions.
cop = oplist(zeros(no, dtype=dtype) for _ in range(len(E)))
def new_list(bools, tmp, we):
for bl, w in zip(bools, we):
if bl:
yield tmp * w
else:
yield 0.0
for e, s in zip(eig, state):
# calculate contribution from this state
we = distribution(E - e)
bools = we >= tol
if np.any(bools):
tmp = (s.conj() * s).real
cop += new_list(bools, tmp, we)
# Now recreate the full size (in sparse form)
idx = np.arange(no)
def tosize(diag, idx):
return csr_matrix((diag, (idx, idx)), shape=M.shape)
cop = oplist(tosize(d, idx) for d in cop)
elif issparse(M):
# create the new list
cop0 = M.multiply(0.0).real
cop = oplist(cop0.copy() for _ in range(len(E)))
del cop0
# split M, then we will rejoin afterwards
Ms = [M[:, i * no : (i + 1) * no] for i in range(n_s)]
def new_list(bools, tmp, we):
for bl, w in zip(bools, we):
if bl:
yield tmp.multiply(w)
else:
yield 0.0
for e, s in zip(eig, state):
# calculate contribution from this state
we = distribution(E - e)
bools = we >= tol
if np.any(bools):
s = np.outer(s.conj(), s)
tmp = hstack([m.multiply(s).real for m in Ms])
cop += new_list(bools, tmp, we)
else:
old_shape = M.shape
Ms = M.reshape(old_shape[0], n_s, old_shape[0])
cop = oplist(zeros(old_shape, dtype=dtype) for _ in range(len(E)))
for e, s in zip(eig, state):
we = distribution(E - e)
# expand the state and do multiplication
s = np.outer(s.conj(), s)[:, None, :]
cop += we.reshape(-1, 1, 1) * (Ms * s).real.reshape(old_shape)
return cop
[docs]
@set_module("sisl.physics.electron")
def spin_moment(state, S=None, project=False):
r""" Spin magnetic moment (spin texture) and optionally orbitally resolved moments
This calculation only makes sense for non-colinear calculations.
The returned quantities are given in this order:
- Spin magnetic moment along :math:`x` direction
- Spin magnetic moment along :math:`y` direction
- Spin magnetic moment along :math:`z` direction
These are calculated using the Pauli matrices :math:`\boldsymbol\sigma_x`, :math:`\boldsymbol\sigma_y` and :math:`\boldsymbol\sigma_z`:
.. math::
\mathbf{S}_\alpha^x &= \langle \psi_\alpha | \boldsymbol\sigma_x \mathbf S | \psi_\alpha \rangle
\\
\mathbf{S}_\alpha^y &= \langle \psi_\alpha | \boldsymbol\sigma_y \mathbf S | \psi_\alpha \rangle
\\
\mathbf{S}_\alpha^z &= \langle \psi_\alpha | \boldsymbol\sigma_z \mathbf S | \psi_\alpha \rangle
If `project` is true, the above will be the orbitally resolved quantities.
Parameters
----------
state : array_like
vectors describing the electronic states, 2nd dimension contains the states
S : array_like, optional
overlap matrix used in the :math:`\langle\psi|\mathbf S|\psi\rangle` calculation. If `None` the identity
matrix is assumed. The overlap matrix should correspond to the system and :math:`\mathbf k` point the eigenvectors
has been evaluated at.
project: bool, optional
whether the spin-moments will be orbitally resolved or not
Notes
-----
This routine cannot check whether the input eigenvectors originate from a non-colinear calculation.
If a non-polarized eigenvector is passed to this routine, the output will have no physical meaning.
See Also
--------
DOS : total DOS
PDOS : projected DOS
COP : calculate COOP or COHP curves
Returns
-------
numpy.ndarray
spin moments per state with final dimension ``(3, state.shape[0])``, or ``(3, state.shape[0], state.shape[1]//2)`` if project is true
"""
if state.ndim == 1:
return spin_moment(state.reshape(1, -1), S, project)[0]
if S is None:
class S:
__slots__ = []
shape = (state.shape[1] // 2, state.shape[1] // 2)
@staticmethod
def dot(v):
return v
if S.shape[1] == state.shape[1]:
S = S[::2, ::2]
# see PDOS for details related to the spin-box calculations
if project:
s = empty(
[3, state.shape[0], state.shape[1] // 2],
dtype=dtype_complex_to_real(state.dtype),
)
for i in range(len(state)):
cs = conj(state[i]).reshape(-1, 2)
Sstate = S.dot(state[i].reshape(-1, 2))
D1 = (cs * Sstate).real
s[2, i] = D1[:, 0] - D1[:, 1]
D1 = cs[:, 1] * Sstate[:, 0]
D2 = cs[:, 0] * Sstate[:, 1]
s[0, i] = D1.real + D2.real
s[1, i] = D2.imag - D1.imag
else:
s = empty([3, state.shape[0]], dtype=dtype_complex_to_real(state.dtype))
# TODO consider doing this all in a few lines
# TODO Since there are no energy dependencies here we can actually do all
# TODO dot products in one go and then use b-casting rules. Should be much faster
# TODO but also way more memory demanding!
for i in range(len(state)):
cs = conj(state[i]).reshape(-1, 2)
Sstate = S.dot(state[i].reshape(-1, 2))
D = cs.T @ Sstate
s[2, i] = D[0, 0].real - D[1, 1].real
s[0, i] = D[1, 0].real + D[0, 1].real
s[1, i] = D[0, 1].imag - D[1, 0].imag
return s
[docs]
@set_module("sisl.physics.electron")
def spin_contamination(state_alpha, state_beta, S=None, sum: bool = True):
r""" Calculate the spin contamination value between two spin states
This calculation only makes sense for spin-polarized calculations.
The contamination value is calculated using the following formula:
.. math::
S^2_{\alpha,i} &= \sum_j |\langle \psi_j^\beta | \mathbf S | \psi_i^\alpha \rangle|^2
\\
S^2_{\beta,j} &= \sum_i |\langle \psi_i^\alpha | \mathbf S | \psi_j^\beta \rangle|^2
where :math:`\alpha` and :math:`\beta` are different spin-components.
The arrays :math:`S^2_\alpha` and :math:`S^2_\beta` are returned.
Parameters
----------
state_alpha : array_like
vectors describing the electronic states of spin-channel :math:`\alpha`, 2nd dimension contains the states
state_beta : array_like
vectors describing the electronic states of spin-channel :math:`\beta`, 2nd dimension contains the states
S : array_like, optional
overlap matrix used in the :math:`\langle\psi|\mathbf S|\psi\rangle` calculation. If `None` the identity
matrix is assumed. The overlap matrix should correspond to the system and :math:`\mathbf k` point the eigenvectors
have been evaluated at.
sum:
whether the spin-contamination should be summed for all states (a single number returned).
If false, a spin-contamination per state per spin-channel will be returned.
Notes
-----
`state_alpha` and `state_beta` need not have the same number of states.
Returns
-------
~sisl._core.oplist
spin squared expectation value per spin channel :math:`\alpha` and :math:`\beta`.
If `sum` is true, only a single number is returned (not an `~sisl._core.oplist`, otherwise a list for each
state.
"""
if state_alpha.ndim == 1:
if state_beta.ndim == 1:
ret = spin_contamination(
state_alpha.reshape(1, -1), state_beta.reshape(1, -1), S, sum
)
if sum:
return ret
return oplist((ret[0][0], ret[1][0]))
return spin_contamination(state_alpha.reshape(1, -1), state_beta, S, sum)
elif state_beta.ndim == 1:
return spin_contamination(state_alpha, state_beta.reshape(1, -1), S, sum)
if state_alpha.shape[1] != state_beta.shape[1]:
raise ValueError(
"spin_contamination requires alpha and beta states to have same number of orbitals"
)
n_alpha = state_alpha.shape[0]
n_beta = state_beta.shape[0]
if S is None:
S_state_beta = state_beta.T
else:
S_state_beta = S.dot(state_beta.T)
if sum:
Ss = 0
for s in state_alpha:
D = conj(s) @ S_state_beta
Ss += (D @ D.conj()).real
return Ss
else:
Sa = empty([n_alpha], dtype=dtype_complex_to_real(state_alpha.dtype))
Sb = zeros([n_beta], dtype=Sa.dtype)
# Loop alpha...
for i, s in enumerate(state_alpha):
D = conj(s) @ S_state_beta
D = (D * conj(D)).real
Sb += D
Sa[i] = D.sum()
return oplist((Sa, Sb))
# dHk is in [Ang eV]
# velocity units in [Ang/ps]
_velocity_const = 1 / constant.hbar("eV ps")
def _velocity_matrix_non_ortho(
state, dHk, energy, dSk, degenerate, degenerate_dir, dtype
):
r"""For states in a non-orthogonal basis"""
# All matrix elements along the 3 directions
n = state.shape[0]
v = empty([3, n, n], dtype=dtype)
# Decouple the degenerate states
if not degenerate is None:
degenerate_dir = _a.asarrayd(degenerate_dir)
degenerate_dir /= (degenerate_dir**2).sum() ** 0.5
deg_dHk = sum(d * dh for d, dh in zip(degenerate_dir, dHk))
for deg in degenerate:
# Set the average energy
e = np.average(energy[deg])
energy[deg] = e
# Now diagonalize to find the contributions from individual states
# then re-construct the seperated degenerate states
# Since we do this for all directions we should decouple them all
state[deg] = degenerate_decouple(
state[deg],
deg_dHk - sum(d * e * ds for d, ds in zip(degenerate_dir, dSk)),
)
del deg_dHk
# Since they depend on the state energies and dSk we have to loop them individually.
cs = conj(state)
for s, e in enumerate(energy):
# Since dHk *may* be a csr_matrix or sparse, we have to do it like
# this. A sparse matrix cannot be re-shaped with an extra dimension.
v[0, s] = cs @ (dHk[0] - e * dSk[0]).dot(state[s])
v[1, s] = cs @ (dHk[1] - e * dSk[1]).dot(state[s])
v[2, s] = cs @ (dHk[2] - e * dSk[2]).dot(state[s])
v *= _velocity_const
return v
def _velocity_matrix_ortho(state, dHk, degenerate, degenerate_dir, dtype):
r"""For states in an orthogonal basis"""
# All matrix elements along the 3 directions
n = state.shape[0]
v = empty([3, n, n], dtype=dtype)
# Decouple the degenerate states
if not degenerate is None:
degenerate_dir = _a.asarrayd(degenerate_dir)
degenerate_dir /= (degenerate_dir**2).sum() ** 0.5
deg_dHk = sum(d * dh for d, dh in zip(degenerate_dir, dHk))
for deg in degenerate:
# Now diagonalize to find the contributions from individual states
# then re-construct the seperated degenerate states
# Since we do this for all directions we should decouple them all
state[deg] = degenerate_decouple(state[deg], deg_dHk)
del deg_dHk
cs = conj(state)
for s in range(n):
v[0, s] = cs @ dHk[0].dot(state[s])
v[1, s] = cs @ dHk[1].dot(state[s])
v[2, s] = cs @ dHk[2].dot(state[s])
v *= _velocity_const
return v
[docs]
@set_module("sisl.physics.electron")
def berry_curvature(
state, energy, dHk, dSk=None, degenerate=None, degenerate_dir=(1, 1, 1)
):
r"""Calculate the Berry curvature matrix for a set of states (using Kubo)
The Berry curvature is calculated using the following expression
(:math:`\alpha`, :math:`\beta` corresponding to Cartesian directions):
.. math::
\boldsymbol\Omega_{\alpha\beta,i} = - \frac2{\hbar^2}\Im\sum_{j\neq i}
\frac{v^\alpha_{ij} v^\beta_{ji}}
{[\epsilon_j - \epsilon_i]^2}
For details see Eq. (11) in :cite:`Wang2006` or Eq. (2.59) in :cite:`TopInvCourse`.
Parameters
----------
state : array_like
vectors describing the electronic states, 2nd dimension contains the states. In case of degenerate
states the vectors *may* be rotated upon return.
energy : array_like, optional
energies of the states. In case of degenerate
states the eigenvalues of the states will be averaged in the degenerate sub-space.
dHk : list of array_like
Hamiltonian derivative with respect to :math:`\mathbf k`. This needs to be a tuple or
list of the Hamiltonian derivative along the 3 Cartesian directions.
dSk : list of array_like, optional
:math:`\delta \mathbf S_{\mathbf k}` matrix required for non-orthogonal basis.
Same derivative as `dHk`.
NOTE: Using non-orthogonal basis sets are not tested.
degenerate : list of array_like, optional
a list containing the indices of degenerate states. In that case a prior diagonalization
is required to decouple them.
degenerate_dir : (3,), optional
along which direction degenerate states are decoupled.
See Also
--------
velocity : calculate state velocities
velocity_matrix : calculate state velocities between all states
Hamiltonian.dHk : function for generating the Hamiltonian derivatives (`dHk` argument)
Hamiltonian.dSk : function for generating the Hamiltonian derivatives (`dSk` argument)
Returns
-------
numpy.ndarray
Berry flux with final dimension ``(3, 3, state.shape[0])``
"""
if state.ndim == 1:
return berry_curvature(
state.reshape(1, -1), energy, dHk, dSk, degenerate, degenerate_dir
)[0]
# cast dtypes to *any* complex valued data-type that can be expressed
# minimally by a complex64 object
dtype = np.result_type(state.dtype, dHk[0].dtype, np.complex64)
if dSk is None:
v_matrix = _velocity_matrix_ortho(state, dHk, degenerate, degenerate_dir, dtype)
else:
v_matrix = _velocity_matrix_non_ortho(
state, dHk, energy, dSk, degenerate, degenerate_dir, dtype
)
warn(
"berry_curvature calculation for non-orthogonal basis sets are not tested! Do not expect this to be correct!"
)
return _berry_curvature(v_matrix, energy)
# This reverses the velocity unit (squared since Berry curvature is v.v)
_berry_curvature_const = 1 / _velocity_const**2
def _berry_curvature(v_M, energy):
r"""Calculate Berry curvature for a given velocity matrix"""
# All matrix elements along the 3 directions
N = v_M.shape[1]
# For cases where all states are degenerate then we would not be able
# to calculate anything. Hence we need to initialize as zero
# \Omega_{\alpha \beta, n}
sigma = zeros([3, 3, N], dtype=dtype_complex_to_real(v_M.dtype))
for s, e in enumerate(energy):
de = (energy - e) ** 2
# add factor 2 here, but omit the minus sign until later
# where we are forced to use the constant upon return anyways
np.divide(2, de, where=(de != 0), out=de)
# Calculate the berry-curvature
sigma[:, :, s] = ((de * v_M[:, s]) @ v_M[:, :, s].T).imag
# negative here
sigma *= -_berry_curvature_const
return sigma
[docs]
@set_module("sisl.physics.electron")
def conductivity(
bz,
distribution="fermi-dirac",
method="ahc",
degenerate=1.0e-5,
degenerate_dir=(1, 1, 1),
*,
eigenstate_kwargs=None,
):
r"""Electronic conductivity for a given `BrillouinZone` integral
Currently the *only* implemented method is the anomalous Hall conductivity (AHC, see :cite:`Wang2006`)
which may be calculated as:
.. math::
\sigma_{\alpha\beta} = \frac{-e^2}{\hbar}\int\,\mathrm d\mathbf k\sum_i f_i\Omega_{i,\alpha\beta}(\mathbf k)
where :math:`\Omega_{i,\alpha\beta}` and :math:`f_i` is the Berry curvature and occupation
for state :math:`i`.
The conductivity will be averaged by the Brillouin zone volume of the parent. See `BrillouinZone.volume` for details.
Hence for 1D the returned unit will be S/Ang, 2D it will be S/Ang^2 and 3D it will be S/Ang^3.
Parameters
----------
bz : BrillouinZone
containing the integration grid and has the ``bz.parent`` as an instance of Hamiltonian.
distribution : str or func, optional
distribution used to find occupations
method : {"ahc"}
"ahc" calculates the dc anomalous Hall conductivity
degenerate : float, optional
de-couple degenerate states within the given tolerance (in eV)
degenerate_dir : (3,), optional
along which direction degenerate states are decoupled.
eigenstate_kwargs : dict, optional
keyword arguments passed directly to the ``contour.eigenstate`` method.
One should *not* pass a ``k`` or a ``wrap`` keyword argument as they are
already used.
Returns
-------
cond : float
conductivity in units [S/cm^D]. The D is the dimensionality of the system.
See Also
--------
berry_curvature: method used to calculate the Berry-flux for calculating the conductivity
BrillouinZone.volume: volume calculation of the Brillouin zone
"""
from .hamiltonian import Hamiltonian
# Currently we require the conductivity calculation to *only* accept Hamiltonians
if not isinstance(bz.parent, Hamiltonian):
raise SislError(
"conductivity: requires the Brillouin zone object to contain a Hamiltonian!"
)
if isinstance(distribution, str):
distribution = get_distribution(distribution)
if eigenstate_kwargs is None:
eigenstate_kwargs = {}
method = method.lower()
if method == "ahc":
def _ahc(es):
occ = distribution(es.eig)
bc = es.berry_curvature(
degenerate=degenerate, degenerate_dir=degenerate_dir
)
return bc @ occ
vol, dim = bz.volume(ret_dim=True)
if dim == 0:
raise SislError(
f"conductivity: found a dimensionality of 0 which is non-physical"
)
cond = bz.apply.average.eigenstate(**eigenstate_kwargs, wrap=_ahc) * (
-constant.G0 / (4 * np.pi)
)
# Convert the dimensions from S/m^D to S/cm^D
cond /= vol * units(f"Ang^{dim}", f"cm^{dim}")
warn(
"conductivity: be aware that the units are currently not tested, please provide feedback!"
)
else:
raise SislError("conductivity: requires the method to be [ahc]")
return cond
[docs]
@set_module("sisl.physics.electron")
def berry_phase(
contour,
sub=None,
eigvals=False,
closed=True,
method="berry",
*,
eigenstate_kwargs=None,
ret_overlap=False,
):
r""" Calculate the Berry-phase on a loop path
The Berry phase for a single Bloch state is calculated using the discretized formula:
.. math::
\mathbf S = \prod_\alpha^{N-1} \langle \psi_{\mathbf k_\alpha} | \psi_{\mathbf k_{\alpha+1}} \rangle
\\
\phi = - \Im\ln \mathrm{det} \mathbf S
where :math:`\langle \psi_{\mathbf k_\alpha} | \psi_{\mathbf k_{\alpha+1}} \rangle` may be exchanged with an overlap matrix
of the investigated bands. I.e. :math:`\psi` is a manifold set of wavefunctions.
The overlap matrix :math:`\mathbf S` is also known as the global unitary
rotation matrix corresponding to the maximally localized Wannier centers.
If `closed` is true the overlap matrix will also include the circular inner product:
.. math::
\mathbf S^{\mathcal C} = \mathbf S \langle \psi_{\mathbf k_N} | \psi_{\mathbf k_1} \rangle
Parameters
----------
contour : BrillouinZone
containing the closed contour and has the ``contour.parent`` as an instance of Hamiltonian. The
first and last k-point must not be the same.
sub : None or list of int, optional
selected bands to calculate the Berry phase of
eigvals : bool, optional
return the eigenvalues of the product of the overlap matrices
closed : bool, optional
whether or not to include the connection of the last and first points in the loop
Forced true for Zak-phase calculations.
method : {"berry", "zak"}
"berry" will return the usual integral of the Berry connection over the specified contour
"zak" will compute the Zak phase for 1D systems by performing
a closed loop integration, see :cite:`Zak1989`.
Additionally, one may do the Berry-phase calculation using the SVD method of the
overlap matrices. Simply append ":svd" to the chosen method, e.g. "berry:svd".
eigenstate_kwargs : dict, optional
keyword arguments passed directly to the ``contour.eigenstate`` method.
One should *not* pass ``k`` as that is already used.
ret_overlap: bool, optional
optionally return the overlap matrix :math:`\mathbf S`
Notes
-----
The Brillouin zone object *need* not contain a closed discretized contour by doubling the first point.
The implementation is very similar to PythTB, except we are here using the :math:`\mathbf R` gauge
(convention II according to PythTB), see discussion in :pull:`131`.
For systems with band-crossings or degenerate states there is an arbitrariness to the definition
of the Berry phase for *individual* bands. However, the total phase (i.e., sum over filled bands) is
invariant and unaffected by this arbitrariness as long as the filled and empty bands do not intersect,
see :cite:`Resta2000`.
For non-orthogonal basis sets it is not fully known how important the :math:`\delta\mathbf k` spacing is since
it relies on the Lowdin transformation of the states. However, one should be careful about choosing
the correct bands for examination.
The returned angles are _not_ placed in the interval :math:`]-\pi;\pi]` as what `numpy.angle` would do.
This is to allow users to examine the quantities as is.
For more understanding of the Berry-phase and its calculation :cite:`TopInvTut` is a good reference.
Examples
--------
Calculate Berry-phase for first band but using the SVD method
>>> N = 30
>>> kR = 0.01
>>> normal = [0, 0, 1]
>>> origin = [1/3, 2/3, 0]
>>> bz = BrillouinZone.param_circle(H, N, kR, normal, origin)
>>> phase = berry_phase(bz, sub=0)
Calculate the multi-band Berry-phase using the SVD method, thus
ensuring removal of singular vectors.
>>> N = 30
>>> kR = 0.01
>>> normal = [0, 0, 1]
>>> origin = [1/3, 2/3, 0]
>>> bz = BrillouinZone.param_circle(H, N, kR, normal, origin)
>>> phase = berry_phase(bz, method="berry:svd")
"""
from .hamiltonian import Hamiltonian
# Currently we require the Berry phase calculation to *only* accept Hamiltonians
if not isinstance(contour.parent, Hamiltonian):
raise SislError(
"berry_phase: requires the Brillouin zone object to contain a Hamiltonian!"
)
if eigenstate_kwargs is None:
eigenstate_kwargs = {}
if contour.parent.orthogonal:
def _lowdin(state):
pass
else:
gauge = eigenstate_kwargs.get("gauge", "cell")
def _lowdin(state):
"""change state to the lowdin state, assuming everything is in R gauge
So needs to be done before changing gauge"""
S12 = sqrth(
state.parent.Sk(state.info["k"], gauge=gauge, format="array"),
overwrite_a=True,
)
state.state[:, :] = (S12 @ state.state.T).T
method, *opts = method.lower().split(":")
if method == "berry":
pass
elif method == "zak":
closed = True
else:
raise ValueError("berry_phase: requires the method to be [berry, zak]")
# We calculate the final angle from the determinant
_process = dot
if "svd" in opts:
def _process(prd, overlap):
U, _, V = svd_destroy(overlap)
return dot(prd, U @ V)
if sub is None:
def _berry(eigenstates):
# Grab the first one to be able to form a loop
first = next(eigenstates)
_lowdin(first)
# Create a variable to keep track of the previous state
prev = first
# Initialize the consecutive product
# Starting with the identity matrix!
prd = 1
# Loop remaining eigenstates
for second in eigenstates:
_lowdin(second)
prd = _process(prd, prev.inner(second, diag=False))
prev = second
# Complete the loop
if closed:
# Include last-to-first segment
prd = _process(prd, prev.inner(first, diag=False))
return prd
else:
def _berry(eigenstates):
nonlocal sub
first = next(eigenstates)
first.sub(sub, inplace=True)
_lowdin(first)
prev = first
prd = 1
for second in eigenstates:
second.sub(sub, inplace=True)
_lowdin(second)
prd = _process(prd, prev.inner(second, diag=False))
prev = second
if closed:
prd = _process(prd, prev.inner(first, diag=False))
return prd
S = _berry(contour.apply.iter.eigenstate(**eigenstate_kwargs))
# Get the angle of the berry-phase
# When using np.angle the returned value is in ]-pi; pi]
# However, small numerical differences makes wrap-arounds annoying.
# We'll always return the full angle. Then users can them-selves control
# how to convert them.
if eigvals:
ret = -log(la_eigvals(S, overwrite_a=not ret_overlap)).imag
ret = sort(ret)
else:
ret = -log(det(S, overwrite_a=not ret_overlap)).imag
if ret_overlap:
return ret, S
return ret
[docs]
@set_module("sisl.physics.electron")
def wavefunction(v, grid, geometry=None, k=None, spinor=0, spin=None, eta=None):
r"""Add the wave-function (`Orbital.psi`) component of each orbital to the grid
This routine calculates the real-space wave-function components in the
specified grid.
This is an *in-place* operation that *adds* to the current values in the grid.
It may be instructive to check that an eigenstate is normalized:
>>> grid = Grid(...)
>>> wavefunction(state, grid)
>>> (np.absolute(grid.grid) ** 2).sum() * grid.dvolume == 1.
Note: To calculate :math:`\psi(\mathbf r)` in a unit-cell different from the
originating geometry, simply pass a grid with a unit-cell smaller than the originating
supercell.
The wavefunctions are calculated in real-space via:
.. math::
\psi(\mathbf r) = \sum_i\phi_i(\mathbf r) |\psi\rangle_i \exp(i\mathbf k\cdot\mathbf R)
While for non-colinear/spin-orbit calculations the wavefunctions are determined from the
spinor component (`spinor`)
.. math::
\psi_{\alpha/\beta}(\mathbf r) = \sum_i\phi_i(\mathbf r) |\psi_{\alpha/\beta}\rangle_i \exp(i\mathbf k\cdot \mathbf R)
where ``spinor in [0, 1]`` determines :math:`\alpha` or :math:`\beta`, respectively.
Notes
-----
Currently this method only works for `v` being coefficients of the gauge="cell" method. In case
you are passing a `v` with the incorrect gauge you will find a phase-shift according to:
.. math::
\tilde v_j = e^{i\mathbf k\cdot\mathbf r_j} v_j
where :math:`j` is the orbital index and :math:`\mathbf r_j` is the orbital position.
Parameters
----------
v : array_like
coefficients for the orbital expansion on the real-space grid.
If `v` is a complex array then the `grid` *must* be complex as well. The coefficients
must be using the ``R`` gauge.
grid : Grid
grid on which the wavefunction will be plotted.
If multiple eigenstates are in this object, they will be summed.
geometry : Geometry, optional
geometry where the orbitals are defined. This geometry's orbital count must match
the number of elements in `v`.
If this is ``None`` the geometry associated with `grid` will be used instead.
k : array_like, optional
k-point associated with wavefunction, by default the inherent k-point used
to calculate the eigenstate will be used (generally shouldn't be used unless the `EigenstateElectron` object
has not been created via :meth:`~.Hamiltonian.eigenstate`).
spinor : int, optional
the spinor for non-colinear/spin-orbit calculations. This is only used if the
eigenstate object has been created from a parent object with a `Spin` object
contained, *and* if the spin-configuration is non-colinear or spin-orbit coupling.
Default to the first spinor component.
spin : Spin, optional
specification of the spin configuration of the orbital coefficients. This only has
influence for non-colinear wavefunctions where `spinor` choice is important.
eta : bool, optional
Display a console progressbar.
"""
# Decipher v from State type
if isinstance(v, State):
if geometry is None:
geometry = v._geometry()
if k is None:
k = v.info.get("k", k)
elif not np.allclose(k, v.info.get("k", k)):
raise ValueError(
f"wavefunction: k passed and k in info does not match: {k} and {v.info.get('k')}"
)
v = v.state
if geometry is None:
geometry = grid.geometry
if geometry is None:
raise SislError(
"wavefunction: did not find a usable Geometry through keywords or the Grid!"
)
# We cannot move stuff since outside stuff may rely on exact coordinates.
# If people have out-liers, they should do it them-selves.
# We'll do this and warn if they are dissimilar.
dxyz = geometry.lattice.cell2length(1e-6).sum(0)
dxyz = (
geometry.move(dxyz).translate2uc(axes=(0, 1, 2)).move(-dxyz).xyz - geometry.xyz
)
if not np.allclose(dxyz, 0):
info(
f"wavefunction: coordinates may be outside your primary unit-cell. "
"Translating all into the primary unit cell could disable this information"
)
# In case the user has passed several vectors we sum them to plot the summed state
if v.ndim == 2:
if v.shape[0] > 1:
info(
f"wavefunction: summing {v.shape[0]} different state coefficients, will continue silently!"
)
v = v.sum(0)
if spin is None:
if len(v) // 2 == geometry.no:
# We can see from the input that the vector *must* be a non-colinear calculation
v = v.reshape(-1, 2)[:, spinor]
info(
"wavefunction: assumes the input wavefunction coefficients to originate from a non-colinear calculation!"
)
elif spin.kind > Spin.POLARIZED:
# For non-colinear cases the user selects the spinor component.
v = v.reshape(-1, 2)[:, spinor]
if len(v) != geometry.no:
raise ValueError(
"wavefunction: require wavefunction coefficients corresponding to number of orbitals in the geometry."
)
# Check for k-points
k = _a.asarrayd(k)
kl = k.dot(k) ** 0.5
has_k = kl > 0.000001
if has_k:
info("wavefunction: k != Gamma is currently untested!")
# Check that input/grid makes sense.
# If the coefficients are complex valued, then the grid *has* to be
# complex valued.
# Likewise if a k-point has been passed.
is_complex = np.iscomplexobj(v) or has_k
if is_complex and not np.iscomplexobj(grid.grid):
raise SislError(
"wavefunction: input coefficients are complex, while grid only contains real."
)
if is_complex:
psi_init = _a.zerosz
else:
psi_init = _a.zerosd
# Extract sub variables used throughout the loop
shape = _a.asarrayi(grid.shape)
dcell = grid.dcell
ic_shape = grid.lattice.icell * shape.reshape(3, 1)
# Convert the geometry (hosting the wavefunction coefficients) coordinates into
# grid-fractionals X grid-shape to get index-offsets in the grid for the geometry
# supercell.
geom_shape = dot(geometry.cell, ic_shape.T)
# In the following we don't care about division
# So 1) save error state, 2) turn off divide by 0, 3) calculate, 4) turn on old error state
old_err = np.seterr(divide="ignore", invalid="ignore")
addouter = add.outer
def idx2spherical(ix, iy, iz, offset, dc, R):
"""Calculate the spherical coordinates from indices"""
rx = addouter(
addouter(ix * dc[0, 0], iy * dc[1, 0]), iz * dc[2, 0] - offset[0]
).ravel()
ry = addouter(
addouter(ix * dc[0, 1], iy * dc[1, 1]), iz * dc[2, 1] - offset[1]
).ravel()
rz = addouter(
addouter(ix * dc[0, 2], iy * dc[1, 2]), iz * dc[2, 2] - offset[2]
).ravel()
# Total size of the indices
n = rx.shape[0]
# Reduce our arrays to where the radius is "fine"
idx = indices_le(rx**2 + ry**2 + rz**2, R**2)
rx = rx[idx]
ry = ry[idx]
rz = rz[idx]
xyz_to_spherical_cos_phi(rx, ry, rz)
return n, idx, rx, ry, rz
# Figure out the max-min indices with a spacing of 1 radian
# calculate based on the minimum length of the grid-spacing
rad1 = pi / 180
theta, phi = ogrid[-pi:pi:rad1, 0:pi:rad1]
cphi, sphi = cos(phi), sin(phi)
ctheta_sphi = cos(theta) * sphi
stheta_sphi = sin(theta) * sphi
del sphi
nrxyz = (theta.size, phi.size, 3)
del theta, phi, rad1
# First we calculate the min/max indices for all atoms
rxyz = _a.emptyd(nrxyz)
rxyz[..., 0] = ctheta_sphi
rxyz[..., 1] = stheta_sphi
rxyz[..., 2] = cphi
# Reshape
rxyz.shape = (-1, 3)
idx = dot(rxyz, ic_shape.T)
idxm = idx.min(0)
idxM = idx.max(0)
del ctheta_sphi, stheta_sphi, cphi, idx, rxyz, nrxyz
# Fast loop (only per specie)
origin = grid.lattice.origin
idx_mm = _a.emptyd([geometry.na, 2, 3])
all_negative_R = True
for atom, ia in geometry.atoms.iter(True):
if len(ia) == 0:
continue
R = atom.maxR()
all_negative_R = all_negative_R and R < 0.0
# Now do it for all the atoms to get indices of the middle of
# the atoms
# The coordinates are relative to origin, so we need to shift (when writing a grid
# it is with respect to origin)
idx = dot(geometry.xyz[ia, :] - origin, ic_shape.T)
# Get min-max for all atoms
idx_mm[ia, 0, :] = idxm * R + idx
idx_mm[ia, 1, :] = idxM * R + idx
if all_negative_R:
raise SislError(
"wavefunction: Cannot create wavefunction since no atoms have an associated basis-orbital on a real-space grid"
)
# Now we have min-max for all atoms
# When we run the below loop all indices can be retrieved by looking
# up in the above table.
# Before continuing, we can easily clean up the temporary arrays
del origin, idx
arangei = _a.arangei
# In case this grid does not have a Geometry associated
# We can *perhaps* easily attach a geometry with the given
# atoms in the unit-cell
lattice = grid.lattice.copy()
# Find the periodic directions
pbc = [
bc == BC.PERIODIC or geometry.nsc[i] > 1
for i, bc in enumerate(grid.lattice.boundary_condition[:, 0])
]
if grid.geometry is None:
# Create the actual geometry that encompass the grid
ia, xyz, _ = geometry.within_inf(lattice, periodic=pbc)
if len(ia) > 0:
grid.set_geometry(Geometry(xyz, geometry.atoms[ia], lattice=lattice))
# Instead of looping all atoms in the supercell we find the exact atoms
# and their supercell indices.
# plus some tolerance
add_R = _a.fulld(3, geometry.maxR()) + 1.0e-6
# Calculate the required additional vectors required to increase the fictitious
# supercell by add_R in each direction.
# For extremely skewed lattices this will be way too much, hence we make
# them square.
o = lattice.to.Cuboid(orthogonal=True)
lattice = Lattice(o._v + np.diag(2 * add_R), origin=o.origin - add_R)
# Retrieve all atoms within the grid supercell
# (and the neighbors that connect into the cell)
# Note that we cannot pass the "moved" origin because then ISC would be wrong
IA, XYZ, ISC = geometry.within_inf(lattice, periodic=pbc)
# We need to revert the grid supercell origin as that is not subtracted in the `within_inf` returned
# coordinates (and the below loop expects positions with respect to the origin of the plotting
# grid).
XYZ -= grid.lattice.origin
phk = k * 2 * np.pi
phase = 1
# Retrieve progressbar
eta = progressbar(len(IA), "wavefunction", "atom", eta)
# Loop over all atoms in the grid-cell
for ia, xyz, isc in zip(IA, XYZ, ISC):
# Get current atom
atom = geometry.atoms[ia]
# Extract maximum R
R = atom.maxR()
if R <= 0.0:
warn(
f"wavefunction: Atom '{atom}' does not have a wave-function, skipping atom."
)
eta.update()
continue
# Get indices in the supercell grid
idx = (isc.reshape(3, 1) * geom_shape).sum(0)
idxm = floor(idx_mm[ia, 0, :] + idx).astype(int32)
idxM = ceil(idx_mm[ia, 1, :] + idx).astype(int32) + 1
# Fast check whether we can skip this point
if (
idxm[0] >= shape[0]
or idxm[1] >= shape[1]
or idxm[2] >= shape[2]
or idxM[0] <= 0
or idxM[1] <= 0
or idxM[2] <= 0
):
eta.update()
continue
# Truncate values
if idxm[0] < 0:
idxm[0] = 0
if idxM[0] > shape[0]:
idxM[0] = shape[0]
if idxm[1] < 0:
idxm[1] = 0
if idxM[1] > shape[1]:
idxM[1] = shape[1]
if idxm[2] < 0:
idxm[2] = 0
if idxM[2] > shape[2]:
idxM[2] = shape[2]
# Now idxm/M contains min/max indices used
# Convert to spherical coordinates
n, idx, r, theta, phi = idx2spherical(
arangei(idxm[0], idxM[0]),
arangei(idxm[1], idxM[1]),
arangei(idxm[2], idxM[2]),
xyz,
dcell,
R,
)
# Get initial orbital
io = geometry.a2o(ia)
if has_k:
phase = exp(1j * phk.dot(isc))
# Allocate a temporary array where we add the psi elements
psi = psi_init(n)
# Loop on orbitals on this atom, grouped by radius
for os in atom.iter(True):
# Get the radius of orbitals (os)
oR = os[0].R
if oR <= 0.0:
warn(
f"wavefunction: Orbital(s) '{os}' does not have a wave-function, skipping orbital!"
)
# Skip these orbitals
io += len(os)
continue
# Downsize to the correct indices
if R - oR < 1e-6:
idx1 = idx
r1 = r
theta1 = theta
phi1 = phi
else:
idx1 = indices_le(r, oR)
# Reduce arrays
r1 = r[idx1]
theta1 = theta[idx1]
phi1 = phi[idx1]
idx1 = idx[idx1]
# Loop orbitals with the same radius
for o in os:
# Evaluate psi component of the wavefunction and add it for this atom
psi[idx1] += o.psi_spher(r1, theta1, phi1, cos_phi=True) * (
v[io] * phase
)
io += 1
# Clean-up
del idx1, r1, theta1, phi1, idx, r, theta, phi
# Convert to correct shape and add the current atom contribution to the wavefunction
psi.shape = idxM - idxm
grid.grid[idxm[0] : idxM[0], idxm[1] : idxM[1], idxm[2] : idxM[2]] += psi
# Clean-up
del psi
# Step progressbar
eta.update()
eta.close()
# Reset the error code for division
np.seterr(**old_err)
[docs]
class _electron_State:
# pylint: disable=E1101
__slots__ = []
def __is_nc(self, spin: Optional[Spin] = None):
"""Internal routine to check whether this is a non-colinear calculation"""
if spin is not None:
return not spin.is_diagonal
try:
return not self.parent.spin.is_diagonal
except Exception:
return False
[docs]
def Sk(self, format=None):
r"""Retrieve the overlap matrix corresponding to the originating parent structure.
When ``self.parent`` is a Hamiltonian this will return :math:`\mathbf S(\mathbf k)` for the
:math:`\mathbf k`-point these eigenstates originate from.
Parameters
----------
format : str, optional
the returned format of the overlap matrix. This only takes effect for
non-orthogonal parents.
"""
if format is None:
format = self.info.get("format", "csr")
if isinstance(self.parent, SparseOrbitalBZSpin):
# Calculate the overlap matrix
if not self.parent.orthogonal:
opt = {
"k": self.info.get("k", (0, 0, 0)),
"dtype": self.dtype,
"format": format,
}
for key in ("gauge",):
val = self.info.get(key, None)
if not val is None:
opt[key] = val
return self.parent.Sk(**opt)
n = m = self.shape[1]
if "sc:" in format:
m = n * self.parent.n_s
return _FakeMatrix(n, m)
[docs]
@deprecate_argument(
"sum",
"projection",
"argument sum has been deprecated in favor of projection",
"0.15",
"0.16",
)
def norm2(self, projection: Literal["sum", "orbitals", "basis", "atoms"] = "sum"):
r"""Return a vector with the norm of each state :math:`\langle\psi|\mathbf S|\psi\rangle`
:math:`\mathbf S` is the overlap matrix (or basis), for orthogonal basis
:math:`\mathbf S \equiv \mathbf I`.
Parameters
----------
projection :
whether to compute the norm per state as a single number or as orbital-/atom-resolved quantity
See Also
--------
inner: used method for calculating the squared norm.
Returns
-------
numpy.ndarray
the squared norm for each state
"""
return self.inner(matrix=self.Sk(), projection=projection)
[docs]
def spin_moment(self, project=False):
r"""Calculate spin moment from the states
This routine calls `~sisl.physics.electron.spin_moment` with appropriate arguments
and returns the spin moment for the states.
See `~sisl.physics.electron.spin_moment` for details.
Parameters
----------
project : bool, optional
whether the moments are orbitally resolved or not
"""
return spin_moment(self.state, self.Sk(), project=project)
[docs]
def wavefunction(self, grid, spinor=0, eta=None):
r"""Expand the coefficients as the wavefunction on `grid` *as-is*
See `~sisl.physics.electron.wavefunction` for argument details, the arguments not present
in this method are automatically passed from this object.
"""
spin = getattr(self.parent, "spin", None)
if isinstance(self.parent, Geometry):
geometry = self.parent
else:
geometry = getattr(self.parent, "geometry", None)
if not isinstance(grid, Grid):
# probably the grid is a Real, or a tuple that denotes the shape
# at least this makes it easier to parse
grid = Grid(grid, geometry=geometry, dtype=self.dtype)
# Ensure we are dealing with the R gauge
self.change_gauge("R")
# Retrieve k
k = self.info.get("k", _a.zerosd(3))
wavefunction(
self.state, grid, geometry=geometry, k=k, spinor=spinor, spin=spin, eta=eta
)
[docs]
@set_module("sisl.physics.electron")
class CoefficientElectron(Coefficient):
r"""Coefficients describing some physical quantity related to electrons"""
__slots__ = []
[docs]
@set_module("sisl.physics.electron")
class StateElectron(_electron_State, State):
r"""A state describing a physical quantity related to electrons"""
__slots__ = []
[docs]
@set_module("sisl.physics.electron")
class StateCElectron(_electron_State, StateC):
r"""A state describing a physical quantity related to electrons, with associated coefficients of the state"""
__slots__ = []
[docs]
def velocity(self, *args, **kwargs):
r"""Calculate velocity for the states
This routine calls ``derivative(1, *args, **kwargs)`` and returns the velocity for the states.
Note that the coefficients associated with the `StateCElectron` *must* correspond
to the energies of the states.
Notes
-----
The states and energies for the states *may* have changed after calling this routine.
This is because of the velocity un-folding for degenerate modes. I.e. calling
`PDOS` after this method *may* change the result.
The velocities are calculated without the Berry curvature contribution see Eq. (2) in :cite:`Wang2006`.
The missing contribution may be added in later editions, for completeness sake, it is:
.. math::
\delta \mathbf v = - \mathbf k\times \Omega_i(\mathbf k)
where :math:`\Omega_i` is the Berry curvature for state :math:`i`.
See Also
--------
derivative : for details of the implementation
"""
v = self.derivative(1, *args, **kwargs)
v *= _velocity_const
return v
[docs]
def berry_curvature(self, *args, **kwargs):
r"""Calculate Berry curvature for the states
This routine calls ``derivative(1, *args, **kwargs, matrix=True)`` and
returns the Berry curvature for the states.
Note that the coefficients associated with the `StateCElectron` *must* correspond
to the energies of the states.
See Also
--------
derivative : for details of the velocity matrix calculation implementation
sisl.physics.electron.berry_curvature : for details of the Berry curvature implementation
"""
v = self.derivative(1, *args, **kwargs, matrix=True)
return _berry_curvature(v, self.c)
[docs]
def effective_mass(self, *args, **kwargs):
r"""Calculate effective mass tensor for the states, units are (ps/Ang)^2
This routine calls ``derivative(2, *args, **kwargs)`` and
returns the effective mass for all states.
Note that the coefficients associated with the `StateCElectron` *must* correspond
to the energies of the states.
Notes
-----
Since some directions may not be periodic there will be zeros. This routine will
invert elements where the values are different from 0.
It is not advisable to use a `sub` before calculating the effective mass
since the 1st order perturbation uses the energy differences and the 1st derivative
matrix for correcting the curvature.
The returned effective mass is given in the Voigt notation.
For :math:`\Gamma` point calculations it may be beneficial to pass `dtype=np.complex128`
to the `eigenstate` argument to ensure their complex values. This is necessary for the
degeneracy decoupling.
See Also
--------
derivative: for details of the implementation
"""
ieff = self.derivative(2, *args, **kwargs)[1].real
np.divide(_velocity_const**2, ieff, where=(ieff != 0), out=ieff)
return ieff
[docs]
@set_module("sisl.physics.electron")
class EigenvalueElectron(CoefficientElectron):
r"""Eigenvalues of electronic states, no eigenvectors retained
This holds routines that enable the calculation of density of states.
"""
__slots__ = []
@property
def eig(self):
"""Eigenvalues"""
return self.c
[docs]
def occupation(self, distribution="fermi_dirac"):
r"""Calculate the occupations for the states according to a distribution function
Parameters
----------
distribution : str or func, optional
distribution used to find occupations
Returns
-------
numpy.ndarray
``len(self)`` with occupation values
"""
if isinstance(distribution, str):
distribution = get_distribution(distribution)
return distribution(self.eig)
[docs]
def DOS(self, E, distribution="gaussian"):
r"""Calculate DOS for provided energies, `E`.
This routine calls `sisl.physics.electron.DOS` with appropriate arguments
and returns the DOS.
See `~sisl.physics.electron.DOS` for argument details.
"""
return DOS(E, self.eig, distribution)
[docs]
@set_module("sisl.physics.electron")
class EigenvectorElectron(StateElectron):
r"""Eigenvectors of electronic states, no eigenvalues retained
This holds routines that enable the calculation of spin moments.
"""
__slots__ = []
[docs]
@set_module("sisl.physics.electron")
class EigenstateElectron(StateCElectron):
r"""Eigen states of electrons with eigenvectors and eigenvalues.
This holds routines that enable the calculation of (projected) density of states,
spin moments (spin texture).
"""
__slots__ = []
@property
def eig(self):
r"""Eigenvalues for each state"""
return self.c
[docs]
def occupation(self, distribution="fermi_dirac"):
r"""Calculate the occupations for the states according to a distribution function
Parameters
----------
distribution : str or func, optional
distribution used to find occupations
Returns
-------
numpy.ndarray
``len(self)`` with occupation values
"""
if isinstance(distribution, str):
distribution = get_distribution(distribution)
return distribution(self.eig)
[docs]
def DOS(self, E, distribution="gaussian"):
r"""Calculate DOS for provided energies, `E`.
This routine calls `sisl.physics.electron.DOS` with appropriate arguments
and returns the DOS.
See `~sisl.physics.electron.DOS` for argument details.
"""
return DOS(E, self.c, distribution)
[docs]
def PDOS(self, E, distribution="gaussian"):
r"""Calculate PDOS for provided energies, `E`.
This routine calls `~sisl.physics.electron.PDOS` with appropriate arguments
and returns the PDOS.
See `~sisl.physics.electron.PDOS` for argument details.
"""
return PDOS(
E,
self.c,
self.state,
self.Sk(),
distribution,
getattr(self.parent, "spin", None),
)
[docs]
def COP(self, E, M, *args, **kwargs):
r"""Calculate COP for provided energies, `E` using matrix `M`
This routine calls `~sisl.physics.electron.COP` with appropriate arguments.
"""
return COP(E, self.c, self.state, M, *args, **kwargs)
[docs]
def COOP(self, E, *args, **kwargs):
r"""Calculate COOP for provided energies, `E`.
This routine calls `~sisl.physics.electron.COP` with appropriate arguments.
"""
format = self.info.get("format", "csr")
Sk = self.Sk(format=f"sc:{format}")
return COP(E, self.c, self.state, Sk, *args, **kwargs)
[docs]
def COHP(self, E, *args, **kwargs):
r"""Calculate COHP for provided energies, `E`.
This routine calls `~sisl.physics.electron.COP` with appropriate arguments.
"""
format = self.info.get("format", "csr")
Hk = self.parent.Hk(format=f"sc:{format}")
return COP(E, self.c, self.state, Hk, *args, **kwargs)