Overview
This document describes the most important objects used for a DFT calculation. More information can be found in the code.
PAW
This object is the central object for a GPAW calculation:
+----------+
|GPAWLogger| +-----------+
+----------+ --->|Hamiltonian|
^ / +-----------+
| ---- +------+
| / ---->|Setups|
+-----+ +------+ / +------+
|Atoms|<-------------| GPAW |-----
+-----+ +------+ \
/ | \ \ +-----------+
+-------------+ / | --- ----------->|Occupations|
|WaveFunctions|<-- v \ +-----------+
+-------------+ +-------+ \ +-------+
|Density| -->|SCFLoop|
+-------+ +-------+
The implementation is in gpaw/calculator.py. The
GPAW
class doesn’t do any part of the actual
calculation - it only handles the logic of parsing the input
parameters and setting up the necessary objects for doing the actual
work (see figure above).
A GPAW instance has the following attributes: atoms
,
parameters
, wfs
, density
, setups
,
hamiltonian
, scf
, log
, timer
,
occupations
, initialized
, world
and observers
.
The GPAW
inherits from:
ase.calculators.calculator.Calculator
This implements the ASE calculator interface.
Note
GPAW uses atomic units internally (\(\hbar=e=m=1\)) and ASE uses
Ångström and eV (units
).
Generating a GPAW instance from scratch
When a GPAW instance is created from scratch:
calc = GPAW(mode='fd', xc='LDA', nbands=7)
the GPAW object is almost empty. In order to start a calculation, one will have to do something like:
atoms = Atoms(...)
atoms.calc = calc
atoms.get_potential_energy()
ASE will then arrange to call the calculate()
method with the correct arguments. This will trigger:
A call to the
initialize()
method, which will set up the objects needed for a calculation:Density
,Hamiltonian
,WaveFunctions
,Setups
and a few more (see figure above).A call to the
set_positions()
method, which will initialize everything that depends on the atomic positions:Pass on the atomic positions to the wave functions, Hamiltonian and density objects (call their
set_positions()
methods).Make sure the wave functions are initialized.
Reset the
SCFLoop
.
Generating a GPAW instance from a restart file
When a GPAW instance is created like this:
calc = GPAW('restart.gpw')
the initialize()
method is called first, so that
the parts read from the file can be placed inside the objects where they
belong: the effective pseudo potential and the total energy are put in the
Hamiltonian, the pseudo density is put in the density object and so on.
After a restart, everything should be as before the restart file was written. However, there are a few exceptions:
The wave functions are only read when needed … XXX
Atom centered functions (\(\tilde{p}_i^a\), \(\bar{v}^a\), \(\tilde{n}_c^a\) and \(\hat{g}_{\ell m}^a\)) are not initialized. … XXX
WaveFunctions
We currently have two representations for the wave functions: uniform
3-d grids and expansions in atom centered basis functions as
implemented in the two classes
FDWaveFunctions
and
LCAOWaveFunctions
. Both inherit from the
WaveFunctions
class, so the wave
functions object will always have a
GridDescriptor
, an
Eigensolver
, a
Setups
object and a list of KPoint
objects.
+--------------+ +-----------+
|GridDescriptor| |Eigensolver|
+--------------+ +-----------+
^ ^
|gd |
\ |
+------+ +-------------+ kpt_u +------+
|Setups|<-------|WaveFunctions|-------->|KPoint|+
+------+ +-------------+ +------+|+
^ +------+|
/_\ +------+
|
|
--------------------------------
| |
+-----------------+ +-----------------+
|LCAOWaveFunctions| | FDWaveFunctions |
+-----------------+ +-----------------+
| | / | |
v |tci | |kin |pt
+--------------+ | v | v
|BasisFunctions| | +-------+ | +----------+
+--------------+ | |Overlap| | |Projectors|
v +-------+ | +----------+
+------------------+ v
|TwoCenterIntegrals| +---------------------+
+------------------+ |KineticEnergyOperator|
+---------------------+
Attributes of the wave function object: gd
, nspins
,
nbands
, mynbands
, dtype
, world
,
kpt_comm
, band_comm
, gamma
, bzk_kc
,
ibzk_kc
, weight_k
, symmetry
, kpt_comm
,
rank_a
, nibzkpts
, kpt_u
, setups
,
ibzk_qc
, eigensolver
, and timer
.
Exchange-correlation functionals module
The gpaw.xc
module contains all the code for XC functionals in
GPAW:
+--------------+
| XCFunctional |
+--------------+
^ ^
/_\ /_\
| |
+-------+ | +------------------------+
| LDA | ----|vdW-DF/HybridXC/SIC/GLLB|
+-------+ +------------------------+
^
/_\
|
+---+
|GGA|
+---+
^
/_\
|
+----+
|MGGA|
+----+
An XCFunctional
object is usually created
using the gpaw.xc.XC()
function:
- gpaw.xc.XC(kernel, parameters=None, atoms=None, collinear=True, xp=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/gpaw/envs/latest/lib/python3.12/site-packages/numpy/__init__.py'>)[source]
Create XCFunctional object.
- kernel: XCKernel object, dict or str
Kernel object or name of functional.
- parameters: ndarray
Parameters for BEE functional.
Recognized names are: LDA, PW91, PBE, revPBE, RPBE, BLYP, HCTH407, TPSS, M06-L, revTPSS, vdW-DF, vdW-DF2, EXX, PBE0, B3LYP, BEE, GLLBSC. One can also use equivalent libxc names, for example GGA_X_PBE+GGA_C_PBE is equivalent to PBE, and LDA_X to the LDA exchange. In this way one has access to all the functionals defined in libxc. See xc_funcs.h for the complete list.
Warning - if an MGGA from libxc is used, libxc should be compiled with –disable-fhc. Otherwise the calcualtions won’t converge
Example:
# Implementation of PBE from LibXC:
from gpaw.xc import XC
xc = XC('PBE')
# alternative call:
from gpaw.xc.libxc import LibXC
from gpaw.xc.gga import GGA
xc = GGA(LibXC('PBE'))
# or, explicitly:
xc = GGA(LibXC('GGA_X_PBE+GGA_C_PBE'))
In this example, calling the calculate
method of the xc
object passing in a GridDescriptor
, an
input density array and an output array for the potential, the
GGA
object will calculate the gradient of the
density and pass that and the density on to the libxc kernel.
Refer to Exchange-Correlation functional for other examples.
GPAW also has a few non-libxc kernels that one can use like this:
from gpaw.xc.kernel import XCKernel
xc = XC(XCKernel('PBE'))
Naming convention for arrays
A few examples:
name
shape
spos_c
(3,)
Scaled position vector
nt_sG
(2, 24, 24, 24)
Pseudo-density array \(\tilde{n}_\sigma(\vec{r})\) (
t
means tilde): two spins, 24*24*24 grid points.
cell_cv
(3, 3)
Unit cell vectors.
Commonly used indices:
index
description
a
Atom number
c
Unit cell axis-index (0, 1, 2)
v
xyz-index (0, 1, 2)
k
k-point index
q
k-point index (local, i.e. it starts at 0 on each processor)
s
Spin index (\(\sigma\))
u
Combined spin and k-point index (local)
G
Three indices into the coarse 3D grid
g
Three indices into the fine 3D grid
M
LCAO orbital index (\(\mu\))
n
Principal quantum number or band number
l
Angular momentum quantum number (s, p, d, …)
m
Magnetic quantum number (0, 1, …, 2*l - 1)
L
l
andm
(L = l**2 + m
)
j
Valence orbital number (
n
andl
)
i
Valence orbital number (
n
,l
andm
)
q
j1
andj2
pair
p
i1
andi2
pair
r
CPU-rank
Array names and their definition
name in the code |
definition |
wfs.kpt_u[u].P_ani |
\(\langle\tilde{p}_i^a|\tilde{\psi}_{\sigma\mathbf{k}n} \rangle\) |
density.D_asp |
\(D_{s i_1i_2}^a\) |
|
\(\Delta H_{s i_1i_2}^a\) |
setup.Delta_pL |
\(\Delta_{Li_1i_2}^a\) |
setup.M_pp |
\(\Delta C_{i_1i_2i_3i_4}^a\) eq. (C2) in [1] or eq. (47) in [2] |
wfs.kpt_u[u].psit_nG |
\(\tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r})\) |
setup.pt_j |
\(\tilde{p}_j^a(r)\) |
wfs.pt |
\(\tilde{p}_i^a(\mathbf{r}-\mathbf{R}^a)\) |
The Setup
instances are stored in the
Setups
list, shared by the wfs, density, and
Hamiltonian instances. E.g. paw.wfs.setups, paw.density.setups, or
paw.hamiltonian.setups.
Parallelization over spins, k-points domains and states
When using parallelization over spins, k-points, bands and domains, four different MPI communicators are used:
- mpi.world
Communicator containing all processors.
- domain_comm
One domain_comm communicator contains the whole real space domain for a selection of the spin/k-point pairs and bands.
- kpt_comm
One kpt_comm communicator contains all k-points and spin for a selection of bands over part of the real space domain.
- band_comm
One band_comm communicator contains all bands for a selection of k-points and spins over part of the real space domain.
These communicators constitute MPI groups, of which the latter three
are subsets of the world
communicator. The number of members in
the a communicator group is signified by comm.size
. Within each
group, every element (i.e. processor) is assigned a unique index
comm.rank
into the list of processor ids in the group. For instance,
a domain_comm rank of zero signifies that the processor is first in
the group, hence it functions as a domain master.
For an example on how to use an MPI communicator to perform simple data communication, please refer to parallelization.py.
To investigate the way GPAW distributes calculated quantities across the
various MPI groups, simulating an MPI run can be done using gpaw-mpisim
:
$ gpaw-mpisim -v --dry-run=4 --spins=2 --kpoints=4 --bands=3 --domain-decomposition=2,1,1
Simulating: world.size = 4
parsize_c = (2, 1, 1)
parsize_bands = 1
nspins = 2
nibzkpts = 4
nbands = 3
world: rank=0, ranks=None
kpt_comm : rank=0, ranks=[0 2], mynks=4, kpt_u=[0^,1^,2^,3^]
band_comm : rank=0, ranks=[0], mynbands=3, mybands=[0, 1, 2]
domain_comm : rank=0, ranks=[0 1]
world: rank=1, ranks=None
kpt_comm : rank=0, ranks=[1 3], mynks=4, kpt_u=[0^,1^,2^,3^]
band_comm : rank=0, ranks=[1], mynbands=3, mybands=[0, 1, 2]
domain_comm : rank=1, ranks=[0 1]
world: rank=2, ranks=None
kpt_comm : rank=1, ranks=[0 2], mynks=4, kpt_u=[0v,1v,2v,3v]
band_comm : rank=0, ranks=[2], mynbands=3, mybands=[0, 1, 2]
domain_comm : rank=0, ranks=[2 3]
world: rank=3, ranks=None
kpt_comm : rank=1, ranks=[1 3], mynks=4, kpt_u=[0v,1v,2v,3v]
band_comm : rank=0, ranks=[3], mynbands=3, mybands=[0, 1, 2]
domain_comm : rank=1, ranks=[2 3]
For the case of a \(\Gamma\)-point calculation without band-parallelization, all parallel communication is done in the one domain_comm communicator, which in this case is equal to mpi.world.