Coverage for gpaw/utilities/sic.py: 100%

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

14 Helper object for applying non-self-consistent self-interaction 

15 corrections. 

16 

17 :param paw: 

18 Converged GPAW calculator 

19 :type paw: :class:`gpaw.calculator.GPAW` 

20 :param \**sic_parameters: 

21 Keyword arguments used to initialize the :class:`gpaw.xc.sic.SIC` 

22 object, defaults to :py:attr:`~sic_defaults` 

23 """ 

24 # NOTE: I'm only guessing what this part of the code does. Pls get someone 

25 # more competent to review. 

26 

27 def __init__(self, paw, **sic_parameters): 

28 self.paw = paw 

29 self.sic_parameters = dict(self.sic_defaults, **sic_parameters) 

30 

31 def calculate(self): 

32 ESIC = 0 

33 xc = self.paw.hamiltonian.xc 

34 assert xc.type == 'LDA' 

35 

36 # Calculate the contribution from the core orbitals 

37 for a in self.paw.density.D_asp: 

38 setup = self.paw.density.setups[a] 

39 # TODO: Use XC which has been used to calculate the actual 

40 # calculation. 

41 # TODO: Loop over setups, not atoms. 

42 print('Atom core SIC for ', setup.symbol) 

43 print('%10s%10s%10s' % ('E_xc[n_i]', 'E_Ha[n_i]', 'E_SIC')) 

44 g = Generator(setup.symbol, xcname='LDA', nofiles=True, txt=None) 

45 g.run(**parameters[setup.symbol]) 

46 njcore = g.njcore 

47 for f, l, e, u in zip(g.f_j[:njcore], g.l_j[:njcore], 

48 g.e_j[:njcore], g.u_j[:njcore]): 

49 # Calculate orbital density 

50 # NOTE: It's spherically symmetrized! 

51 # n = np.dot(self.f_j, 

52 assert l == 0, ('Not tested for l>0 core states') 

53 na = np.where(abs(u) < 1e-160, 0, u)**2 / (4 * pi) 

54 na[1:] /= g.r[1:]**2 

55 na[0] = na[1] 

56 nb = np.zeros(g.N) 

57 v_sg = np.zeros((2, g.N)) 

58 vHr = np.zeros(g.N) 

59 Exc = xc.calculate_spherical(g.rgd, np.array([na, nb]), v_sg) 

60 hartree(0, na * g.r * g.dr, g.r, vHr) 

61 EHa = 2 * pi * np.dot(vHr * na * g.r, g.dr) 

62 print('{:10.2f}{:10.2f}{:10.2f}'.format(Exc * Hartree, 

63 EHa * Hartree, 

64 -f * (EHa + 

65 Exc) * Hartree)) 

66 ESIC += -f * (EHa + Exc) 

67 

68 sic = SIC(**self.sic_parameters) 

69 sic.initialize(self.paw.density, self.paw.hamiltonian, self.paw.wfs) 

70 sic.set_positions(self.paw.spos_ac) 

71 

72 print('Valence electron sic ') 

73 print('%10s%10s%10s%10s%10s%10s' % ('spin', 'k-point', 'band', 

74 'E_xc[n_i]', 'E_Ha[n_i]', 'E_SIC')) 

75 assert len(self.paw.wfs.kpt_u) == 1, 'Not tested for bulk calculations' 

76 

77 for s, spin in sic.spin_s.items(): 

78 spin.initialize_orbitals() 

79 spin.update_optimal_states() 

80 spin.update_potentials() 

81 

82 n = 0 

83 for xc, c in zip(spin.exc_m, spin.ecoulomb_m): 

84 print('%10i%10i%10i%10.2f%10.2f%10.2f' % 

85 (s, 0, n, -xc * Hartree, -c * Hartree, 

86 2 * (xc + c) * Hartree)) 

87 n += 1 

88 

89 ESIC += spin.esic 

90 

91 print('Total correction for self-interaction energy:') 

92 print('%10.2f eV' % (ESIC * Hartree)) 

93 print('New total energy:') 

94 total = (ESIC * Hartree + self.paw.get_potential_energy() + 

95 self.paw.get_reference_energy()) 

96 print('%10.2f eV' % total) 

97 return total 

98 

99 sic_defaults = dict(finegrid=True, 

100 coulomb_factor=1, 

101 xc_factor=1)