Coverage for gpaw/test/ext_potential/test_external_pw.py: 95%

1import pytest 

2from ase import Atoms 

3from ase.units import Bohr, Ha 

4 

5from gpaw import GPAW, PW 

6from gpaw.external import ConstantElectricField, ConstantPotential 

7 

8 

9@pytest.mark.old_gpaw_only 

10def test_stark_pw(): 

11 h = Atoms('H', pbc=(1, 1, 0), magmoms=[1]) 

12 h.center(vacuum=3.0) 

13 field = 0.2 

14 params = dict(mode=PW(300)) 

15 

16 h.calc = GPAW(**params) 

17 

18 e0 = h.get_potential_energy() 

19 eig0 = h.calc.get_eigenvalues()[0] 

20 

21 h.calc = GPAW(**params, 

22 external=ConstantElectricField(field), 

23 poissonsolver={'dipolelayer': 'xy'}) 

24 

25 e = h.get_potential_energy() 

26 eig = h.calc.get_eigenvalues()[0] 

27 dip = h.get_dipole_moment()[2] 

28 

29 to_au = Ha / Bohr**2 

30 

31 a1 = -2 * (e - e0) / field**2 * to_au 

32 a2 = -(eig - eig0) / field**2 * to_au 

33 a3 = dip / field * to_au 

34 

35 aref = 6.02 

36 

37 for a in [a1, a2, a3]: 

38 print(a) 

39 assert a == pytest.approx(aref, abs=0.1) 

40 

41 

42@pytest.mark.old_gpaw_only 

43def test_ext_potential_external_pw(): 

44 ConstantPotential() 

45 

46 sc = 2.9 

47 R = 0.7 # approx. experimental bond length 

48 R = 1.0 

49 a = 2 * sc 

50 c = 3 * sc 

51 H2 = Atoms('HH', [(a / 2, a / 2, (c - R) / 2), 

52 (a / 2, a / 2, (c + R) / 2)], 

53 cell=(a, a, c), 

54 pbc=False) 

55 

56 txt = None 

57 

58 convergence = {'eigenstates': 1.e-6 * 40 * 1.5**3, 

59 'density': 1.e-2, 

60 'energy': 0.1} 

61 

62 # without potential 

63 if True: 

64 if txt: 

65 print('\n################## no potential') 

66 c00 = GPAW(mode=PW(200), nbands=-1, 

67 convergence=convergence, 

68 txt=txt) 

69 c00.calculate(H2) 

70 c00.get_eigenvalues() 

71 

72 # 0 potential 

73 if True: 

74 if txt: 

75 print('\n################## 0 potential') 

76 cp0 = ConstantPotential(0.0) 

77 c01 = GPAW(mode=PW(200), nbands=-2, external=cp0, 

78 convergence=convergence, 

79 txt=txt) 

80 c01.calculate(H2) 

81 

82 # 1 potential 

83 if True: 

84 if txt: 

85 print('################## 1 potential') 

86 cp1 = ConstantPotential(-1.0) 

87 c1 = GPAW(mode=PW(200), nbands=-2, external=cp1, 

88 convergence=convergence, 

89 txt=txt) 

90 c1.calculate(H2) 

91 

92 for i in range(c00.get_number_of_bands()): 

93 f00 = c00.get_occupation_numbers()[i] 

94 if f00 > 0.01: 

95 e00 = c00.get_eigenvalues()[i] 

96 e1 = c1.get_eigenvalues()[i] 

97 print('Eigenvalues no pot, expected, error=', 

98 e00, e1 + 1, e00 - e1 - 1) 

99 assert e00 == pytest.approx(e1 + 1., abs=0.02) 

100 

101 E_c00 = c00.get_potential_energy() 

102 E_c1 = c1.get_potential_energy() 

103 DeltaE = E_c00 - E_c1 

104 assert DeltaE == pytest.approx(0, abs=0.002)