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

1import pytest 

2from ase import Atom, Atoms 

3from gpaw import GPAW, Davidson, Mixer 

4from gpaw.xc.hybrid import HybridXC 

5 

6 

7@pytest.mark.mgga 

8@pytest.mark.libxc 

9def test_xc_atomize(in_tmp_dir, gpaw_new): 

10 def xc(name): 

11 return {'name': name, 'stencil': 1} 

12 

13 a = 6. # Size of unit cell (Angstrom) 

14 c = a / 2 

15 # Hydrogen atom: 

16 atom = Atoms([Atom('H', (c, c, c), magmom=1)], 

17 cell=(a, a, a), pbc=False) 

18 

19 # gpaw calculator: 

20 calc = GPAW(mode='fd', 

21 gpts=(32, 32, 32), 

22 nbands=1, 

23 xc=xc('PBE'), 

24 txt='H.txt', 

25 eigensolver=Davidson(12), 

26 mixer=Mixer(0.5, 5), 

27 parallel=dict(kpt=1), 

28 convergence=dict(eigenstates=3.3e-8)) 

29 atom.calc = calc 

30 

31 e1 = atom.get_potential_energy() 

32 de1t = calc.get_xc_difference(xc('TPSS')) 

33 de1m = calc.get_xc_difference(xc('M06-L')) 

34 if not gpaw_new: 

35 de1x = calc.get_xc_difference( 

36 HybridXC('EXX', stencil=1, finegrid=True)) 

37 de1xb = calc.get_xc_difference( 

38 HybridXC('EXX', stencil=1, finegrid=False)) 

39 

40 # Hydrogen molecule: 

41 d = 0.74 # Experimental bond length 

42 molecule = Atoms([Atom('H', (c - d / 2, c, c)), 

43 Atom('H', (c + d / 2, c, c))], 

44 cell=(a, a, a), pbc=False) 

45 

46 molecule.calc = calc.new(txt='H2.txt') 

47 e2 = molecule.get_potential_energy() 

48 de2t = molecule.calc.get_xc_difference(xc('TPSS')) 

49 de2m = molecule.calc.get_xc_difference(xc('M06-L')) 

50 

51 print('hydrogen atom energy: %5.2f eV' % e1) 

52 print('hydrogen molecule energy: %5.2f eV' % e2) 

53 print('atomization energy: %5.2f eV' % (2 * e1 - e2)) 

54 print('atomization energy TPSS: %5.2f eV' % 

55 (2 * (e1 + de1t) - (e2 + de2t))) 

56 print('atomization energy M06-L: %5.2f eV' % 

57 (2 * (e1 + de1m) - (e2 + de2m))) 

58 PBETPSSdifference = (2 * e1 - e2) - (2 * (e1 + de1t) - (e2 + de2t)) 

59 PBEM06Ldifference = (2 * e1 - e2) - (2 * (e1 + de1m) - (e2 + de2m)) 

60 print(PBETPSSdifference) 

61 print(PBEM06Ldifference) 

62 # TPSS value is from JCP 120 (15) 6898, 2004 

63 # e.g. Table VII: DE(PBE - TPSS) = (104.6-112.9)*kcal/mol 

64 # EXX value is from PRL 77, 3865 (1996) 

65 assert PBETPSSdifference == pytest.approx(-0.3599, abs=0.04) 

66 assert PBEM06Ldifference == pytest.approx(-0.169, abs=0.01) 

67 

68 energy_tolerance = 0.002 

69 assert e1 == pytest.approx(-1.081638, abs=energy_tolerance) 

70 assert e2 == pytest.approx(-6.726356, abs=energy_tolerance) 

71 

72 if not gpaw_new: 

73 de2x = molecule.calc.get_xc_difference( 

74 HybridXC('EXX', stencil=1, finegrid=True)) 

75 de2xb = molecule.calc.get_xc_difference( 

76 HybridXC('EXX', stencil=1, finegrid=False)) 

77 

78 print('atomization energy EXX: %5.2f eV' % 

79 (2 * (e1 + de1x) - (e2 + de2x))) 

80 print('atomization energy EXX: %5.2f eV' % 

81 (2 * (e1 + de1xb) - (e2 + de2xb))) 

82 PBEEXXdifference = (2 * e1 - e2) - (2 * (e1 + de1x) - (e2 + de2x)) 

83 PBEEXXbdifference = (2 * e1 - e2) - (2 * (e1 + de1xb) - (e2 + de2xb)) 

84 print(PBEEXXdifference) 

85 print(PBEEXXbdifference) 

86 assert PBEEXXdifference == pytest.approx(0.91, abs=0.005) 

87 assert PBEEXXbdifference == pytest.approx(0.91, abs=0.005)