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.
Mg-setup: 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.
Mg-setup: 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
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.
The finite grid spacing in real-space DFT calculations leads to the so-called
egg-box effectwhere the total energy of a system depends on the positions of all atoms in the system relative to the grid points
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]
atoms.set_pbc([1, 1, 0]) for non orthogonal cellatoms.set_pbc([0, 0, 0]) and atoms.set_cell([a, b, c]) at the same timedav is stable and efficient for most cases.rmm-diis is very sensitive to the initial guess of the wavefunction.$$\rho^{in} _ {i+1} = \sum^{nmaxold} \alpha_i (\rho^{in} _ {i} + \beta (\rho^{out} _ {i} - \rho^{in} _ {i}))$$
XC_FAMILY_LDA, XC_FAMILY_GGA, XC_FAMILY_MGGA, XC_FAMILY_HYB_GGA in gpaw/c/xc/libxc.c which only checks is_gga() and is_mgga().XC_FAMILY_HYB_LDA, XC_FAMILY_HYB_MGGA, XC_FAMILY_LCA in libxc/examples/xcinfo.f90This 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 ase.build import molecule from ase.parallel import parprint from gpaw import GPAW atoms = molecule('CH3CH2OH') atoms.center(5) calc = GPAW(mode={'name': 'pw', 'ecut':400}, kpts=(1, 1, 1), xc='PBE', txt='gpaw.txt') atoms.calc = calc parprint(atoms.get_potential_energy()) parprint(atoms.get_forces()) calc.set(xc={'name': 'PBE0', 'backend': 'pw'}) parprint(atoms.get_potential_energy()) parprint(atoms.get_forces())
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.
PBE single point calculation first, then switch to optB88-vdW for the structure optimizationmix = Mixer(beta=0.05, nmaxold=5, weight=50), use smaller beta and nmaxoldcalc = GPAW(mode='fd', xc='PBE', basis = 'dzp', kpts=(1, 1, 1), #mixer=Mixer(beta=0.05, nmaxold=5, weight=50), txt='gpaw.txt') atoms.set_calculator(calc) atoms.get_potential_energy() atoms.calc.set(xc='optB88-vdW') opt = BFGS(atoms, trajectory='opt.traj', logfile='opt.log') opt.run(fmax=0.01)
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/init.py and xc/vdw.py to use optB86b-vdW.
# gpaw/xc/__init__.py ... 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/vdw.py ... 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.
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'))
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).
# Tested and working on v23.9.1 from ase.build import molecule from ase.parallel import parprint from gpaw import GPAW from gpaw.wannier import calculate_overlaps atoms = molecule('C6H6') atoms.center(vacuum=2.5) calc = GPAW(mode={'name': 'pw', 'ecut': 400}, kpts=(1, 1, 1), xc='PBE', maxiter=200, txt='gpaw.txt', ) atoms.calc = calc atoms.get_potential_energy() # 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) parprint(n_wan_center) new_atoms = wan.centers_as_atoms() new_atoms.write('benzene-wannier.traj') parprint(new_atoms)
In ASE, there is also a wannier module. It is working nicely with finite-difference mode of GPAW
from ase import Atoms from ase.build import molecule from ase.parallel import parprint from gpaw import GPAW from ase.dft.wannier import Wannier atoms = molecule('C6H6') atoms.center(vacuum=2.5) #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 atoms.get_potential_energy() # Make wannier functions of occupied space only wan = Wannier(nwannier=15, calc=calc) wan.localize() centers = wan.get_centers() n_wan_center = len(centers) parprint(n_wan_center) new_atoms = Atoms(f'X{n_wan_center}', positions=centers) + atoms new_atoms.write('benzene-wannier.traj') parprint(new_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 "calc_opt_gpaw_in_one.py", 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/wannier.py", 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/calculator.py", 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/calculator.py", line 1975, in get_wannier_integrals rank=1 L24: Z_nn = self.wfs.gd.wannier_matrix(kpt_u[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/grid_descriptor.py", 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).
The first step is to get a wavefunction file from GPAW (plane wave is also working)
# Tested and working on v23.9.1 from ase.build import molecule from gpaw import GPAW atoms = molecule('C6H6') atoms.center(vacuum=3.5) calc = GPAW(mode={'name': 'pw'}, xc='PBE', txt='log-step1.txt') atoms.calc = calc atoms.get_potential_energy() 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)]) #print(Nw) w90 = Wannier90(calc, seed=seed, bands=range(15), # define num_bands in .win file orbitals_ai=orbitals_ai) # 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 print_projectors('C') ''' 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 '''
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 https://gitlab.com/gpaw/gpaw.git 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 setup.py build_ext 2>&1 | tee gpaw_compile.txt pip install .
GPAW v19.8.1 uses the following customize.py
compiler = 'mpicc' mpicompiler = 'mpicc' mpilinker = 'mpicc' libraries = [] library_dirs = [] include_dirs = [] extra_link_args = [] extra_compile_args = ['-xSSE4.2', '-O3', '-ipo', '-no-prec-div', '-fPIC'] # 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: libraries.remove('xc') # - 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: libraries.append('xc')