1from math import sqrt, pi
2
3import numpy as np
4import pytest
5
6from gpaw.setup import create_setup
7from gpaw.grid_descriptor import GridDescriptor
8from gpaw.lfc import LFC
9from gpaw.xc import XC
10
11n = 60  # 40 /8 * 10
12a = 10.0
13
14
15@pytest.mark.ci
16def test_multipole():
17    gd = GridDescriptor((n, n, n), (a, a, a))
18    c_LL = np.identity(9, float)
19    a_Lg = gd.zeros(9)
20    xc = XC('LDA')
21    for soft in [False]:
22        s = create_setup('Cu', xc, lmax=2)
23        ghat_l = s.ghat_l
24        ghat_Lg = LFC(gd, [ghat_l])
25        ghat_Lg.set_positions([(0.54321, 0.5432, 0.543)])
26        a_Lg[:] = 0.0
27        ghat_Lg.add(a_Lg, {0: c_LL} if ghat_Lg.my_atom_indices else {})
28        for l in range(3):
29            for m in range(2 * l + 1):
30                L = l**2 + m
31                a_g = a_Lg[L]
32                Q0 = gd.integrate(a_g) / sqrt(4 * pi)
33                Q1_m = -gd.calculate_dipole_moment(a_g) / sqrt(4 * pi / 3)
34                print(Q0)
35                if l == 0:
36                    Q0 -= 1.0
37                    Q1_m[:] = 0.0
38                elif l == 1:
39                    Q1_m[(m + 1) % 3] -= 1.0
40                print(soft, l, m, Q0, Q1_m)
41                assert abs(Q0) < 2e-6
42                assert np.alltrue(abs(Q1_m) < 3e-5)
43        b_Lg = np.reshape(a_Lg, (9, -1))
44        S_LL = np.inner(b_Lg, b_Lg)
45        gd.comm.sum(S_LL)
46        S_LL.ravel()[::10] = 0.0
47        print(max(abs(S_LL).ravel()))
48        assert max(abs(S_LL).ravel()) < 3e-4
49