1from math import pi
2
3import numpy as np
4from ase.units import Hartree
5
6from gpaw.xc.sic import SIC
7from gpaw.atom.generator import Generator
8from gpaw.atom.configurations import parameters
9from gpaw.utilities import hartree
10
11
12class NSCFSIC:
13    def __init__(self, paw):
14        self.paw = paw
15
16    def calculate(self):
17        ESIC = 0
18        xc = self.paw.hamiltonian.xc
19        assert xc.type == 'LDA'
20
21        # Calculate the contribution from the core orbitals
22        for a in self.paw.density.D_asp:
23            setup = self.paw.density.setups[a]
24            # TODO: Use XC which has been used to calculate the actual
25            # calculation.
26            # TODO: Loop over setups, not atoms.
27            print('Atom core SIC for ', setup.symbol)
28            print('%10s%10s%10s' % ('E_xc[n_i]', 'E_Ha[n_i]', 'E_SIC'))
29            g = Generator(setup.symbol, xcname='LDA', nofiles=True, txt=None)
30            g.run(**parameters[setup.symbol])
31            njcore = g.njcore
32            for f, l, e, u in zip(g.f_j[:njcore], g.l_j[:njcore],
33                                  g.e_j[:njcore], g.u_j[:njcore]):
34                # Calculate orbital density
35                # NOTE: It's spherically symmetrized!
36                # n = np.dot(self.f_j,
37                assert l == 0, ('Not tested for l>0 core states')
38                na = np.where(abs(u) < 1e-160, 0, u)**2 / (4 * pi)
39                na[1:] /= g.r[1:]**2
40                na[0] = na[1]
41                nb = np.zeros(g.N)
42                v_sg = np.zeros((2, g.N))
43                vHr = np.zeros(g.N)
44                Exc = xc.calculate_spherical(g.rgd, np.array([na, nb]), v_sg)
45                hartree(0, na * g.r * g.dr, g.r, vHr)
46                EHa = 2 * pi * np.dot(vHr * na * g.r, g.dr)
47                print(('%10.2f%10.2f%10.2f' % (Exc * Hartree, EHa * Hartree,
48                                               -f * (EHa + Exc) * Hartree)))
49                ESIC += -f * (EHa + Exc)
50
51        sic = SIC(finegrid=True, coulomb_factor=1, xc_factor=1)
52        sic.initialize(self.paw.density, self.paw.hamiltonian, self.paw.wfs)
53        sic.set_positions(self.paw.spos_ac)
54
55        print('Valence electron sic ')
56        print('%10s%10s%10s%10s%10s%10s' % ('spin', 'k-point', 'band',
57                                            'E_xc[n_i]', 'E_Ha[n_i]', 'E_SIC'))
58        assert len(self.paw.wfs.kpt_u) == 1, 'Not tested for bulk calculations'
59
60        for s, spin in sic.spin_s.items():
61            spin.initialize_orbitals()
62            spin.update_optimal_states()
63            spin.update_potentials()
64
65            n = 0
66            for xc, c in zip(spin.exc_m, spin.ecoulomb_m):
67                print(('%10i%10i%10i%10.2f%10.2f%10.2f' %
68                       (s, 0, n, -xc * Hartree, -c * Hartree,
69                        2 * (xc + c) * Hartree)))
70                n += 1
71
72            ESIC += spin.esic
73
74        print('Total correction for self-interaction energy:')
75        print('%10.2f eV' % (ESIC * Hartree))
76        print('New total energy:')
77        total = (ESIC * Hartree + self.paw.get_potential_energy() +
78                 self.paw.get_reference_energy())
79        print('%10.2f eV' % total)
80        return total
81