Coverage for gpaw/xc/tools.py: 93%

1import numpy as np 

2from ase.units import Ha 

3 

4from gpaw.xc import XC 

5from gpaw.utilities import pack_hermitian, unpack_density, unpack_hermitian 

6 

7 

8def vxc(gs, xc=None, coredensity=True, n1=0, n2=0): 

9 """Calculate XC-contribution to eigenvalues.""" 

10 

11 if n2 <= 0: 

12 n2 += gs.bd.nbands 

13 

14 ham = gs.hamiltonian 

15 dens = gs.density 

16 

17 if xc is None: 

18 xc = ham.xc 

19 elif isinstance(xc, str): 

20 xc = XC(xc) 

21 

22 if dens.nt_sg is None: 

23 dens.interpolate_pseudo_density() 

24 

25 thisisatest = not True 

26 

27 if xc.orbital_dependent: 

28 gs.get_xc_difference(xc) 

29 

30 # Calculate XC-potential: 

31 vxct_sg = ham.finegd.zeros(gs.nspins) 

32 xc.calculate(dens.finegd, dens.nt_sg, vxct_sg) 

33 vxct_sG = ham.restrict_and_collect(vxct_sg) 

34 if thisisatest: 

35 vxct_sG[:] = 1 

36 

37 # ... and PAW corrections: 

38 dvxc_asii = {} 

39 for a, D_sp in dens.D_asp.items(): 

40 dvxc_sp = np.zeros_like(D_sp) 

41 

42 pawdata = gs.pawdatasets.by_atom[a] 

43 if pawdata.hubbard_u is not None: 

44 _, dHU_sii = pawdata.hubbard_u.calculate(pawdata, 

45 unpack_density(D_sp)) 

46 dvxc_sp += pack_hermitian(dHU_sii) 

47 xc.calculate_paw_correction(pawdata, D_sp, dvxc_sp, a=a, 

48 addcoredensity=coredensity) 

49 

50 dvxc_asii[a] = [unpack_hermitian(dvxc_p) for dvxc_p in dvxc_sp] 

51 if thisisatest: 

52 dvxc_asii[a] = [pawdata.dO_ii] 

53 

54 vxc_un = np.empty((gs.kd.mynk * gs.nspins, gs.bd.mynbands)) 

55 for u, vxc_n in enumerate(vxc_un): 

56 kpt = gs.kpt_u[u] 

57 vxct_G = vxct_sG[kpt.s] 

58 for n in range(gs.bd.mynbands): 

59 if n1 <= n + gs.bd.beg < n2: 

60 psit_G = gs.get_wave_function_array(u, n) 

61 vxc_n[n] = gs.gd.integrate((psit_G * psit_G.conj()).real, 

62 vxct_G, global_integral=False) 

63 

64 for a, dvxc_sii in dvxc_asii.items(): 

65 m1 = max(n1, gs.bd.beg) - gs.bd.beg 

66 m2 = min(n2, gs.bd.end) - gs.bd.beg 

67 if m1 < m2: 

68 P_ni = kpt.P_ani[a][m1:m2] 

69 vxc_n[m1:m2] += ((P_ni @ dvxc_sii[kpt.s]) * 

70 P_ni.conj()).sum(1).real 

71 

72 gs.gd.comm.sum(vxc_un) 

73 vxc_skn = gs.kd.collect(vxc_un, broadcast=True) 

74 

75 if xc.orbital_dependent: 

76 vxc_skn += xc.exx_skn 

77 

78 vxc_skn = gs.bd.collect(vxc_skn.T.copy(), broadcast=True).T 

79 return vxc_skn[:, :, n1:n2] * Ha