Coverage for gpaw/test/xc/test_lb94.py: 100%

1import pytest 

2from ase import Atoms 

3from gpaw.mpi import world 

4from gpaw import GPAW, FermiDirac, Mixer, MixerSum, Davidson 

5from gpaw.atom.all_electron import AllElectron 

6from gpaw.atom.generator import Generator 

7from gpaw.atom.configurations import parameters 

8 

9 

10@pytest.mark.slow 

11def test_xc_lb94(in_tmp_dir, add_cwd_to_setup_paths): 

12 ref1 = 'R. v. Leeuwen PhysRevA 49, 2421 (1994)' 

13 ref2 = 'Gritsenko IntJQuanChem 76, 407 (2000)' 

14 # HOMO energy in mHa for closed shell atoms 

15 e_HOMO_cs = {'He': 851, 'Be': 321, 'Ne': 788, 

16 'Ar': 577, 'Kr': 529, 'Xe': 474, 

17 'Mg': 281 + 8} 

18 

19 txt = None 

20 

21 print('--- Comparing LB94 with', ref1) 

22 print('and', ref2) 

23 

24 print('**** all electron calculations') 

25 print('atom [refs] -e_homo diff all in mHa') 

26 if world.rank == 0: 

27 for atom in e_HOMO_cs.keys(): 

28 ae = AllElectron(atom, 'LB94', txt=txt) 

29 ae.run() 

30 e_homo = int(ae.e_j[-1] * 10000 + .5) / 10. 

31 diff = e_HOMO_cs[atom] + e_homo 

32 print('%2s %8g %6.1f %4.1g' % 

33 (atom, e_HOMO_cs[atom], -e_homo, diff)) 

34 assert abs(diff) < 6 

35 world.barrier() 

36 

37 setups = {} 

38 

39 print('**** 3D calculations') 

40 print('atom [refs] -e_homo diff all in mHa') 

41 

42 for atom in sorted(e_HOMO_cs): 

43 e_ref = e_HOMO_cs[atom] 

44 # generate setup for the atom 

45 if world.rank == 0 and atom not in setups: 

46 g = Generator(atom, 'LB94', nofiles=True, txt=txt) 

47 g.run(**parameters[atom]) 

48 setups[atom] = 1 

49 world.barrier() 

50 

51 SS = Atoms(atom, cell=(7, 7, 7), pbc=False) 

52 SS.center() 

53 c = GPAW(mode='fd', h=.3, xc='LB94', 

54 eigensolver=Davidson(3), 

55 mixer=Mixer(0.5, 7, 50.0), nbands=-2, txt=txt) 

56 c.calculate(SS) 

57 # find HOMO energy 

58 eps_n = c.get_eigenvalues(kpt=0, spin=0) / 27.211 

59 f_n = c.get_occupation_numbers(kpt=0, spin=0) 

60 for e, f in zip(eps_n, f_n): 

61 if f < 0.99: 

62 break 

63 e_homo = e 

64 e_homo = int(e_homo * 10000 + .5) / 10. 

65 diff = e_ref + e_homo 

66 print('%2s %8g %6.1f %4.1f' % (atom, e_ref, -e_homo, diff)) 

67 assert abs(diff) < 7 

68 

69 # HOMO energy in mHa and magn. mom. for open shell atoms 

70 e_HOMO_os = {'He': [851, 0], # added for cross check 

71 'H': [440, 1], 

72 'N': [534 - 23, 3], 

73 'Na': [189 + 17, 1], 

74 'P': [385 - 16, 3]} 

75 

76 for atom in sorted(e_HOMO_os): 

77 e_ref = e_HOMO_os[atom][0] 

78 magmom = e_HOMO_os[atom][1] 

79 # generate setup for the atom 

80 if world.rank == 0 and atom not in setups: 

81 g = Generator(atom, 'LB94', nofiles=True, txt=txt) 

82 g.run(**parameters[atom]) 

83 setups[atom] = 1 

84 world.barrier() 

85 

86 SS = Atoms(atom, magmoms=[magmom], cell=(7, 7, 7), pbc=False) 

87 SS.center() 

88 # fine grid needed for convergence! 

89 c = GPAW(mode='fd', 

90 h=0.2, 

91 xc='LB94', 

92 nbands=-2, 

93 spinpol=True, 

94 hund=True, 

95 eigensolver=Davidson(3), 

96 mixer=MixerSum(0.5, 7, 50.0), 

97 occupations=FermiDirac(0.0, fixmagmom=True), 

98 txt=txt) 

99 c.calculate(SS) 

100 # find HOMO energy 

101 eps_n = c.get_eigenvalues(kpt=0, spin=0) / 27.211 

102 f_n = c.get_occupation_numbers(kpt=0, spin=0) 

103 for e, f in zip(eps_n, f_n): 

104 if f < 0.99: 

105 break 

106 e_homo = e 

107 e_homo = int(e_homo * 10000 + .5) / 10. 

108 diff = e_ref + e_homo 

109 print('%2s %8g %6.1f %4.1f' % (atom, e_ref, -e_homo, diff)) 

110 assert abs(diff) < 15