Introduction to GPAW internals

This guide will contain graphs showing the relationship between objects that build up a DFT calculation engine.

Hint

Here is a simple graph showing the relations between the classes A, B and C:

../_images/abc.svg

Here, the object of type A has an attribute a of type int and an attribute b of type B or C, where C inherits from B.

DFT components

The components needed for a DFT calculation are created by a “builder” that can be made with the builder() function, an ASE ase.Atoms object and some input parameters:

>>> from ase import Atoms
>>> atoms = Atoms('Li', cell=[2, 2, 2], pbc=True)
>>> from gpaw.new.builder import builder
>>> params = {'mode': 'pw', 'kpts': (5, 5, 5)}
>>> b = builder(atoms, params)
../_images/builder.svg

As seen in the figure above, there are builders for each of the modes: PW, FD and LCAO (builders for TB and ATOM modes are not shown).

The InputParameters object takes care of user parameters:

  • checks for errors

  • does normalization

  • handles backwards compatibility and deprecation warnings

Normally, you will not need to create a DFT-components builder yourself. It will happen automatically when you create a DFT-calculation object like this:

>>> from gpaw.new.calculation import DFTCalculation
>>> calculation = DFTCalculation.from_parameters(atoms, params)

or when you create an ASE-calculator interface:

>>> from gpaw.new.ase_interface import GPAW
>>> atoms.calc = GPAW(**params, txt='li.txt')

Full picture

The ase.Atoms object has an gpaw.new.ase_interface.ASECalculator object attached created with the gpaw.new.ase_interface.GPAW() function:

>>> atoms = Atoms('H2',
...               positions=[(0, 0, 0), (0, 0, 0.75)],
...               cell=[2, 2, 3],
...               pbc=True)
>>> atoms.calc = GPAW(mode='pw', txt='h2.txt')
>>> atoms.calc
ASECalculator(mode: {'name': 'pw'})

The atoms.calc object manages a gpaw.new.calculation.DFTCalculation object that does the actual work. When we do this:

>>> e = atoms.get_potential_energy()

the gpaw.new.ase_interface.ASECalculator.get_potential_energy() method gets called (atoms.calc.get_potential_energy(atoms)) and the following will happen:

  • create gpaw.new.calculation.DFTCalculation object if not already done

  • update positions/unit cell if they have changed

  • start SCF loop and converge if needed

  • calculate energy

  • store a copy of the atoms

../_images/code.svg

DFT-calculation object

An instance of the gpaw.new.calculation.DFTCalculation class has the following attributes:

state

gpaw.new.calculation.DFTState

scf_loop

gpaw.new.scf.SCFLoop

pot_calc

gpaw.new.pot_calc.PotentialCalculator

and a the gpaw.new.calculation.DFTState object has these attributes:

density

gpaw.new.density.Density

ibzwfs

gpaw.new.ibzwfs.IBZWaveFunctions

potential

gpaw.new.potential.Potential

Naming convention for arrays

Commonly used indices:

index

description

a

Atom number

c

Unit cell axis-index (0, 1, 2)

v

xyz-index (0, 1, 2)

K

BZ k-point index

k

IBZ k-point index

q

IBZ k-point index (local, i.e. it starts at 0 on each processor)

s

Spin index (\(\sigma\))

s

Symmetry index

u

Combined spin and k-point index (local)

R

Three indices into the coarse 3D grid

r

Three indices into the fine 3D grid

G

Index of plane-wave coefficient (wave function expansion, ecut)

g

Index of plane-wave coefficient (densities, 2 * ecut)

h

Index of plane-wave coefficient (compensation charges, 8 * ecut)

X

R or G

x

r, g or h

x

Zero or more extra dimensions

M

LCAO orbital index (\(\mu\))

n

Band number

n

Principal quantum number

l

Angular momentum quantum number (s, p, d, …)

m

Magnetic quantum number (0, 1, …, 2*`ell` - 1)

L

l and m (L = l**2 + m)

j

Valence orbital number (n and l)

i

Valence orbital number (n, l and m)

q

j1 and j2 pair

p

i1 and i2 pair

r

CPU-rank

Examples:

density.D_asii

\(D_{\sigma,i_1,i_2}^a\)

AtomArrays

density.nt_sR

\(\tilde{n}_\sigma(\mathbf{r})\)

UGArray

ibzwfs.wfs_qs[q][s].P_ain

\(P_{\sigma \mathbf{k} in}^a\)

AtomArrays

ibzwfs.wfs_qs[q][s].psit_nX

\(\tilde{\psi}_{\sigma \mathbf{k} n}(\mathbf{r})\)

UGArray | PWArray

ibzwfs.wfs_qs[q][s].pt_aX

\(\tilde{p}_{\sigma \mathbf{k} i}^a(\mathbf{r}-\mathbf{R}^a)\)

AtomCenteredFunctions

Domain descriptors

GPAW has two different container types for storing one or more functions in a unit cell (wave functions, electron densities, …):

../_images/da.svg

Uniform grids

A uniform grid can be created with the UGDesc class:

>>> import numpy as np
>>> from gpaw.core import UGDesc
>>> a = 4.0
>>> n = 20
>>> grid = UGDesc(cell=a * np.eye(3),
...               size=(n, n, n))

Given a UGDesc object, one can create UGArray objects like this

>>> func_R = grid.empty()
>>> func_R.data.shape
(20, 20, 20)
>>> func_R.data[:] = 1.0
>>> grid.zeros((3, 2)).data.shape
(3, 2, 20, 20, 20)

Here are the methods of the UGDesc class:

size()

Size of uniform grid.

global_shape()

Actual size of uniform grid.

phase_factor_cd()

Phase factor for block-boundary conditions.

new()

Create new uniform grid description.

empty()

Create new UGArray object.

from_data()

blocks()

Yield views of blocks of data.

xyz()

Create array of (x, y, z) coordinates.

atom_centered_functions()

Create UGAtomCenteredFunctions object.

transformer()

Create transformer from one grid to another.

eikr()

Plane wave.

from_cell_and_grid_spacing()

Create UGDesc from grid-spacing.

fft_plans()

Create FFTW-plans.

ranks_from_fractional_positions()

ecut_max()

and the UGArray class:

new()

Create new UniforGridFunctions object of same kind.

xy()

Extract x, y values along line.

scatter_from()

Scatter data from rank-0 to all ranks.

gather()

Gather data from all ranks to rank-0.

fft()

Do FFT.

norm2()

Calculate integral over cell of absolute value squared.

integrate()

Integral of self or self times cc(other).

to_pbc_grid()

Convert to UniformGrud with pbc=(True, True, True).

multiply_by_eikr()

Multiply by \(exp(ik.r)\).

interpolate()

Interpolate to finer grid.

fft_restrict()

Restrict to coarser grid.

abs_square()

Add weighted absolute square of data to output array.

symmetrize()

Make data symmetric.

randomize()

Insert random numbers between -0.5 and 0.5 into data.

moment()

Calculate moment of data.

scaled()

Create new scaled UGArray object.

add_ked()

redist()

Plane waves

A set of plane-waves are characterized by a cutoff energy and a uniform grid:

>>> from gpaw.core import PWDesc
>>> pw = PWDesc(ecut=100, cell=grid.cell)
>>> func_G = pw.empty()
>>> func_R.fft(out=func_G)
PWArray(pw=PWDesc(ecut=100 <coefs=1536/1536>, cell=[4.0, 4.0, 4.0], pbc=[True, True, True], comm=0/1, dtype=float64), dims=())
>>> G = pw.reciprocal_vectors()
>>> G.shape
(1536, 3)
>>> G[0]
array([0., 0., 0.])
>>> func_G.data[0]
(1+0j)
>>> func_G.ifft(out=func_R)
UGArray(grid=UGDesc(size=[20, 20, 20], cell=[4.0, 4.0, 4.0], pbc=[True, True, True], comm=0/1, dtype=float64), dims=())
>>> round(func_R.data[0, 0, 0], 15)
1.0

Here are the methods of the PWDesc class:

itemsize()

int([x]) -> integer

global_shape()

Tuple with one element: number of plane waves.

reciprocal_vectors()

Returns reciprocal lattice vectors, G + k, in xyz coordinates.

kinetic_energies()

Kinetic energy of plane waves.

empty()

Create new PlaneWaveExpanions object.

from_data()

new()

Create new plane-wave expansion description.

indices()

Return indices into FFT-grid.

cut()

Cut out G-vectors with (G+k)^2/2<E_kin.

paste()

Paste G-vectors with (G+k)^2/2<E_kin into 3-D FFT grid and

map_indices()

Map from one (distributed) set of plane waves to smaller global set.

atom_centered_functions()

Create PlaneWaveAtomCenteredFunctions object.

and the PWArray class:

new()

Create new PWArray object of same kind.

copy()

Create a copy (surprise!).

matrix()

Matrix view of data.

ifft()

Do inverse FFT(s) to uniform grid(s).

interpolate()

gather()

Gather coefficients on master.

gather_all()

Gather coefficients from self[r] on rank r.

scatter_from()

Scatter data from rank-0 to all ranks.

scatter_from_all()

Scatter all coefficients from rank r to self on other cores.

integrate()

Integral of self or self time cc(other).

norm2()

Calculate integral over cell.

abs_square()

Add weighted absolute square of self to output array.

to_pbc_grid()

randomize()

Insert random numbers between -0.5 and 0.5 into data.

moment()

morph()

add_ked()

Atoms-arrays

../_images/aa.svg

Block boundary conditions

Matrix elements

>>> psit_nG = pw.zeros(5)
>>> def T(psit_nG):
...     """Kinetic energy operator."""
...     out = psit_nG.new()
...     out.data[:] = psit_nG.desc.ekin_G * psit_nG.data
...     return out
>>> H_nn = psit_nG.matrix_elements(psit_nG, function=T)

Same as:

>>> Tpsit_nG = T(psit_nG)
>>> psit_nG.matrix_elements(Tpsit_nG, symmetric=True)
Matrix(float64: 5x5)

but faster.

Atom-centered functions

../_images/acf.svg
# creates: acf_example.png
import numpy as np
import matplotlib.pyplot as plt
from gpaw.core import UGDesc


alpha = 4.0
rcut = 2.0
l = 0
gauss = (l, rcut, lambda r: (4 * np.pi)**0.5 * np.exp(-alpha * r**2))
grid = UGDesc(cell=[4.0, 2.0, 2.0], size=[40, 20, 20])
pos = [[0.25, 0.5, 0.5], [0.75, 0.5, 0.5]]
acf_aR = grid.atom_centered_functions([[gauss], [gauss]], pos)
coef_as = acf_aR.empty(dims=(2,))
coef_as[0] = [[1], [-1]]
coef_as[1] = [[2], [1]]
print(coef_as.data, coef_as[0])
f_sR = grid.zeros(2)
acf_aR.add_to(f_sR, coef_as)
x = grid.xyz()[:, 10, 10, 0]
y1, y2 = f_sR.data[:, :, 10, 10]
ax = plt.subplot(1, 1, 1)
ax.plot(x, y1, 'o-')
ax.plot(x, y2, 'x-')
# plt.show()
plt.savefig('acf_example.png')
../_images/acf_example.png

Matrix object

Here are the methods of the Matrix class:

new()

Create new matrix of same shape and dtype.

copy()

Create a copy.

multiply()

BLAS matrix-multiplication with other matrix.

redist()

Redistribute to other BLACS layout.

gather()

Gather the Matrix on the root rank.

inv()

Inplace inversion.

invcholesky()

In-place inverse of Cholesky decomposition.

eigh()

Calculate eigenvectors and eigenvalues.

eighg()

Solve generalized eigenvalue problem.

complex_conjugate()

Inplace complex conjugation.

add_hermitian_conjugate()

Add hermitian conjugate to myself.

tril2full()

Fill in upper triangle from lower triangle.

add_to_diagonal()

Add list of numbers or single number to diagonal of matrix.

to_cpu()

A simple example that we can run with MPI on 4 cores:

from gpaw.core.matrix import Matrix
from gpaw.mpi import world
a = Matrix(5, 5, dist=(world, 2, 2, 2))
a.data[:] = world.rank
print(world.rank, a.data.shape)

Here, we have created a 5x5 Matrix of floats distributed on a 2x2 BLACS grid with a block size of 2 and we then print the shapes of the ndarrays, which looks like this (in random order):

1 (2, 3)
2 (3, 2)
3 (2, 2)
0 (3, 3)

Let’s create a new matrix b and redistribute from a to b:

b = a.new(dist=(None, 1, 1, None))
a.redist(b)
if world.rank == 0:
    print(b.array)

This will output:

[[ 0.  0.  2.  2.  0.]
 [ 0.  0.  2.  2.  0.]
 [ 1.  1.  3.  3.  1.]
 [ 1.  1.  3.  3.  1.]
 [ 0.  0.  2.  2.  0.]]

Matrix-matrix multiplication works like this:

c = a.multiply(a, opb='T')

API

Core

class gpaw.core.UGDesc(*, cell, size, pbc=(True, True, True), kpt=None, comm=MPIComm(size=1, rank=0), decomp=None, dtype=None)[source]

Description of 3D uniform grid.

Parameters:
  • cell (ArrayLike1D | ArrayLike2D) – Unit cell given as three floats (orthorhombic grid), six floats (three lengths and the angles in degrees) or a 3x3 matrix (units: bohr).

  • size (ArrayLike1D) – Number of grid points along axes.

  • pbc – Periodic boundary conditions flag(s).

  • comm (MPIComm) – Communicator for domain-decomposition.

  • kpt (Vector | None) – K-point for Block-boundary conditions specified in units of the reciprocal cell.

  • decomp (Sequence[Sequence[int]] | None) – Decomposition of the domain.

  • dtype – Data-type (float or complex).

atom_centered_functions(functions, positions, *, qspiral_v=None, atomdist=None, integral=None, cut=False, xp=None)[source]

Create UGAtomCenteredFunctions object.

blocks(data)[source]

Yield views of blocks of data.

ecut_max()[source]
eikr(kpt_c=None)[source]

Plane wave.

\[e^{i\mathbf{k}\cdot \mathbf{r}}\]
Parameters:

kpt_c (Optional[Union[Sequence[float], ndarray]]) – k-point in units of the reciprocal cell. Defaults to the UGDesc objects own k-point.

empty(dims=(), comm=MPIComm(size=1, rank=0), xp=<module 'numpy' from '/tmp/gpaw-docs/venv/lib/python3.10/site-packages/numpy/__init__.py'>)[source]

Create new UGArray object.

Parameters:
  • dims (int | tuple[int, ...]) – Extra dimensions.

  • comm (_Communicator) – Distribute dimensions along this communicator.

fft_plans(flags=0, xp=<module 'numpy' from '/tmp/gpaw-docs/venv/lib/python3.10/site-packages/numpy/__init__.py'>, dtype=None)[source]

Create FFTW-plans.

classmethod from_cell_and_grid_spacing(cell, grid_spacing, pbc=(True, True, True), kpt=None, comm=MPIComm(size=1, rank=0), dtype=None)[source]

Create UGDesc from grid-spacing.

from_data(data)[source]
global_shape()[source]

Actual size of uniform grid.

new(*, kpt=None, dtype=None, comm='inherit', size=None, pbc=None, decomp=None)[source]

Create new uniform grid description.

property phase_factor_cd

Phase factor for block-boundary conditions.

ranks_from_fractional_positions(fracpos_ac)[source]
property size

Size of uniform grid.

transformer(other, stencil_range=3, xp=<module 'numpy' from '/tmp/gpaw-docs/venv/lib/python3.10/site-packages/numpy/__init__.py'>)[source]

Create transformer from one grid to another.

(for interpolation and restriction).

xyz()[source]

Create array of (x, y, z) coordinates.

class gpaw.core.PWDesc(*, ecut, cell, kpt=None, comm=MPIComm(size=1, rank=0), dtype=None)[source]

Description of plane-wave basis.

Parameters:
  • ecut (float) – Cutoff energy for kinetic energy of plane waves (units: hartree).

  • cell (ArrayLike1D | ArrayLike2D) – Unit cell given as three floats (orthorhombic grid), six floats (three lengths and the angles in degrees) or a 3x3 matrix (units: bohr).

  • comm (MPIComm) – Communicator for distribution of plane-waves.

  • kpt (Vector | None) – K-point for Block-boundary conditions specified in units of the reciprocal cell.

  • dtype – Data-type (float or complex).

atom_centered_functions(functions, positions, *, qspiral_v=None, atomdist=None, integral=None, cut=False, xp=None)[source]

Create PlaneWaveAtomCenteredFunctions object.

cut(array_Q)[source]

Cut out G-vectors with (G+k)^2/2<E_kin.

empty(dims=(), comm=MPIComm(size=1, rank=0), xp=None)[source]

Create new PlaneWaveExpanions object.

Parameters:
  • dims (int | tuple[int, ...]) – Extra dimensions.

  • comm (_Communicator) – Distribute dimensions along this communicator.

from_data(data)[source]
global_shape()[source]

Tuple with one element: number of plane waves.

indices(shape)[source]

Return indices into FFT-grid.

itemsize: int = 16
kinetic_energies()[source]

Kinetic energy of plane waves.

\[|\mathbf{G}+\mathbf{k}|^{2} / 2\]
map_indices(other)[source]

Map from one (distributed) set of plane waves to smaller global set.

Say we have 9 G-vector on two cores:

5 3 4             . 3 4           0 . .
2 0 1 -> rank=0:  2 0 1  rank=1:  . . .
8 6 7             . . .           3 1 2

and we want a mapping to these 5 G-vectors:

  3
2 0 1
  4

On rank=0: the return values are:

[0, 1, 2, 3], [[0, 1, 2, 3], [4]]

and for rank=1:

[1], [[0, 1, 2, 3], [4]]
new(*, ecut=None, kpt=None, dtype=None, comm='inherit')[source]

Create new plane-wave expansion description.

paste(coef_G, array_Q)[source]

Paste G-vectors with (G+k)^2/2<E_kin into 3-D FFT grid and zero-pad.

reciprocal_vectors()[source]

Returns reciprocal lattice vectors, G + k, in xyz coordinates.

class gpaw.core.atom_centered_functions.AtomCenteredFunctions(functions, fracpos_ac, atomdist=None, xp=None)[source]
add_to(functions, coefs=1.0)[source]

Add atom-centered functions multiplied by coefs to functions.

property atomdist
derivative(functions, out=None)[source]

Calculate derivatives of integrals with respect to atom positions.

empty(dims=(), comm=MPIComm(size=1, rank=0))[source]

Create AtomsArray for coefficients.

integrate(functions, out=None)[source]

Calculate integrals of atom-centered functions multiplied by functions.

property layout
move(fracpos_ac, atomdist)[source]

Move atoms to new positions.

new(desc, atomdist)[source]
stress_contribution(a, c=1.0)[source]
class gpaw.core.UGArray(grid, dims=(), comm=MPIComm(size=1, rank=0), data=None, xp=None)[source]

Object for storing function(s) on a uniform grid.

Parameters:
  • grid (UGDesc) – Description of uniform grid.

  • dims (int | tuple[int, ...]) – Extra dimensions.

  • comm (MPIComm) – Distribute dimensions along this communicator.

  • data (np.ndarray | None) – Data array for storage.

abs_square(weights, out=None)[source]

Add weighted absolute square of data to output array.

add_ked(occ_n, taut_R)[source]

Add weighted absolute square of gradient of data to output array.

fft(plan=None, pw=None, out=None)[source]

Do FFT.

\[\int d\mathbf{r} e^{i\mathbf{G}\cdot \mathbf{r}} f(\mathbf{r})\]
fft_restrict(plan1=None, plan2=None, grid=None, out=None)[source]

Restrict to coarser grid.

Parameters:
gather(out=None, broadcast=False)[source]

Gather data from all ranks to rank-0.

integrate(other=None, skip_sum=False)[source]

Integral of self or self times cc(other).

interpolate(plan1=None, plan2=None, grid=None, out=None)[source]

Interpolate to finer grid.

Parameters:
moment(center_v=None)[source]

Calculate moment of data.

multiply_by_eikr(kpt_c=None)[source]

Multiply by \(exp(ik.r)\).

new(data=None, zeroed=False)[source]

Create new UniforGridFunctions object of same kind.

Parameters:

data – Array to use for storage.

norm2()[source]

Calculate integral over cell of absolute value squared.

\[\int |a(\mathbf{r})|^{2} d\mathbf{r}\]
randomize(seed=None)[source]

Insert random numbers between -0.5 and 0.5 into data.

redist(domain, comm1, comm2)[source]
scaled(s, v=1.0)[source]

Create new scaled UGArray object.

Unit cell axes are multiplied by \(s\) and data by \(v\).

scatter_from(data=None)[source]

Scatter data from rank-0 to all ranks.

symmetrize(rotation_scc, translation_sc)[source]

Make data symmetric.

to_pbc_grid()[source]

Convert to UniformGrud with pbc=(True, True, True).

xy(*axes)[source]

Extract x, y values along line.

Useful for plotting:

x, y = grid.xy(0, ..., 0)
plt.plot(x, y)
class gpaw.core.arrays.DistributedArrays(dims, myshape, comm, domain_comm, data, dv, dtype, xp=None)[source]
abs_square(weights, out)[source]

Add weighted absolute square of data to output array.

See also XKCD: 849.

add_ked(weights, out)[source]

Add weighted absolute square of gradient of data to output array.

copy()[source]
desc: DomainType
flat()[source]
gather(out=None, broadcast=False)[source]
gathergather()[source]
interpolate(plan1=None, plan2=None, grid=None, out=None)[source]
property matrix: Matrix
matrix_elements(other, *, out=None, symmetric='_default', function=None, domain_sum=True, cc=False)[source]
new(data=None)[source]
redist(domain, comm1, comm2)[source]
scatter_from(data=None)[source]
to_xp(xp)[source]
class gpaw.core.atom_arrays.AtomArrays(layout, dims=(), comm=MPIComm(size=1, rank=0), data=None)[source]

AtomArrays object.

Parameters:
  • layout (AtomArraysLayout) – Layout-description.

  • dims (int | tuple[int, ...]) – Extra dimensions.

  • comm (MPIComm) – Distribute dimensions along this communicator.

  • data (np.ndarray | None) – Data array for storage.

block_diag_multiply(block_diag_matrix_axii, out_ani, index=None)[source]

Multiply by block diagonal matrix.

with A, B and C refering to self, block_diag_matrix_axii and out_ani:

\[\sum^{}_{i} A^{a}_{ni} B^{a}_{ij} \rightarrow C^{a}_{nj}\]

If index is not None, block_diag_matrix_axii must have an extra dimension: \(B_{ij}^{ax}\) and x=index is used.

gather(broadcast: Literal[False] = False, copy: bool = False) gpaw.core.atom_arrays.AtomArrays | None[source]
gather(broadcast: Literal[True], copy: bool = False) AtomArrays

Gather all atoms on master.

get(a)[source]
items()[source]
keys()[source]
property matrix: Matrix
moved(atomdist)[source]
new(*, layout=None, data=None, xp=None)[source]

Create new AtomArrays object of same kind.

Parameters:
  • layout – Layout-description.

  • data – Array to use for storage.

redist(atomdist, comm1, comm2)[source]
scatter_from(data=None)[source]
to_cpu()[source]
to_full()[source]

Convert \(N(N+1)/2\) vectors to \(N\times N\) matrices.

>>> a = AtomArraysLayout([6]).empty()
>>> a[0][:] = [1, 2, 3, 4, 5, 6]
>>> a.to_full()[0]
array([[1., 2., 4.],
       [2., 3., 5.],
       [4., 5., 6.]])
to_lower_triangle()[source]

Convert \(N*N\) matrices to \(N*(N+1)/2\) vectors.

>>> a = AtomArraysLayout([(3, 3)]).empty()
>>> a[0][:] = [[11, 12, 13],
...            [12, 22, 23],
...            [13, 23, 33]]
>>> a.to_lower_triangle()[0]
array([11., 12., 22., 13., 23., 33.])
to_xp(xp)[source]
values()[source]
class gpaw.core.atom_arrays.AtomArraysLayout(shapes, atomdist=MPIComm(size=1, rank=0), dtype=<class 'float'>, xp=None)[source]

Description of layout of atom arrays.

Parameters:
  • shapes (Sequence[int | tuple[int, ...]]) – Shapse of arrays - one for each atom.

  • atomdist (AtomDistribution | MPIComm) – Distribution of atoms.

  • dtype – Data-type (float or complex).

empty(dims=(), comm=MPIComm(size=1, rank=0))[source]

Create new AtomArrays object.

Parameters:
  • dims (int | tuple[int, ...]) – Extra dimensions.

  • comm (_Communicator) – Distribute dimensions along this communicator.

new(atomdist=None, dtype=None, xp=None)[source]

Create new AtomsArrayLayout object with new atomdist.

sizes()[source]

Compute array sizes for all ranks.

>>> AtomArraysLayout([3, 4]).sizes()
([{0: 3, 1: 4}], array([7]))
zeros(dims=(), comm=MPIComm(size=1, rank=0))[source]
class gpaw.core.atom_arrays.AtomDistribution(ranks, comm=MPIComm(size=1, rank=0))[source]

Atom-distribution.

Parameters:
  • ranks (ArrayLike1D) – List of ranks, one rank per atom.

  • comm (MPIComm) – MPI-communicator.

classmethod from_atom_indices(atom_indices, comm=MPIComm(size=1, rank=0), *, natoms=None)[source]

Create distribution from atom indices.

>>> AtomDistribution.from_atom_indices([0, 1, 2]).rank_a
array([0, 0, 0])
classmethod from_number_of_atoms(natoms, comm=MPIComm(size=1, rank=0))[source]

Distribute atoms evenly.

>>> AtomDistribution.from_number_of_atoms(3).rank_a
array([0, 0, 0])
gather()[source]
class gpaw.core.PWArray(pw, dims=(), comm=MPIComm(size=1, rank=0), data=None, xp=None)[source]

Object for storing function(s) as a plane-wave expansions.

Parameters:
  • pw (PWDesc) – Description of plane-waves.

  • dims (int | tuple[int, ...]) – Extra dimensions.

  • comm (MPIComm) – Distribute plane-waves along this communicator.

  • data (np.ndarray | None) – Data array for storage.

abs_square(weights, out, _slow=False)[source]

Add weighted absolute square of self to output array.

With \(a_n(G)\) being self and \(w_n\) the weights:

\[out(\mathbf{r}) \leftarrow out(\mathbf{r}) + \sum^{}_{n} |FFT^{-1} [a_{n} (\mathbf{G})]|^{2} w_{n}\]
add_ked(occ_n, taut_R)[source]

Add weighted absolute square of gradient of data to output array.

copy()[source]

Create a copy (surprise!).

gather(out=None, broadcast=False)[source]

Gather coefficients on master.

gather_all(out)[source]

Gather coefficients from self[r] on rank r.

On rank r, an array of all G-vector coefficients will be returned. These will be gathered from self[r] on all the cores.

ifft(*, plan=None, grid=None, out=None, periodic=False)[source]

Do inverse FFT(s) to uniform grid(s).

Parameters:
  • plan – Plan for inverse FFT.

  • grid – Target grid.

  • out – Target UGArray object.

integrate(other=None)[source]

Integral of self or self time cc(other).

interpolate(plan1=None, plan2=None, grid=None, out=None)[source]
property matrix: Matrix

Matrix view of data.

moment()[source]
morph(pw)[source]
new(data=None)[source]

Create new PWArray object of same kind.

Parameters:

data – Array to use for storage.

norm2(kind='normal')[source]

Calculate integral over cell.

For kind=’normal’ we calculate:

\[\int |a(\mathbf{r})|^{2} d\mathbf{r} = \sum^{}_{G} |c_{G} |^{2} V,\]

where V is the volume of the unit cell.

And for kind=’kinetic’:

\[\frac{1}{2} \sum^{}_{G} |c_{G} |^{2} G^{2} V,\]
randomize(seed=None)[source]

Insert random numbers between -0.5 and 0.5 into data.

scatter_from(data=None)[source]

Scatter data from rank-0 to all ranks.

scatter_from_all(a_G)[source]

Scatter all coefficients from rank r to self on other cores.

to_pbc_grid()[source]
class gpaw.core.plane_waves.Empty(dims)[source]
class gpaw.core.matrix.Matrix(M, N, dtype=None, data=None, dist=None, xp=None)[source]

Matrix object.

Parameters:
  • M (int) – Rows.

  • N (int) – Columns.

  • dtype – Data type (float or complex).

  • dist (MatrixDistribution | tuple | None) – BLACS distribution given as (communicator, rows, columns, blocksize) tuple. Default is None meaning no distribution.

  • data (ArrayLike2D | None) – Numpy ndarray to use for storage. By default, a new ndarray will be allocated.

add_hermitian_conjugate(scale=1.0)[source]

Add hermitian conjugate to myself.

add_to_diagonal(d)[source]

Add list of numbers or single number to diagonal of matrix.

complex_conjugate()[source]

Inplace complex conjugation.

copy()[source]

Create a copy.

eigh(S=None, *, cc=False, scalapack=(None, 1, 1, None), limit=None)[source]

Calculate eigenvectors and eigenvalues.

Matrix must be symmetric/hermitian and stored in lower half. If S is given, solve a generalized eigenvalue problem.

Parameters:
  • cc (bool) – Complex conjugate matrix before finding eigenvalues.

  • scalapack (tuple) – BLACS distribution for ScaLapack to use. Default is to do serial diagonalization.

  • limit (Optional[int]) – Number of eigenvector and values to find. Defaults to all.

eighg(L, comm2=MPIComm(size=1, rank=0))[source]

Solve generalized eigenvalue problem.

With \(H\) being self, we solve for the eigenvectors \(C\) and the eigenvalues \(Λ\) (a diagonal matrix):

\[HC = SCΛ,\]

where \(L\) is a lower triangle matrix such that:

\[LSL^{†} = 1\cdot\]

The solution has these three steps:

\[\tilde{H} = LHL^{†} , \tilde{H}\tilde{C} = \tilde{C}Λ, C = L^{†} \tilde{C}\cdot\]

Note that \(H\) must be the full matrix not just half of it!

gather(root=0, broadcast=False)[source]

Gather the Matrix on the root rank.

Returns a new Matrix distributed so that all data is on the root rank

inv(uplo='L')[source]

Inplace inversion.

invcholesky()[source]

In-place inverse of Cholesky decomposition.

Calculate a lower triangle matrix \(L\) where:

\[LSL^{†} = 1,\]

and \(S\) is self. Only the lower part of \(S\) is used.

>>> S = Matrix(2, 2, data=[[1.0, np.nan],
...                        [0.1, 1.0]])
>>> S.invcholesky()
>>> S.data
array([[ 1.        , -0.        ],
       [-0.10050378,  1.00503782]])
multiply(other, alpha=1.0, opa='N', opb='N', out=None, beta=0.0, symmetric=False)[source]

BLAS matrix-multiplication with other matrix.

new(dist='inherit', data=None)[source]

Create new matrix of same shape and dtype.

Default is to use same BLACS distribution. Use dist to use another distribution.

redist(other)[source]

Redistribute to other BLACS layout.

to_cpu()[source]
tril2full()[source]

Fill in upper triangle from lower triangle.

For a real matrix:

a ? ?    a b d
b c ? -> b c e
d e f    d e f

For a complex matrix, the complex conjugate of the lower part will be inserted into the upper part.

class gpaw.core.matrix.MatrixDistribution[source]
blocksize: int | None
columns: int
comm: _Communicator
desc: ndarray
eighg(H, L)[source]
full_shape: tuple[int, int]
matrix(dtype=None, data=None)[source]
multiply(alpha, a, opa, b, opb, beta, c, symmetric)[source]
my_row_range()[source]

Return indices for range of my rows.

>>> Matrix(2, 2).dist.my_row_range()
(0, 2)
new(M, N)[source]
rows: int
shape: tuple[int, int]

DFT

class gpaw.new.calculation.DFTCalculation(state, setups, scf_loop, pot_calc, log)[source]
converge(convergence=None, maxiter=None, steps=99999999999999999, calculate_forces=None)[source]

Converge to self-consistent solution of Kohn-Sham equation.

densities()[source]
dipole()[source]
electrostatic_potential()[source]
energies()[source]
forces(silent=False)[source]

Calculate atomic forces.

classmethod from_parameters(atoms, params, comm=None, log=None, builder=None)[source]

Create DFTCalculation object from parameters and atoms.

iconverge(convergence=None, maxiter=None, calculate_forces=None)[source]
magmoms()[source]
move_atoms(atoms)[source]
new(atoms, params, log=None)[source]

Create new DFTCalculation object.

stress()[source]
write_converged()[source]
class gpaw.new.calculation.DFTState(ibzwfs, density, potential)[source]

State of a Kohn-Sham calculation.

move(fracpos_ac, atomdist)[source]
class gpaw.new.density.Density(nt_sR, taut_sR, D_asii, charge, delta_aiiL, delta0_a, N0_aii, l_aj, nct_aX, tauct_aX)[source]
calculate_compensation_charge_coefficients()[source]
calculate_dipole_moment(fracpos_ac)[source]
calculate_magnetic_moments()[source]
classmethod from_data_and_setups(nt_sR, taut_sR, D_asii, charge, setups, nct_aX, tauct_aX)[source]
classmethod from_superposition(*, grid, nct_aX, tauct_aX, atomdist, setups, basis_set, magmom_av, ncomponents, charge=0.0, hund=False, mgga=False)[source]
move(fracpos_ac, atomdist)[source]
property nct_R
new(new_grid, pw, fracpos_ac, atomdist)[source]
normalize()[source]
redist(grid, xdesc, atomdist, comm1, comm2)[source]
symmetrize(symmetries)[source]
property tauct_R
update(ibzwfs, ked=False)[source]
update_ked(ibzwfs, symmetrize=True)[source]
write(writer)[source]
gpaw.new.builder.builder(atoms, params, comm=None)[source]

Create DFT-components builder.

  • pw

  • lcao

  • fd

  • tb

  • atom

class gpaw.new.ibzwfs.IBZWaveFunctions(ibz, *, nelectrons, ncomponents, wfs_qs, kpt_comm=MPIComm(size=1, rank=0), kpt_band_comm=MPIComm(size=1, rank=0), comm=MPIComm(size=1, rank=0))[source]

Collection of wave function objects for k-points in the IBZ.

add_to_density(nt_sR, D_asii)[source]

Compute density and add to nt_sR and D_asii.

add_to_ked(taut_sR)[source]
calculate_occs(occ_calc, fixed_fermi_level=False)[source]
forces(potential)[source]
get_all_eigs_and_occs()[source]
get_all_electron_wave_function(band, kpt=0, spin=0, grid_spacing=0.05, skip_paw_correction=False)[source]
get_eigs_and_occs(k=0, s=0)[source]
get_homo_lumo(spin=None)[source]

Return HOMO and LUMO eigenvalues.

get_max_shape(global_shape=False)[source]

Find the largest wave function array shape.

For a PW-calculation, this shape could depend on k-point.

get_wfs(*, kpt=0, spin=0, n1=0, n2=0)[source]
make_sure_wfs_are_read_from_gpw_file()[source]
property mode
move(fracpos_ac, atomdist)[source]
orthonormalize(work_array_nX=None)[source]
write(writer, skip_wfs)[source]

Write fermi-level(s), eigenvalues, occupation numbers, …

… k-points, symmetry information, projections and possibly also the wave functions.

write_summary(log)[source]
class gpaw.new.potential.Potential(vt_sR, dH_asii, dedtaut_sR, energies, vHt_x=None)[source]
dH(P_ani, out_ani, spin)[source]
move(atomdist)[source]

Move atoms inplace.

redist(grid, desc, atomdist, comm1, comm2)[source]
class gpaw.new.pot_calc.PotentialCalculator(xc, poisson_solver, setups, *, fracpos_ac, external_potential=None, soc=False)[source]
calculate(density, ibzwfs=None, vHt_x=None, kpt_band_comm=None)[source]
calculate_charges(vHt_x)[source]
calculate_pseudo_potential(density, ibzwfs, vHt_x)[source]
restrict(a_r, a_R=None)[source]
class gpaw.new.scf.SCFLoop(hamiltonian, occ_calc, eigensolver, mixer, comm, convergence, maxiter)[source]
iterate(state, pot_calc, convergence=None, maxiter=None, calculate_forces=None, log=None)[source]
class gpaw.new.input_parameters.InputParameters(params, warn=True)[source]
basis: Any
charge: float
convergence: dict[str, Any]
eigensolver: dict[str, Any]
experimental: dict[str, Any]
external: dict[str, Any]
gpts: Union[None, Sequence[int]]
h: float | None
hund: bool
items()[source]
kpts: dict[str, Any]
magmoms: Any
mode: dict[str, Any]
nbands: None | int | str
parallel: dict[str, Any]
poissonsolver: dict[str, Any]
setups: Any
soc: bool
spinpol: bool
symmetry: dict[str, Any]
xc: dict[str, Any]
class gpaw.new.pwfd.wave_functions.PWFDWaveFunctions(psit_nX, *, spin, q, k, setups, fracpos_ac, atomdist, weight=1.0, ncomponents=1, qspiral_v=None)[source]
property P_ani: AtomArrays

PAW projections.

\[\langle \tilde{p}^{a}_{i}|\tilde{ψ}_{n} \rangle\]
add_to_density(nt_sR, D_asii)[source]
add_to_ked(taut_sR)[source]
array_shape(global_shape=False)[source]
collect(n1=0, n2=0)[source]

Collect range of bands to master of band and domain communicators.

dipole_matrix_elements(center_v=None)[source]

Calculate dipole matrix-elements.

\[\mathbf{μ}_{mn} = \int d\mathbf{r} \tilde{𝜓}_{m} \tilde{𝜓}_{n} \mathbf{r} + \sum^{}_{aij} P^{a}_{im} P^{a}_{jn} Δ\mathbf{μ}^{a}_{ij}\]
Parameters:

center_v (Optional[Union[Sequence[float], ndarray]]) – Center of molecule. Defaults to center of cell.

Returns:

matrix elements in atomic units.

Return type:

Array3D

force_contribution(potential, F_av)[source]
classmethod from_wfs(wfs, psit_nX, fracpos_ac=None, atomdist=None)[source]
gather_wave_function_coefficients()[source]
morph(desc, fracpos_ac, atomdist)[source]
move(fracpos_ac, atomdist)[source]
orthonormalize(work_array_nX=None)[source]

Orthonormalize wave functions.

Computes the overlap matrix:

\[S_{mn} = \int \tilde{ψ}_{m}(\mathbf{r})^{*} \tilde{ψ}_{n}(\mathbf{r}) d\mathbf{r} + \sum^{}_{aij} (P^{a}_{im} )^{*} P^{a}_{jn} ΔS^{a}_{ij}\]

With \(LSL^\dagger=1\), we update the wave functions and projections inplace like this:

\[Ψ_{m} \leftarrow \sum^{}_{n} L^{*}_{mn} Ψ_{n} ,\]

and:

\[P^{a}_{mi} \leftarrow \sum^{}_{n} L^{*}_{mn} P^{a}_{ni} \cdot\]
property pt_aiX: AtomCenteredFunctions

PAW projector functions.

\[\tilde{p}^{a}_{i} (\mathbf{r})\]
receive(rank, comm)[source]
send(rank, comm)[source]
subspace_diagonalize(Ht, dH, work_array=None, Htpsit_nX=None, scalapack_parameters=(None, 1, 1, None))[source]

Ht(in, out):

\[\tilde{H} = \hat{T} + \tilde{v}\]

dH:

\[\langle \tilde{𝜓}|\tilde{p} \rangle ΔH^{a}_{ij} \langle \tilde{p}_{j}|\tilde{𝜓} \rangle\]
to_uniform_grid_wave_functions(grid, basis)[source]
class gpaw.new.ase_interface.ASECalculator(params, *, log, calculation=None, atoms=None)[source]

This is the ASE-calculator frontend for doing a GPAW calculation.

band_structure()[source]

Create band-structure object for plotting.

calculate(atoms, properties=None, system_changes=None)[source]
calculate_property(atoms, prop)[source]

Calculate (if not already calculated) a property.

The prop string must be one of

  • energy

  • forces

  • stress

  • magmom

  • magmoms

  • dipole

converge()[source]

Iterate to self-consistent solution.

Will also calculate “cheap” properties: energy, magnetic moments and dipole moment.

converge_wave_functions()[source]
create_new_calculation(atoms)[source]
create_new_calculation_from_old(atoms)[source]
property density
diagonalize_full_hamiltonian(nbands=None, scalapack=None, expert=None)[source]
dos(soc=False, theta=0.0, phi=0.0, shift_fermi_level=True)[source]

Create DOS-calculator.

Default is to shift_fermi_level to 0.0 eV. For soc=True, angles can be given in degrees.

fixed_density(txt='-', **kwargs)[source]
get_all_electron_density(spin=None, gridrefinement=1, broadcast=True, skip_core=False)[source]
get_atomic_electrostatic_potentials()[source]
get_atoms()[source]
get_bz_k_points()[source]
get_dipole_moment(atoms=None)[source]
get_effective_potential(spin=0)[source]
get_eigenvalues(kpt=0, spin=0, broadcast=True)[source]
get_electrostatic_corrections()[source]
get_electrostatic_potential()[source]
get_fermi_level()[source]
get_fermi_levels()[source]
get_forces(atoms=None)[source]
get_homo_lumo(spin=None)[source]
get_ibz_k_points()[source]
get_magnetic_moment(atoms=None)[source]
get_magnetic_moments(atoms=None)[source]
get_non_collinear_magnetic_moment(atoms=None)[source]
get_non_collinear_magnetic_moments(atoms=None)[source]
get_number_of_bands()[source]
get_number_of_electrons()[source]
get_number_of_grid_points()[source]
get_number_of_iterations()[source]
get_number_of_spins()[source]
get_occupation_numbers(kpt=0, spin=0, broadcast=True)[source]
get_orbital_magnetic_moments()[source]

Return the orbital magnetic moment vector for each atom.

get_potential_energy(atoms=None, force_consistent=False)[source]
get_property(name, atoms=None, allow_calculation=True)[source]
get_pseudo_density(spin=None, gridrefinement=1, broadcast=True)[source]
get_pseudo_wave_function(band, kpt=0, spin=None, periodic=False, broadcast=True)[source]
get_reference_energy()[source]
get_stress(atoms=None)[source]
get_xc_difference(xcparams)[source]

Calculate non-selfconsistent XC-energy difference.

get_xc_functional()[source]
gs_adapter()[source]
property hamiltonian
implemented_properties = ['energy', 'free_energy', 'forces', 'stress', 'dipole', 'magmom', 'magmoms']
initialize(atoms)[source]
property initialized
move_atoms(atoms)[source]
name = 'gpaw'
new(**kwargs)[source]
property parameters
property results
property setups
property spos_ac
property symmetry
property wfs
property world
write(filename, mode='')[source]

Write calculator object to a file.

Parameters:
  • filename – File to be written

  • mode – Write mode. Use mode='all' to include wave functions in the file.

gpaw.new.ase_interface.GPAW(filename=None, txt='?', communicator=None, **kwargs)[source]

Create ASE-compatible GPAW calculator.

FFTW

Python wrapper for FFTW3 library

class gpaw.fftw.FFTPlans(size_c, dtype, empty=<function empty>)[source]
class gpaw.fftw.FFTPlan(in_R, out_R, sign, flags=0)[source]

FFT 3d transform.

class gpaw.fftw.FFTWPlan(in_R, out_R, sign, flags=0)[source]

FFTW3 3d transform.

class gpaw.fftw.FFTWPlans(size_c, dtype, flags=0)[source]

FFTW3 3d transforms.

fft()[source]

Do FFT from tmp_R to tmp_Q.

>>> plans = create_plans([4, 1, 1], float)
>>> plans.tmp_R[:, 0, 0] = [1, 0, 1, 0]
>>> plans.fft()
>>> plans.tmp_Q[:, 0, 0]
array([2.+0.j, 0.+0.j, 2.+0.j, 0.+0.j])
ifft()[source]

Do inverse FFT from tmp_Q to tmp_R.

>>> plans = create_plans([4, 1, 1], complex)
>>> plans.tmp_Q[:, 0, 0] = [0, 1j, 0, 0]
>>> plans.ifft()
>>> plans.tmp_R[:, 0, 0]
array([ 0.+1.j, -1.+0.j,  0.-1.j,  1.+0.j])
class gpaw.fftw.NumpyFFTPlan(in_R, out_R, sign, flags=0)[source]

Numpy fallback.

class gpaw.fftw.NumpyFFTPlans(size_c, dtype, empty=<function empty>)[source]

Numpy fallback.

fft()[source]

Do FFT from tmp_R to tmp_Q.

>>> plans = create_plans([4, 1, 1], float)
>>> plans.tmp_R[:, 0, 0] = [1, 0, 1, 0]
>>> plans.fft()
>>> plans.tmp_Q[:, 0, 0]
array([2.+0.j, 0.+0.j, 2.+0.j, 0.+0.j])
ifft()[source]

Do inverse FFT from tmp_Q to tmp_R.

>>> plans = create_plans([4, 1, 1], complex)
>>> plans.tmp_Q[:, 0, 0] = [0, 1j, 0, 0]
>>> plans.ifft()
>>> plans.tmp_R[:, 0, 0]
array([ 0.+1.j, -1.+0.j,  0.-1.j,  1.+0.j])
gpaw.fftw.check_fft_size(n, factors=[2, 3, 5, 7])[source]

Check if n is an efficient fft size.

Efficient means that n can be factored into small primes (2, 3, 5, 7).

>>> check_fft_size(17)
False
>>> check_fft_size(18)
True
gpaw.fftw.create_plans(size_c, dtype, flags=0, xp=<module 'numpy' from '/tmp/gpaw-docs/venv/lib/python3.10/site-packages/numpy/__init__.py'>)[source]

Create plan-objects for FFT and inverse FFT.

gpaw.fftw.empty(shape, dtype=<class 'float'>)[source]

numpy.empty() equivalent with 16 byte alignment.

gpaw.fftw.get_efficient_fft_size(N, n=1, factors=[2, 3, 5, 7])[source]

Return smallest efficient fft size.

Must be greater than or equal to N and divisible by n.

>>> get_efficient_fft_size(17)
18
gpaw.fftw.have_fftw()[source]

Did we compile with FFTW?

BLAS

gpaw.utilities.blas.mmm(alpha, a, opa, b, opb, beta, c)[source]

Matrix-matrix multiplication using dgemm or zgemm.

For opa=’N’ and opb=’N’, we have:

\[c \leftarrow αab + βc\cdot\]

Use ‘T’ to transpose matrices and ‘C’ to transpose and complex conjugate matrices.

gpaw.utilities.blas.rk(alpha, a, beta, c, trans='c')[source]

Rank-k update of a matrix.

For trans='c' the following operation is performed:

\[c \leftarrow αaa^{†} + βc,\]

and for trans='t' we get:

\[c \leftarrow αa^{†} a + βc\]

If the a array has more than 2 dimensions then the 2., 3., … axes are combined.

Only the lower triangle of c will contain sensible numbers.