Coverage for gpaw/test/utilities/test_ewald.py: 100%

1import numpy as np 

2from gpaw.utilities.ewald import Ewald 

3 

4 

5def test_utilities_ewald(): 

6 verbose = True 

7 

8 # References for Madelung constants with homogenous background charges: 

9 # M.P. Marder, "Condensed Matter Physics", J. Wiley, (2000), 

10 # scaled to the Wigner-Seitz radius. 

11 # References for Madelung constants of ionic crystals: 

12 # Y. Sakamoto, The J. of Chem. Phys., vol. 28, (1958), p. 164, 

13 # scaled to the nearest neighbor distance. 

14 

15 if 1: # fcc 

16 cell = np.array([[0., .5, .5], [.5, .0, .5], [.5, .5, .0]]) 

17 basis = np.array([[0., 0., 0.]]) 

18 charges = np.array([1.]) 

19 r = np.array([0.0, 0.0, 0.0]) 

20 charges = np.array([1]) 

21 ewald = Ewald(cell) 

22 e = ewald.get_electrostatic_potential(r, basis, charges, 

23 excludefroml0=0) 

24 a_this = -e * ewald.get_wigner_seitz_radius(sum(charges)) 

25 a_ref = 1.79175 

26 if verbose: 

27 print('Madelung energy, fcc:', a_this, a_this - a_ref) 

28 assert abs(a_this - a_ref) < 1e-5 

29 

30 if 1: # hcp 

31 a = 1. 

32 c = np.sqrt(8. / 3.) * a 

33 cell = np.array([[.5 * a, -0.5 * np.sqrt(3) * a, 0], 

34 [.5 * a, 0.5 * np.sqrt(3) * a, 0], [0., 0., c]]) 

35 basis = np.array([[.5 * a, .5 / np.sqrt(3) * a, 0.], 

36 [.5 * a, -.5 / np.sqrt(3) * a, 0.5 * c]]) 

37 r = np.array([.5 * a, .5 / np.sqrt(3) * a, 0.]) 

38 charges = np.array([1., 1.]) 

39 ewald = Ewald(cell) 

40 e = ewald.get_electrostatic_potential(r, basis, charges, 

41 excludefroml0=0) 

42 a_this = -e * ewald.get_wigner_seitz_radius(sum(charges)) 

43 a_ref = 1.79168 

44 if verbose: 

45 print('Madelung energy, hcp:', a_this, a_this - a_ref) 

46 assert abs(a_this - a_ref) < 1e-5 

47 

48 if 1: # simple cubic 

49 cell = np.diag([1., 1., 1.]) 

50 basis = np.array([[0., 0., 0.]]) 

51 charges = np.array([1.]) 

52 r = np.array([0.0, 0.0, 0.0]) 

53 charges = np.array([1]) 

54 ewald = Ewald(cell) 

55 e = ewald.get_electrostatic_potential(r, basis, charges, 

56 excludefroml0=0) 

57 a_this = -e * ewald.get_wigner_seitz_radius(1) 

58 a_ref = 1.76012 

59 if verbose: 

60 print('Madelung energy, sc:', a_this, a_this - a_ref) 

61 assert abs(a_this - a_ref) < 1e-5 

62 

63 if 1: # NaCl 

64 cell = np.array([[0., .5, .5], [.5, .0, .5], [.5, .5, .0]]) 

65 basis = np.array([[0., 0., 0.], [.5, .5, .5]]) 

66 r = np.array([0.0, 0.0, 0.0]) 

67 charges = np.array([1, -1]) 

68 ewald = Ewald(cell) 

69 e = ewald.get_electrostatic_potential(r, basis, charges, 

70 excludefroml0=0) 

71 a_this = -0.5 * e 

72 a_ref = 1.7475645946331822 

73 if verbose: 

74 print('Madelung constant, NaCl (B1):', a_this, a_this - a_ref) 

75 assert abs(a_this - a_ref) < 1e-13 

76 

77 if 0: # NaCl 

78 cell = np.diag((1., 1., 1.)) 

79 basis = np.array([[0., 0., 0.], [.5, .5, .0], 

80 [.5, .0, .5], [.0, .5, .5], 

81 [.5, 0., 0.], [.0, .5, .0], 

82 [.0, .0, .5], [.5, .5, .5]]) 

83 r = np.array([0.0, 0.0, 0.0]) 

84 charges = np.array([1, 1, 1, 1, -1, -1, -1, -1]) 

85 ewald = Ewald(cell) 

86 e = ewald.get_electrostatic_potential(r, basis, charges, 

87 excludefroml0=0) 

88 a_this = -0.5 * e 

89 a_ref = 1.7475645946331822 

90 if verbose: 

91 print('Madelung constant, NaCl (B1):', a_this, a_this - a_ref) 

92 assert abs(a_this - a_ref) < 1e-13 

93 

94 if 1: # CsCl 

95 cell = np.diag((1., 1., 1.)) 

96 basis = np.array([[0., 0., 0.], [.5, .5, .5]]) 

97 r = np.array([0.0, 0.0, 0.0]) 

98 charges = np.array([1, -1]) 

99 ewald = Ewald(cell) 

100 e = ewald.get_electrostatic_potential(r, basis, charges, 

101 excludefroml0=0) 

102 a_this = -.5 * np.sqrt(3.) * e 

103 a_ref = 1.76267477307099 

104 if verbose: 

105 print('Madelung constant, CsCl (B2):', a_this, a_this - a_ref) 

106 assert abs(a_this - a_ref) < 1e-13 

107 

108 if 1: # Zincblende, cubic ZnS 

109 cell = np.array([[.0, .5, .5], [.5, .0, .5], [.5, .5, .0]]) 

110 basis = np.array([[0., 0., 0.], [.25, .25, .25]]) 

111 r = np.array([0.0, 0.0, 0.0]) 

112 charges = np.array([1, -1]) 

113 ewald = Ewald(cell) 

114 e = ewald.get_electrostatic_potential(r, basis, charges, 

115 excludefroml0=0) 

116 a_this = -0.25 * np.sqrt(3) * e 

117 a_ref = 1.63805505338879 

118 if verbose: 

119 print('Madelung constant, Zincblende (B3):', 

120 a_this, a_this - a_ref) 

121 assert abs(a_this - a_ref) < 1e-13 

122 

123 if 1: # Ideal Wurtzite, ZnS, (B4) 

124 a = 1. 

125 c = np.sqrt(8. / 3.) * a 

126 u = 3. / 8. 

127 cell = np.array([[.5 * a, -0.5 * np.sqrt(3) * a, 0], 

128 [.5 * a, 0.5 * np.sqrt(3) * a, 0], [0., 0., c]]) 

129 basis = np.array([[.5 * a, .5 / np.sqrt(3) * a, 0.], 

130 [.5 * a, -.5 / np.sqrt(3) * a, 0.5 * c], 

131 [.5 * a, .5 / np.sqrt(3) * a, u * c], 

132 [.5 * a, -.5 / np.sqrt(3) * a, (.5 + u) * c]]) 

133 r = np.array([.5 * a, .5 / np.sqrt(3) * a, 0.]) 

134 charges = np.array([1., 1., -1, -1]) 

135 ewald = Ewald(cell) 

136 e = ewald.get_electrostatic_potential(r, basis, charges, 

137 excludefroml0=0) 

138 a_this = -1. * u * c * e 

139 a_ref = 1.64132162737 

140 if verbose: 

141 print('Madelung constant, Ideal Wurtzite (B4):', a_this, 

142 a_this - a_ref) 

143 assert abs(a_this - a_ref) < 1e-11