Coverage for gpaw/test/test_dipole.py: 96%

1""" 

2Test the dipole correction code by comparing this system: 

3 

4 H 

5z1 O z2 

6 H 

7 

8 

9(where z1 and z2 denote points where the potential is probed) 

10 

11Expected potential: 

12 

13 ----- 

14 / 

15 / 

16---- 

17 

18to this system: 

19 

20 H H 

21z1 O z2 O 

22 H H 

23 

24 

25Expected potential: 

26 

27 ------- 

28 / \ 

29 / \ 

30----- ------ 

31 

32The height of the two potentials are tested to be the same. 

33 

34Enable if-statement in the bottom for nice plots 

35""" 

36import numpy as np 

37from ase.build import molecule 

38from ase.units import Hartree 

39from gpaw import GPAW, Mixer 

40from gpaw.mpi import rank 

41 

42 

43def test_dipole(): 

44 system1 = molecule('H2O') 

45 system1.set_pbc((True, True, False)) 

46 system1.cell = 4.0 * np.array([[1.0, -1.5, 0.0], [1.0, 1.0, 0.0], 

47 [0., 0., 1.]]) 

48 system1.center(vacuum=10.0, axis=2) 

49 

50 system2 = system1.copy() 

51 system2.positions *= [1.0, 1.0, -1.0] 

52 system2 += system1 

53 system2.center(vacuum=6.0, axis=2) 

54 

55 convergence = dict(density=1e-6) 

56 

57 def kw(): 

58 return dict(mode='lcao', 

59 convergence=convergence, 

60 mixer=Mixer(0.5, 7, 50.0), 

61 h=0.25, 

62 xc='oldLDA') 

63 

64 calc1 = GPAW(poissonsolver={'dipolelayer': 'xy'}, **kw()) 

65 

66 system1.calc = calc1 

67 system1.get_potential_energy() 

68 v1 = calc1.get_effective_potential() 

69 

70 calc2 = GPAW(**kw()) 

71 

72 system2.calc = calc2 

73 system2.get_potential_energy() 

74 v2 = calc2.get_effective_potential() 

75 

76 def get_avg(v): 

77 nx, ny, nz = v.shape 

78 vyz = v.sum(axis=0) / nx 

79 vz = vyz.sum(axis=0) / ny 

80 return vz, vyz 

81 

82 if rank == 0: 

83 vz1, vyz1 = get_avg(v1) 

84 vz2, vyz2 = get_avg(v2) 

85 

86 # Compare values that are not quite at the end of the array 

87 # (at the end of the array things can "oscillate" a bit) 

88 dvz1 = vz1[-3] - vz1[3] 

89 dvz2 = vz2[3] - vz2[len(vz2) // 2] 

90 print(dvz1, dvz2) 

91 

92 err1 = abs(dvz1 - dvz2) 

93 

94 try: 

95 correction = calc1.hamiltonian.poisson.correction 

96 except AttributeError: 

97 correction = (calc1.dft.pot_calc.poisson_solver. 

98 solver.correction) 

99 

100 correction_err = abs(2.0 * correction * Hartree + dvz1) 

101 print('correction error %s' % correction_err) 

102 assert correction_err < 3e-5 

103 

104 # Comparison to what the values were when this test was last modified: 

105 ref_value = 2.07342988218 

106 err2 = abs(dvz1 - ref_value) 

107 

108 if 0: 

109 import matplotlib.pyplot as plt 

110 plt.imshow(vyz1) 

111 plt.figure() 

112 plt.imshow(vyz2) 

113 plt.figure() 

114 plt.plot(vz1) 

115 plt.plot(vz2) 

116 plt.show() 

117 

118 print('Ref value of previous calculation', ref_value) 

119 print('Value in this calculation', dvz1) 

120 

121 # fine grid needed to achieve convergence! 

122 print('Error', err1, err2) 

123 assert err1 < 4e-3, err1 

124 assert err2 < 2e-4, err2