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

1import pytest 

2from ase.build import molecule 

3from gpaw import GPAW, Davidson, Mixer, PoissonSolver 

4from gpaw.utilities.tools import split_formula 

5 

6 

7@pytest.mark.mgga 

8def test_tpss(in_tmp_dir): 

9 cell = [6., 6., 7.] 

10 

11 # Reference from J. Chem. Phys. Vol 120 No. 15, 15 April 2004, page 6898 

12 tpss_de = { 

13 'Li2': 22.5} 

14 tpss_old = { 

15 'Li2': 20.9} 

16 

17 exp_bonds_dE = { 

18 'Li2': (2.673, 24.4)} 

19 

20 systems = ['Li2'] 

21 

22 # Add atoms 

23 for formula in systems: 

24 temp = split_formula(formula) 

25 for atom in temp: 

26 if atom not in systems: 

27 systems.append(atom) 

28 energies = {} 

29 

30 # Calculate energies 

31 for formula in systems: 

32 loa = molecule(formula) 

33 loa.set_cell(cell) 

34 loa.center() 

35 calc = GPAW(mode='fd', 

36 h=0.3, 

37 eigensolver=Davidson(8), 

38 parallel=dict(kpt=1), 

39 mixer=Mixer(0.5, 5), 

40 nbands=-2, 

41 poissonsolver=PoissonSolver('fd', relax='GS'), 

42 xc='oldPBE', 

43 txt=formula + '.txt') 

44 if len(loa) == 1: 

45 calc = calc.new(hund=True) 

46 else: 

47 pos = loa.get_positions() 

48 pos[1, :] = pos[0, :] + [0.0, 0.0, exp_bonds_dE[formula][0]] 

49 loa.set_positions(pos) 

50 loa.center() 

51 loa.calc = calc 

52 energy = loa.get_potential_energy() 

53 diff = calc.get_xc_difference('TPSS') 

54 energies[formula] = (energy, energy + diff) 

55 print(formula, energy, energy + diff) 

56 

57 # calculate atomization energies 

58 print('formula\tGPAW\tRef\tGPAW-Ref\tGPAW-exp') 

59 mae_ref, mae_exp, mae_pbe, count = 0.0, 0.0, 0.0, 0 

60 for formula in tpss_de.keys(): 

61 atoms_formula = split_formula(formula) 

62 de_tpss = -1.0 * energies[formula][1] 

63 de_pbe = -1.0 * energies[formula][0] 

64 for atom_formula in atoms_formula: 

65 de_tpss += energies[atom_formula][1] 

66 de_pbe += energies[atom_formula][0] 

67 

68 de_tpss *= 627.5 / 27.211 

69 de_pbe *= 627.5 / 27.211 

70 mae_ref += abs(de_tpss - tpss_de[formula]) 

71 mae_exp += abs(de_tpss - exp_bonds_dE[formula][1]) 

72 mae_pbe += abs(de_pbe - exp_bonds_dE[formula][1]) 

73 count += 1 

74 out = ("%s\t%.1f\t%.1f\t%.1f\t%.1f kcal/mol" % 

75 (formula, de_tpss, tpss_de[formula], 

76 de_tpss - tpss_de[formula], 

77 de_tpss - exp_bonds_dE[formula][1])) 

78 print(out) 

79 

80 # comparison to gpaw revision 5450 version value 

81 # in kcal/mol (note the grid:0.3 Ang) 

82 assert de_tpss == pytest.approx(tpss_old[formula], abs=0.15)