Coverage for gpaw/test/generic/test_hydrogen.py: 100%

45 statements  

« prev     ^ index     » next       coverage.py v7.7.1, created at 2025-07-14 00:18 +0000

1from math import log 

2 

3import numpy as np 

4import pytest 

5from ase import Atoms 

6from ase.db import connect 

7from ase.io import iread, read, write 

8from ase.units import Bohr 

9 

10from gpaw import GPAW, FermiDirac 

11 

12 

13@pytest.mark.old_gpaw_only 

14def test_generic_hydrogen(in_tmp_dir): 

15 a = 4.0 

16 h = 0.2 

17 hydrogen = Atoms('H', 

18 [(a / 2, a / 2, a / 2)], 

19 cell=(a, a, a)) 

20 

21 params = dict(mode='fd', 

22 h=h, 

23 nbands=1, 

24 convergence={'energy': 1e-7}) 

25 hydrogen.calc = GPAW(**params, 

26 txt='h.txt') 

27 e1 = hydrogen.get_potential_energy() 

28 assert e1 == pytest.approx(0.526939, abs=0.001) 

29 

30 dens = hydrogen.calc.density 

31 c = dens.gd.find_center(dens.nt_sG[0]) * Bohr 

32 assert abs(c - a / 2).max() == pytest.approx(0, abs=1e-13) 

33 

34 kT = 0.001 

35 hydrogen.calc = GPAW(**params, 

36 occupations=FermiDirac(width=kT)) 

37 e2 = hydrogen.get_potential_energy() 

38 assert e1 == pytest.approx(e2 + log(2) * kT, abs=3.0e-7) 

39 

40 # Test ase.db a bit: 

41 for name in ['h.json', 'h.db']: 

42 con = connect(name, append=False) 

43 con.write(hydrogen) 

44 id = con.write(hydrogen, foo='bar', data={'abc': [1, 2, 3]}) 

45 assert id == 2 

46 assert con.reserve(foo='bar') is None 

47 row = con.get(foo='bar') 

48 assert row.energy == e2 

49 assert sum(row.data.abc) == 6 

50 del con[1] 

51 assert con.reserve(x=42) == 3 

52 

53 write('x' + name, hydrogen) 

54 write('xx' + name, [hydrogen, hydrogen]) 

55 

56 assert read(name + '@foo=bar')[0].get_potential_energy() == e2 

57 for n, h in zip([1, 0], iread(name + '@:')): 

58 assert n == len(h) 

59 

60 # Test parsing of GPAW's text output: 

61 h = read('h.txt') 

62 error = abs(h.calc.get_eigenvalues() - 

63 hydrogen.calc.get_eigenvalues()).max() 

64 assert error < 1e-5, error 

65 

66 # Test get_electrostatic_potential() method 

67 v = hydrogen.calc.get_electrostatic_potential() 

68 print(v.shape, np.ptp(v))