Coverage for gpaw/test/pseudopotential/test_hgh_h2o.py: 74%

1"""Test of HGH pseudopotential implementation. 

2This is the canonical makes-sure-nothing-breaks test, which checks that the 

3numbers do not change from whatever they were before. 

4 

5The test runs an HGH calculation on a misconfigured H2O molecule, such that 

6the forces are nonzero. 

7 

8Energy is compared to a previous calculation; if it differs significantly, 

9that is considered an error. 

10 

11Forces are compared to a previous finite-difference result. 

12""" 

13 

14 

15import numpy as np 

16from ase.build import molecule 

17 

18from gpaw import GPAW, PoissonSolver 

19import pytest 

20 

21 

22def test_pseudopotential_hgh_h2o(): 

23 mol = molecule('H2O') 

24 mol.rattle(0.2) 

25 mol.center(vacuum=2.0) 

26 calc = GPAW(mode='fd', 

27 nbands=6, 

28 poissonsolver=PoissonSolver('fd'), 

29 gpts=(32, 40, 40), 

30 setups='hgh', 

31 convergence=dict(eigenstates=1e-9, 

32 density=1e-5, 

33 energy=0.3e-5), 

34 txt='-') 

35 mol.calc = calc 

36 e = mol.get_potential_energy() 

37 F_ac = mol.get_forces() 

38 

39 F_ac_ref = np.array([[7.33077718, 3.81069249, -6.07405156], 

40 [-0.9079617, -1.18203514, 3.43777589], 

41 [-0.61642527, -0.41889306, 2.332415]]) 

42 

43 eref = -468.527121304 

44 

45 eerr = abs(e - eref) 

46 

47 print('energy', e) 

48 print('ref energy', eref) 

49 print('energy error', eerr) 

50 

51 print('forces') 

52 print(F_ac) 

53 

54 print('ref forces') 

55 print(F_ac_ref) 

56 

57 ferr = np.abs(F_ac - F_ac_ref).max() 

58 print('max force error', ferr) 

59 

60 fdcheck = False 

61 

62 if fdcheck: 

63 from gpaw.test import calculate_numerical_forces 

64 F_ac_fd = calculate_numerical_forces(mol) 

65 

66 print('Self-consistent forces') 

67 print(F_ac) 

68 print('FD forces') 

69 print(F_ac_fd) 

70 print() 

71 print(repr(F_ac_fd)) 

72 print() 

73 err = np.abs(F_ac - F_ac_fd).max() 

74 print('max err', err) 

75 

76 wfs = calc.wfs 

77 gd = wfs.gd 

78 psit_nG = wfs.kpt_u[0].psit_nG 

79 try: 

80 dH_asp = calc.hamiltonian.dH_asp 

81 except AttributeError: 

82 from gpaw.utilities import pack_hermitian 

83 dH_asii = calc.dft.potential.dH_asii 

84 dH_asp = {a: pack_hermitian(dH_sii) for a, dH_sii in dH_asii.items()} 

85 

86 assert eerr < 1e-3, 'energy changed from reference' 

87 assert ferr < 0.015, 'forces do not match FD check' 

88 

89 # Sanity check. In HGH, the atomic Hamiltonian is constant. 

90 # Also the projectors should be normalized 

91 for a, dH_sp in dH_asp.items(): 

92 dH_p = dH_sp[0] 

93 K_p = wfs.setups[a].K_p 

94 # B_ii = wfs.setups[a].B_ii 

95 # assert np.abs(B_ii.diagonal() - 1).max() < 1e-3 

96 # print 'B_ii' 

97 # print wfs.setups[a].B_ii 

98 # Actually, H2O might not be such a good test, since there'll only 

99 # be one element in the atomic Hamiltonian for O and zero for H. 

100 # print 'dH_p', dH_p 

101 # print 'K_p', K_p 

102 

103 assert np.abs(dH_p - K_p).max() < 1e-10, 'atomic Hamiltonian changed' 

104 

105 # h_ii = setup.data.h_ii 

106 # print 'h_ii', h_ii 

107 # print 'dH_ii', dH_ii 

108 

109 # Sanity check: HGH is normconserving 

110 for psit_G in psit_nG: 

111 norm = gd.integrate(psit_G**2) # around 1e-15 ! Surprisingly good. 

112 assert abs(1 - norm) < 1e-10, 'Not normconserving' 

113 

114 energy_tolerance = 0.0003 

115 assert e == pytest.approx(eref, abs=energy_tolerance)