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), 'txt': None}
>>> 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)

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'}, txt: 'h2.txt')

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

Do-it-yourself example

Let’s try to build a DFT calculation without any of the shortcuts (gpaw.new.builder.builder(), gpaw.new.ase_interface.GPAW() or gpaw.new.calculation.DFTCalculation.from_parameters()).

So, instead of simple three-line example above where we calculate the energy of an H2 molecule, we do it the hard way:

from ase import Atoms
from ase.units import Bohr
from gpaw.setup import Setups
from gpaw.new.xc import XCFunctional
from gpaw.new.density import Density
from gpaw.core import UniformGrid
from gpaw.new.basis import create_basis
from gpaw.new.symmetry import create_symmetries_object
from gpaw.new.brillouin import BZPoints
from gpaw.core.atom_arrays import AtomDistribution
from gpaw.mpi import world

atoms = Atoms('H2',
              positions=[(0, 0, 0), (0, 0, 0.75)],
              cell=[2, 2, 3],
              pbc=True)

xc = XCFunctional('LDA')
setups = Setups([1, 1], {}, {}, xc)
grid = UniformGrid(cell=atoms.cell / Bohr,
                   size=[10, 10, 15],
                   pbc=[True, True, True])
symmetries = create_symmetries_object(atoms)
ibz = symmetries.reduce(BZPoints([[0, 0, 0]]))
basis = create_basis(
    ibz,
    nspins=1,
    pbc_c=atoms.pbc,
    grid=grid,
    setups=setups,
    dtype=float,
    fracpos_ac=atoms.get_scaled_positions())
nct_R = grid.zeros()
atomdist = AtomDistribution([0, 0], world)
density = Density.from_superposition(grid, nct_R, atomdist, setups, basis)
pot_calc = ...
potential = ...
wfs = ...
ibzwfs = ...
state = ...
scf_loop = ...
for _ in scf_loop.iterate(state, pot_calc):
    ...

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})\)

UniformGridFunctions

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})\)

UniformGridFunctions | PlaneWaveExpansions

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 UniformGrid class:

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

Given a UniformGrid object, one can create UniformGridFunctions 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 UniformGrid class:

size()

Size of uniform grid.

global_shape()

Actual size of uniform grid.

phase_factors_cd()

new()

Create new uniform grid description.

empty()

Create new UniformGridFunctions object.

blocks()

Yield views of blocks of data.

xyz()

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

atom_centered_functions()

Create UniformGridAtomCenteredFunctions object.

transformer()

Create transformer from one grid to another.

eikr()

Plane wave.

from_cell_and_grid_spacing()

classmethod(function) -> method

fft_plans()

Create FFTW-plans.

ranks_from_fractional_positions()

and the UniformGridFunctions 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 time 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.

Plane waves

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

>>> from gpaw.core import PlaneWaves
>>> pw = PlaneWaves(ecut=100, cell=grid.cell)
>>> func_G = pw.empty()
>>> func_R.fft(out=func_G)
PlaneWaveExpansions(pw=PlaneWaves(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)
UniformGridFunctions(grid=UniformGrid(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 PlaneWaves class:

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.

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 PlaneWaveExpansions class:

new()

Create new PlaneWaveExpansions object of same kind.

copy()

Create a copy (surprise!).

matrix()

Matrix view of data.

ifft()

Do inverse FFT to uniform grid.

interpolate()

Do inverse FFT to uniform grid.

gather()

Gather coefficients on master.

scatter_from()

Scatter data from rank-0 to all ranks.

integrate()

Integral of self or self time cc(other).

norm2()

Calculate integral over cell.

abs_square()

Add weighted absolute square of data to output array.

to_pbc_grid()

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 UniformGrid


alpha = 4.0
rcut = 2.0
l = 0
gauss = (l, rcut, lambda r: (4 * np.pi)**0.5 * np.exp(-alpha * r**2))
grid = UniformGrid(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()

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.

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.UniformGrid(*, cell, size, pbc=(True, True, True), kpt=None, comm=<gpaw.mpi.SerialCommunicator object>, decomp=None, dtype=None)[source]

Description of 3-D 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.

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

  • pbc – Periodic boundary conditions flag(s).

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

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

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

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

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

Create UniformGridAtomCenteredFunctions object.

blocks(data)[source]

Yield views of blocks of data.

eikr(kpt_c=None)[source]

Plane wave.

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

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

empty(dims=(), comm=<gpaw.mpi.SerialCommunicator object>)[source]

Create new UniformGridFunctions object.

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

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

fft_plans(flags=0)[source]

Create FFTW-plans.

classmethod from_cell_and_grid_spacing(cell, grid_spacing, pbc=(True, True, True), kpt=None, comm=<gpaw.mpi.SerialCommunicator object>, dtype=None)[source]

Create UniformGrid from grid-spacing.

global_shape()[source]

Actual size of uniform grid.

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

Create new uniform grid description.

property phase_factors_cd
ranks_from_fractional_positions(fracpos_ac)[source]
property size

Size of uniform grid.

transformer(other, stencil_range=3)[source]

Create transformer from one grid to another.

(for interpolation and restriction).

xyz()[source]

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

class gpaw.core.PlaneWaves(*, ecut, cell, kpt=None, comm=<gpaw.mpi.SerialCommunicator object>, dtype=None)[source]

Description of plane-wave basis.

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

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

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

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

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

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

Create PlaneWaveAtomCenteredFunctions object.

cut(array_Q)[source]

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

empty(dims=(), comm=<gpaw.mpi.SerialCommunicator object>)[source]

Create new PlaneWaveExpanions object.

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

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

global_shape()[source]

Tuple with one element: number of plane waves.

indices(shape)[source]

Return indices into FFT-grid.

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, 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)[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=<gpaw.mpi.SerialCommunicator object>)[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.

class gpaw.core.uniform_grid.UniformGridFunctions(grid, dims=(), comm=<gpaw.mpi.SerialCommunicator object>, data=None)[source]

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

Parameters
  • grid – Description of uniform grid.

  • dims – Extra dimensions.

  • comm – Distribute dimensions along this communicator.

  • data – Data array for storage.

abs_square(weights, out=None)[source]

Add weighted absolute square of data to output array.

desc: DomainType
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)[source]

Integral of self or self time 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)[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()[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.

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)[source]
abs_square(weights, out=None)[source]

Add weighted absolute square of data to output array.

See also XKCD:849.

desc: DomainType
flat()[source]
gather(out=None, broadcast=False)[source]
property matrix: gpaw.core.matrix.Matrix
matrix_elements(other, *, out=None, symmetric=None, function=None, domain_sum=True, cc=False)[source]
new(data=None)[source]
class gpaw.core.atom_arrays.AtomArrays(layout, dims=(), comm=<gpaw.mpi.SerialCommunicator object>, 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) – Data array for storage.

gather(broadcast=False, copy=False)[source]

Gather all atoms on master.

get(a)[source]
items()[source]
keys()[source]
property matrix: gpaw.core.matrix.Matrix
new(layout=None, data=None)[source]

Create new AtomArrays object of same kind.

Parameters
  • layout – Layout-description.

  • data – Array to use for storage.

scatter_from(data=None)[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.])
values()[source]
class gpaw.core.atom_arrays.AtomArraysLayout(shapes, atomdist=<gpaw.mpi.SerialCommunicator object>, dtype=<class 'float'>)[source]

Description of layour of atom arrays.

Parameters
  • shapes (list[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=<gpaw.mpi.SerialCommunicator object>)[source]

Create new AtomArrays object.

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

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

new(atomdist=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]))
class gpaw.core.atom_arrays.AtomDistribution(ranks, comm=<gpaw.mpi.SerialCommunicator object>)[source]

Atom-distribution.

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

  • comm (MPIComm) – MPI-communicator.

classmethod from_atom_indices(atom_indices, comm=<gpaw.mpi.SerialCommunicator object>, *, 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=<gpaw.mpi.SerialCommunicator object>)[source]

Distribute atoms evenly.

>>> AtomDistribution.from_number_of_atoms(3).rank_a
array([0, 0, 0])
gather()[source]
class gpaw.core.plane_waves.PlaneWaveExpansions(pw, dims=(), comm=<gpaw.mpi.SerialCommunicator object>, data=None)[source]

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

Parameters
  • pw – Description of plane-waves.

  • dims – Extra dimensions.

  • comm – Distribute plane-waves along this communicator.

  • data – Data array for storage.

abs_square(weights, out=None)[source]

Add weighted absolute square of data to output array.

copy()[source]

Create a copy (surprise!).

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

Gather coefficients on master.

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

Do inverse FFT to uniform grid.

Parameters
  • plan – Plan for inverse FFT.

  • grid – Target grid.

  • out – Target UniformGridFunctions object.

integrate(other=None)[source]

Integral of self or self time cc(other).

interpolate(*, plan=None, grid=None, out=None)

Do inverse FFT to uniform grid.

Parameters
  • plan – Plan for inverse FFT.

  • grid – Target grid.

  • out – Target UniformGridFunctions object.

property matrix: gpaw.core.matrix.Matrix

Matrix view of data.

new(data=None)[source]

Create new PlaneWaveExpansions 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,\]
scatter_from(data=None)[source]

Scatter data from rank-0 to all ranks.

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)[source]

Matrix object.

Parameters
  • M (int) – Rows.

  • N (int) – Columns.

  • dtype – Data type (float or complex).

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

  • data (ArrayLike2D) – 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.

complex_conjugate()[source]

Inplace complex conjugation.

copy()[source]

Create a copy.

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

Calculate eigenvectors and eigenvalues.

Matrix must be symmetric/hermitian and stored in lower half.

cc: bool

Complex conjugate matrix before finding eigenvalues.

scalapack: tuple

BLACS distribution for ScaLapack to use. Default is to do serial diagonalization.

eighg(L, comm2=<gpaw.mpi.SerialCommunicator object>)[source]

Solve generalized eigenvalue problem.

\[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\]
gather(root=0)[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]

Inverse of Cholesky decomposition.

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

\[LSL^{†} = 1\]

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.

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: MPIComm
desc: Array1D
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)
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)[source]

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

dipole()[source]
energies()[source]
forces()[source]

Return atomic force contributions.

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

Create DFTCalculation object from parameters and atoms.

iconverge(convergence=None, maxiter=None)[source]
magmoms()[source]
move_atoms(atoms)[source]
stress()[source]
write_converged()[source]
class gpaw.new.calculation.DFTState(ibzwfs, density, potential, vHt_x=None, nct_R=None)[source]

State of a Kohn-Sham calculation.

move(fracpos_ac, atomdist, delta_nct_R)[source]
class gpaw.new.density.Density(nt_sR, D_asii, charge, delta_aiiL, delta0_a, N0_aii, l_aj)[source]
calculate_compensation_charge_coefficients()[source]
calculate_dipole_moment(fracpos_ac)[source]
calculate_magnetic_moments()[source]
classmethod from_data_and_setups(nt_sR, D_asii, charge, setups)[source]
classmethod from_superposition(grid, nct_R, atomdist, setups, basis_set, magmoms=None, charge=0.0, hund=False)[source]
move(delta_nct_R)[source]
normalize()[source]
symmetrize(symmetries)[source]
update(nct_R, ibzwfs)[source]
write(writer)[source]
xxxxx_overlap_correction(P_ani, out)[source]
gpaw.new.builder.builder(atoms, params)[source]

Create DFT-components builder.

  • pw

  • lcao

  • fd

  • tb

  • atom

class gpaw.new.ibzwfs.IBZWaveFunctions(ibz, nelectrons, ncomponents, wfs_qs, kpt_comm=<gpaw.mpi.SerialCommunicator object>)[source]

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

add_to_density(nt_sR, D_asii)[source]

Compute density from wave functions and add to nt_sR and D_asii.

calculate_occs(occ_calc, fixed_fermi_level=False)[source]
forces(dH_asii)[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]
is_master()[source]
make_sure_wfs_are_read_from_gpw_file()[source]
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, energies)[source]
dH(P_ani, out_ani, spin)[source]
write(writer)[source]
class gpaw.new.pot_calc.PotentialCalculator(xc, poisson_solver, setups, nct_R, fracpos_ac)[source]
calculate(density, vHt_x=None)[source]
move(fracpos_ac, atomdist, ndensities)[source]
class gpaw.new.scf.SCFLoop(hamiltonian, occ_calc, eigensolver, mixer, world, convergence, maxiter)[source]
iterate(state, pot_calc, convergence=None, maxiter=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]
force_complex_dtype: bool
gpts: None | Sequence[int]
h: float | None
items()[source]
kpts: dict[str, Any]
magmoms: Any
mode: dict[str, Any]
nbands: None | int | float
parallel: dict[str, Any]
poissonsolver: dict[str, Any]
setups: Any
soc: bool
spinpol: bool
symmetry: dict[str, Any]
txt: str | Path | IO[str] | None
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)[source]
property P_ani
add_to_density(nt_sR, D_asii)[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], numpy.ndarray]]) – Center of molecule. Defaults to center of cell.

Returns

matrix elements in atomic units.

Return type

Array3D

force_contribution(dH_asii, F_av)[source]
gather_wave_function_coefficients()[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: gpaw.core.atom_centered_functions.AtomCenteredFunctions
subspace_diagonalize(Ht, dH, work_array=None, Htpsit_nX=None, scalapack_parameters=(None, 1, 1, - 1))[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\]
class gpaw.new.ase_interface.ASECalculator(params, log, calculation=None, atoms=None)[source]

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

calculate(atoms)[source]
calculate_property(atoms, prop)[source]

Calculate (if not already calculated) a property.

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.

create_new_calculation(atoms)[source]
get_atomic_electrostatic_potentials()[source]
get_atoms()[source]
get_bz_k_points()[source]
get_dipole_moment(atoms)[source]
get_eigenvalues(kpt=0, spin=0)[source]
get_fermi_level()[source]
get_forces(atoms)[source]
get_homo_lumo(spin=None)[source]
get_ibz_k_points()[source]
get_magnetic_moment(atoms)[source]
get_magnetic_moments(atoms)[source]
get_number_of_bands()[source]
get_number_of_electrons()[source]
get_number_of_iterations()[source]
get_potential_energy(atoms, force_consistent=False)[source]
get_pseudo_wave_function(n)[source]
get_reference_energy()[source]
get_stress(atoms)[source]
move_atoms(atoms)[source]
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, **kwargs)[source]

Create ASE-compatible GPAW calculator.

FFTW

Python wrapper for FFTW3 library

class gpaw.fftw.FFTPlans(size_c, dtype)[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)[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)[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?