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
144