[gpaw-users] TPSS band structure

[Markus Meinert]

Dear GPAW users and developers,

I am trying to compute the band structure of bcc Fe with TPSS. I can 
easily compute it with PBE, and I can obtain the ground state and DOS 
with TPSS. However, there seems to be a problem with the band structure, 
which I try to compute according to the Si tutorial. Is this the right 
way to do it? I am using gpaw 0.9.0.8711. Here's the full script, that 
computes the PBE ground state, restarts with TPSS, computes the PBE band 
structure and fails for the TPSS band structure:

################
################

import pickle
import numpy as np
from ase import *
from ase.dft.kpoints import ibz_points, get_bandpath
from ase.parallel import parprint
from gpaw import GPAW, FermiDirac, restart

def groundstate(a, k, xctype):
     atoms = Atoms(
                     'Fe',
                     [(0,0,0)],
                     cell=[(-0.5*a, 0.5*a, 0.5*a),(0.5*a, -0.5*a, 
0.5*a),(0.5*a, 0.5*a, -0.5*a)],
                     pbc=True,
                     magmoms=[2.2]
                   )

     atoms.calc = GPAW(
                         kpts=(k, k, k),
                         xc=xctype,
                         spinpol=True,
                         gpts=(16, 16, 16),
                         stencils=(5,3),
                         usesymm=True,
                         nbands=16,
                         occupations=FermiDirac(0.05),
                         parallel={'domain': 1},
                         convergence={'eigenstates': 1e-7},
                         txt='Fe-PBE.txt'
                       )
     atoms.get_potential_energy()
     return atoms

a = 2.861
k = 20

# run the PBE calculation
parprint('PBE ground state run')
system = groundstate(a, k, 'PBE')
system.calc.write('PBE.gpw', mode='all')

# run the TPSS calculation
parprint('TPSS ground state run')
system, calc = restart('PBE.gpw', txt='Fe-TPSS.txt', 
convergence={'eigenstates': 1e-7})
calc.set(xc='TPSS')
calc.get_potential_energy()
system.calc.write('TPSS.gpw', mode='all')

# calculate the band structure for PBE
parprint('band structure run for PBE')
points = ibz_points['bcc']
G = points['Gamma']
P = points['P']
H = points['H']

calc = GPAW('PBE.gpw',
             txt='PBE-bandstructure.txt',
             spinpol=True,
             parallel={'domain': 1},
             fixdensity=True,
             usesymm=None,
             convergence={'bands': 12},
             eigensolver='cg',
             occupations=FermiDirac(0.05)
            )

kpts, x, X = get_bandpath([G, H, P, G], calc.atoms.cell, npoints=100)
calc.set(kpts=kpts)
calc.get_potential_energy()

ef = calc.get_fermi_level()

e_kn_0 = np.array([calc.get_eigenvalues(k, spin=0)-ef for k in 
range(len(kpts))])
e_kn_1 = np.array([calc.get_eigenvalues(k, spin=1)-ef for k in 
range(len(kpts))])
pickle.dump((x, X, e_kn_0, e_kn_1), open('PBE.pckl', 'w'))

# calculate the band structure for TPSS
parprint('band structure run for TPSS')
calc = GPAW('TPSS.gpw',
             txt='TPSS-bandstructure.txt',
             spinpol=True,
             parallel={'domain': 1},
             fixdensity=True,
             usesymm=None,
             convergence={'bands': 12},
             eigensolver='cg',
             occupations=FermiDirac(0.05)
            )

kpts, x, X = get_bandpath([G, H, P, G], calc.atoms.cell, npoints=100)
calc.set(kpts=kpts)
calc.get_potential_energy()

ef = calc.get_fermi_level()

e_kn_0 = np.array([calc.get_eigenvalues(k, spin=0)-ef for k in 
range(len(kpts))])
e_kn_1 = np.array([calc.get_eigenvalues(k, spin=1)-ef for k in 
range(len(kpts))])
pickle.dump((x, X, e_kn_0, e_kn_1), open('TPSS.pckl', 'w'))

#######################
#######################

Traceback (most recent call last):
   File "b.py", line 28, in <module>
     calc.get_potential_energy()
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/aseinterface.py", line 
37, in get_potential_energy
     self.calculate(atoms, converge=True)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/paw.py", line 245, in 
calculate
     self.set_positions(atoms)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/paw.py", line 303, in 
set_positions
     self.wfs.initialize(self.density, self.hamiltonian, spos_ac)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/wavefunctions/fdpw.py", 
line 70, in initialize
     hamiltonian.update(density)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/hamiltonian.py", line 
202, in update
     self.update_pseudo_potential(density)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/hamiltonian.py", line 
467, in update_pseudo_potential
     Exc = self.xc.calculate(self.finegd, density.nt_sg, self.vt_sg)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/xc/lda.py", line 22, in 
calculate
     self.calculate_lda(e_g, n_sg, v_sg)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/xc/gga.py", line 28, in 
calculate_lda
     self.calculate_gga(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/xc/mgga.py", line 51, 
in calculate_gga
     taut_sG = self.wfs.calculate_kinetic_energy_density(self.taugrad_v)
   File "/home/meinert/gpaw/gpaw-0.9.0.8711/gpaw/wavefunctions/fd.py", 
line 95, in calculate_kinetic_energy_density
     for f, psit_G in zip(kpt.f_n, kpt.psit_nG):
TypeError: zip argument #1 must support iteration
GPAW CLEANUP for serial binary: <type 'exceptions.TypeError'> occured. 
Calling sys.exit()

############################
############################


Thanks a lot in advance,
Markus


-- 
Dipl.-Phys. Markus Meinert
---
Room D2-118
e-mail: meinert at physik.uni-bielefeld.de
Phone: +49 521 106 2661
---
Thin Films and Physics of Nanostructures
Department of Physics - D2
Bielefeld University
Universitätsstraße 25
33615 Bielefeld