Table of Contents


2020-12-18 Reference energy

The reference energy is not important while the relative energy matters. One could use separate electrons and nuclei as the reference energy. A better solution is to have isolated atoms as reference.



The setups keyword defines which PAW or pseudopotential file to use for a specific element. For example, Mg has a default 10 valence electrons.

  name: Magnesium
  id: b56b6cab5149d5bf22ad04c2f71a3023
  Z: 12.0
  valence: 10
  core: 2
  charge: 0.0
  file: /comm/groupstacks/gpaw/gpaw-setup/gpaw-setup/0.9.20000/setup/Mg.PBE.gz
  compensation charges: gauss, rc=0.33, lmax=2
  cutoffs: 1.86(filt), 0.54(core),
  valence states:
                energy  radius
    2s(2.00)   -79.840   1.090
    3s(2.00)    -4.705   1.090
    2p(6.00)   -46.613   1.085
    3p(0.00)    -1.330   1.085
    *d           0.000   1.037

  Using partial waves for Mg as LCAO basis

If only 2 valence electrons are needed, one could set setups={'Mg':'2'}. This makes the calculation less accurate but faster. The absolute value of potential energy will be changed after changing the number of valence electrons. After setting this, fmax can still be optimized to a low level such as 0.01 eV/Å during a structure optimization, but with more optimization steps.

  name: Magnesium
  id: 8f9ce7a98bd5637e9449601a5e11df98
  Z: 12.0
  valence: 2
  core: 10
  charge: 0.0
  file: /comm/groupstacks/gpaw/gpaw-setup/gpaw-setup/0.9.20000/setup/Mg.2.PBE.gz
  compensation charges: gauss, rc=0.32, lmax=2
  cutoffs: 1.86(filt), 2.86(core),
  valence states:
                energy  radius
    3s(2.00)    -4.705   1.005
    3p(0.00)    -1.330   1.058
    *s          22.506   1.005

  Using partial waves for Mg as LCAO basis


zero PBC

For the isolated molecules, the cell must be orthogonal if one wants to use pbc=[0, 0, 0].


[gpaw-users] 2D PBC: 2D PBC is only implemented in LCAO and FD modes, but ignored in PW mode. Often, the dipole-layer correction isn't that important.

egg-box effect


gpaw.poisson.BadAxesError: Each axis must be periodic or orthogonal to other axes. But we have pbc=[0 0 0] and orthogonal=[0 0 1]


Density mixing

$$\rho^{in} _ {i+1} = \sum^{nmaxold} \alpha_i (\rho^{in} _ {i} + \beta (\rho^{out} _ {i} - \rho^{in} _ {i}))$$

Exchange correlation

hybrid GGA funcitonal

This is an example of how to use a hybrid GGA funcitonal (PBE0) in the plane wave mode of GPAW. A single point calculation at PBE level is performed first and generates wavefunctions for further PBE0 calculation. In this way, the SCF cycle needs less steps than creating initial wavefunctions by default LCAO scheme and doing calculation directly at PBE0 level.

# working in v23.9.1
from import molecule
from ase.parallel import parprint
from gpaw import GPAW

atoms = molecule('CH3CH2OH')

calc = GPAW(mode={'name': 'pw', 'ecut':400}, kpts=(1, 1, 1), xc='PBE', txt='gpaw.txt')
atoms.calc = calc

calc.set(xc={'name': 'PBE0', 'backend': 'pw'})

vdW functional

optB88-vdW is not so robust in GPAW and less stable than other vdW functionals like (vdW-DF, optPBE-vdW, vdW-DF2). There are two methods to improve the SCF convergence of a optB88-vdW calculation.

calc = GPAW(mode='fd', xc='PBE',
            basis = 'dzp', kpts=(1, 1, 1),
            #mixer=Mixer(beta=0.05, nmaxold=5, weight=50),

opt = BFGS(atoms, trajectory='opt.traj', logfile='opt.log')

These strategies are working well for 'fd' mode, but not so well for the plane wave mode. For example, in the case of coronene on a 7x7 2L graphite, optB88-vdW in the following scheme always failed.

optB86b-vdW is not supported in GPAW (at least until v20.1.0). In principle, one can add something in xc/ and xc/ to use optB86b-vdW.

# gpaw/xc/
        if name in ['vdW-DF', 'vdW-DF2', 'optPBE-vdW', 'optB88-vdW', 'optB86b-vdW',
                    'C09-vdW', 'mBEEF-vdW', 'BEEF-vdW']:
            from gpaw.xc.vdw import VDWFunctional
            return VDWFunctional(name, **kwargs)

# gpaw/xc/
    elif name == 'optB88-vdW':
        kernel = LibXC('GGA_X_OPTB88_VDW+LDA_C_PW')
    elif name == 'optB86b-vdW':
        kernel = LibXC('GGA_X_OPTB86B_VDW+LDA_C_PW')

However, GGA_X_OPTB86B_VDW is only available after LibXC 5.0.0 which currently has some compatibility issues with GPAW.

Custom GGA functional via LibXC

from gpaw.xc.libxc import LibXC
from gpaw.xc.gga import GGA
from gpaw import GPAW
class MyXC(GGA):
    def __init__(self, kernel, setup_name, *args, **kwargs):
        self.setup_name = setup_name
        GGA.__init__(self, kernel, *args, **kwargs)
    def get_setup_name(self):
        return self.setup_name
calc = GPAW(xc=MyXC(LibXC('GGA_X_OPTB88_VDW+LDA_C_PW'), 'revPBE'))


Wannier centers on a benzene molecule

H has one valance electrons (1s) and C has four valance electrons (two 2s and two 2p). Benzene (C6H6) then has 30 valence electrons and 15 occupied bands. By default, GPAW will use 4+1.2*n_occupied_bands (22 bands in the benzene case). Wannier scheme will produce alternating single-double bonds in C-C bonds and need 15 centers (15 bonds which are 6 C-H bonds, 3 C-C bonds and 3 C=C bonds).

GPAW wannier module

# Tested and working on v23.9.1
from import molecule
from ase.parallel import parprint
from gpaw import GPAW
from gpaw.wannier import calculate_overlaps

atoms = molecule('C6H6')

calc = GPAW(mode={'name': 'pw', 'ecut': 400}, kpts=(1, 1, 1),
            xc='PBE', maxiter=200,               
atoms.calc = calc


# n2 is the final number of wannier centers
wan = calculate_overlaps(calc, nwannier=15, n2=15).localize_er(verbose=False)
n_wan_center = len(wan.centers)


new_atoms = wan.centers_as_atoms()


ASE wannier module

In ASE, there is also a wannier module. It is working nicely with finite-difference mode of GPAW

from ase import Atoms
from import molecule
from ase.parallel import parprint
from gpaw import GPAW
from ase.dft.wannier import Wannier

atoms = molecule('C6H6')

#calc = GPAW(mode={'name': 'pw'}, kpts = (1, 1, 1),xc='PBE', txt='gpaw.txt',)
calc = GPAW(mode={'name': 'fd'}, kpts = (1, 1, 1),xc='PBE', txt='gpaw.txt',)
atoms.calc = calc

# Make wannier functions of occupied space only
wan = Wannier(nwannier=15, calc=calc)

centers = wan.get_centers()
n_wan_center = len(centers)

new_atoms = Atoms(f'X{n_wan_center}', positions=centers) + atoms


If the plane-wave mode is used, I got the error (which has been reported in 2016 gpaw-gitlab-issue-#32 and remains unsolved)

rank=1 L14:   File "", line 17, in <module>
rank=1 L15:     wan = Wannier(nwannier=15, calc=calc)
rank=1 L16:           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
rank=1 L17:   File "/comm/swstack/core/python/3.11.1/lib/python3.11/site-packages/ase/dft/", line 329, in __init__
rank=1 L18:     self.Z_dknn[d, k] = calc.get_wannier_localization_matrix(
rank=1 L19:                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
rank=1 L20:   File "/comm/groupstacks/gpaw/gpaw/gpaw-gcc-mkl/23.9.1/lib/python3.11/site-packages/gpaw/", line 1950, in get_wannier_localization_matrix
rank=1 L21:     return self.get_wannier_integrals(spin, kpoint,
rank=1 L22:            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
rank=1 L23:   File "/comm/groupstacks/gpaw/gpaw/gpaw-gcc-mkl/23.9.1/lib/python3.11/site-packages/gpaw/", line 1975, in get_wannier_integrals
rank=1 L24:     Z_nn =[u].psit_nG,
GPAW CLEANUP (node 3): <class 'ValueError'> occurred.  Calling MPI_Abort!
rank=1 L25:            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
rank=1 L26:   File "/comm/groupstacks/gpaw/gpaw/gpaw-gcc-mkl/23.9.1/lib/python3.11/site-packages/gpaw/", line 655, in wannier_matrix
rank=1 L27:     a_nG = (e_G * psit_nG[:nbands].conj()).reshape((nbands, -1))
rank=1 L28:             ~~~~^~~~~~~~~~~~~~~~~~~~~~~~~
rank=1 L29: ValueError: operands could not be broadcast together with shapes (20,22,24) (22,826) 

Here e_G has a shape of (20,22,24) and psit_nG[:nbands].conj() has a shape of (22,826).

GPAW-Wannier90 interface

The first step is to get a wavefunction file from GPAW (plane wave is also working)

# Tested and working on v23.9.1 
from import molecule
from gpaw import GPAW

atoms = molecule('C6H6')

calc = GPAW(mode={'name': 'pw'}, xc='PBE', txt='log-step1.txt')
atoms.calc = calc

calc.write('benzene.gpw', mode='all')

The second step is to generate necessary wannier input files and run a wannier calculation

# Tested and working on v23.9.1
import os
import numpy as np
from gpaw.wannier90 import Wannier90
from gpaw import GPAW

seed = 'benzene'

calc = GPAW(seed + '.gpw', txt='step2.txt')

# atomic orbital indices
orbitals_ai = [
    [1, 2],
    [1, 2, 3],
    [1, 2],
    [1, 2, 3],
    [1, 2],
    [1, 2, 3],
#n_atoms = 12
#Nw = np.sum([len(orbitals_ai[ai]) for ai in range(n_atoms)])

w90 = Wannier90(calc,
                bands=range(15), # define num_bands in .win file

# generate .win file
w90.write_input(num_iter=1000, write_xyz=True, translate_home_cell=True)

# generate .nnkp file
os.system('wannier90.x -pp ' + seed)

w90.write_projections() # .amn file, initial guess A^k_{mn} for the unitary transformation, dependent on orbitals_ai
w90.write_eigenvalues() # .eig file, eigenvalues for each band (up to num_bands as defined in .win file)
w90.write_overlaps() # read .nnkp file and generate .mmn overlap file M^{(k,b)}_{mn}

os.system('wannier90.x ' + seed)

One can get atomic orbital indices via

from gpaw.utilities.dos import print_projectors
 i n l m
 0 2 s_1
 1 2 p_y
 2 2 p_z
 3 2 p_x
 4 * s_1
 5 * p_y
 6 * p_z
 7 * p_x
 8 * d_xy
 9 * d_yz
10 * d_3z^2-r^2
11 * d_xz
12 * d_x^2-y^2

Wannier centers of H2O

Similarly, H2O has 8 valence electrons and needs 4 wannier centers.


virtualenv --system-site-packages py3-ase-3.18.2

source /home/tang/python-virtualenv/py3-ase-3.18.2/bin/activate

git clone -b 19.8.1 gpaw-py3.6.3-v19.8.1-intel-2019

module load intel/2019.5.281
module load mpich/3.2
module load fftw/3.3.6-pl2

cd /home/tang/opt/gpaw/gpaw-py3.6.3-v19.8.1-intel-2019
python build_ext 2>&1 | tee gpaw_compile.txt
pip install .

GPAW v19.8.1 uses the following

compiler = 'mpicc'
mpicompiler = 'mpicc'
mpilinker = 'mpicc'

libraries = []
library_dirs = []
include_dirs = []
extra_link_args = []
extra_compile_args = ['-xSSE4.2',

# Use Intel MKL
libraries += ['mkl_sequential', 'mkl_core', 'mkl_rt', ]

# FFTW3:
fftw = True
if fftw:
    libraries += ['fftw3']

# ScaLAPACK (version 2.0.1+ required):
scalapack = True
if scalapack:
    libraries += ['mkl_scalapack_lp64', 'mkl_blacs_intelmpi_lp64']
    mpi_define_macros += [('GPAW_NO_UNDERSCORE_CBLACS', '1')]
    mpi_define_macros += [('GPAW_NO_UNDERSCORE_CSCALAPACK', '1')]

# - static linking:
if 0:
    xc = '/home/tang/opt/libxc/libxc-4.3.4-intel-2019-gpaw-20.1/'
    include_dirs += [xc + 'include']
    extra_link_args += [xc + 'lib/libxc.a']
    if 'xc' in libraries:

# - dynamic linking (requires rpath or setting LD_LIBRARY_PATH at runtime):
if 1:
    xc = '/home/tang/opt/libxc/libxc-4.3.4-intel-2019-gpaw-20.1/'
    include_dirs += [xc + 'include']
    library_dirs += [xc + 'lib']
    # You can use rpath to avoid changing LD_LIBRARY_PATH:
    extra_link_args += ['-Wl,-rpath={xc}/lib'.format(xc=xc)]
    if 'xc' not in libraries: