[gpaw-users] TPSS band structure

[Lara Ferrighi]

Hei Markus,
  I've never done the band structure myself, but before running a MGGA  
calculation
you need to run something else to create the wave function (for the  
kinetic energy
density needed in MGGAs). My guess would be that you need to run a PBE  
before TPSS
when you scan for Kpoints for the band structure.

Best, Lara

> 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
>