Add new occupation_numbers function, docs and test

parent b1902798
Pipeline #9466451 passed with stage
in 7 minutes and 27 seconds
Miscellaneous objects
=====================
Miscellaneous objects and functions
===================================
.. autofunction:: gpaw.occupations.occupation_numbers
.. autoclass:: gpaw.lfc.NewLocalizedFunctionsCollection
:members:
......
......@@ -41,6 +41,10 @@ Git master branch
* Calculation of vectorial magnetic moments inside PAW spheres based on
spin-orbit spinors.
* Added a simple :func:`gpaw.occupations.occupation_numbers` function for
calculating occupation numbers, fermi-level, mognetic moment, and entropy
from eigenvalues and k-point weights.
Version 1.2.0
=============
......
......@@ -37,7 +37,7 @@ class KPoint:
Rank of the CPU that does the matrix diagonalization of
H_nn and the Cholesky decomposition of S_nn.
"""
def __init__(self, weight, s, k, q, phase_cd):
"""Construct k-point object.
......@@ -79,7 +79,7 @@ class KPoint:
self.k = k # k-point index
self.q = q # local k-point index
self.phase_cd = phase_cd
self.eps_n = None
self.f_n = None
self.P_ani = None
......@@ -89,7 +89,7 @@ class KPoint:
self.C_nM = None # LCAO coefficients for wave functions XXX
self.rho_MM = None
self.S_MM = None
self.T_MM = None
......@@ -154,4 +154,3 @@ class GlobalKPoint(KPoint):
ni = wfs.setups[a].ni
self.P_ani[a] = my_P_ni[:, i:i + ni] # copy?
i += ni
......@@ -88,6 +88,7 @@ def kpts2sizeandoffsets(size=None, density=None, gamma=None, even=None,
return size, offsets
class KPointDescriptor:
"""Descriptor-class for k-points."""
......
......@@ -22,6 +22,70 @@ def create_occupation_number_object(name, **kwargs):
raise ValueError('Unknown occupation number object name: ' + name)
def occupation_numbers(occ, eps_skn, weight_k, nelectrons):
"""Calculate occupation numbers from eigenvalues.
occ: dict
Example: {'name': 'fermi-dirac', 'width': 0.05}.
eps_skn: ndarray, shape=(nspins, nibzkpts, nbands)
Eigenvalues.
weight_k: ndarray, shape=(nibzkpts,)
Weights of k-points in IBZ (must sum to 1).
nelectrons: int or float
Number of electrons.
Returns a tuple containing:
* f_skn (sums to nelectrons)
* fermi-level
* magnetic moment
* entropy as -S*T
"""
from gpaw.grid_descriptor import GridDescriptor
from gpaw.mpi import serial_comm
occ = create_occupation_number_object(**occ)
eps_skn = np.asarray(eps_skn)
weight_k = np.asarray(weight_k)
nspins, nkpts, nbands = eps_skn.shape
f_skn = np.empty_like(eps_skn)
class SimpleNamespace:
"""Same as types.SimpleNamespace from Python 3.3+."""
def __init__(self, **kwargs):
self.__dict__.update(kwargs)
wfs = SimpleNamespace(kpt_u=[],
nvalence=nelectrons,
nspins=nspins,
kptband_comm=serial_comm,
world=serial_comm,
gd=GridDescriptor([4, 4, 4], [1.0, 1.0, 1.0]),
bd=SimpleNamespace(nbands=nbands,
collect=lambda x: x,
comm=serial_comm),
kd=SimpleNamespace(mynks=nspins * nkpts,
comm=serial_comm,
nspins=nspins,
nibzkpts=nkpts,
weight_k=weight_k,
collect=lambda x, broadcast: x))
for s in range(nspins):
for k in range(nkpts):
kpt = SimpleNamespace(s=s,
weight=weight_k[k] * 2 / nspins,
eps_n=eps_skn[s, k],
f_n=f_skn[s, k])
wfs.kpt_u.append(kpt)
occ.calculate(wfs)
return f_skn, occ.fermilevel, occ.magmom, occ.e_entropy
def findroot(func, x, tol=1e-10):
"""Function used for locating Fermi level."""
xmin = -np.inf
......@@ -451,7 +515,7 @@ class SmoothDistribution(ZeroKelvin):
data = np.empty(4)
def f(x, data=data):
data.fill(0.0)
data[:] = 0.0
for kpt in wfs.kpt_u:
if kpt.s in spins:
data += self.distribution(kpt, x)
......
from __future__ import print_function
import numpy as np
from ase.units import Hartree
from gpaw.occupations import FermiDirac, MethfesselPaxton
from gpaw.occupations import FermiDirac, MethfesselPaxton, occupation_numbers
class KPoint:
......@@ -9,7 +9,8 @@ class KPoint:
f_n = np.empty(1)
weight = 0.2
s = 0
k = KPoint()
......@@ -18,7 +19,7 @@ def f(occ, x):
n, dnde, x, S = occ.distribution(k, 0.0)
return n, dnde, S
def test(occ):
print(occ)
for e in [-0.3 / Hartree, 0, 0.1 / Hartree, 1.2 / Hartree]:
......@@ -33,7 +34,19 @@ def test(occ):
assert abs(d - d0) < 3e-5
assert abs(dS - e * dn) < 1e-13
for w in [0.1, 0.5]:
test(FermiDirac(w))
for n in range(4):
test(MethfesselPaxton(w, n))
occ = {'name': 'fermi-dirac', 'width': 0.1}
for eps_skn, weight_k, n in [
([[[0.0, 1.0]]], [1.0], 2),
([[[0.0, 1.0]], [[0.0, 2.0]]], [1.0], 3),
([[[0.0, 1.0, 2.0], [0.0, 2.0, 2.0]]], [0.5, 0.5], 2)]:
f_skn, fl, m, s = occupation_numbers(occ, eps_skn, weight_k, n)
print(f_skn, fl, m, s)
assert abs(f_skn.sum() - n) < 1e-14, f_skn
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment