1from math import sqrt, pi
2
3import pytest
4import numpy as np
5from ase import Atoms
6from ase.units import Bohr, Hartree
7
8from gpaw.mpi import rank, size
9from gpaw import GPAW
10from gpaw.external import ConstantElectricField
11from gpaw.utilities import packed_index
12from gpaw.pair_density import PairDensity
13
14# Three ways to compute the polarizability of hydrogen:
15# 1. Perturbation theory
16# 2. Constant electric field
17# 3. Middle of an electric dipole
18
19# Note: The analytical value for the polarizability
20# is 4.5 a0**3 (e.g. PRA 33, 3671), while the experimental
21# value is 4.6 a0**3 (e.g. PR 133, A629).
22
23
24@pytest.mark.skip(reason='TODO')
25def test_stark_shift():
26    from gpaw.point_charges import PointCharges
27    to_au = Hartree / Bohr**2
28    to_eVA = Hartree / Bohr
29
30    def dipole_op(c, state1, state2, k=0, s=0):
31        # Taken from KSSingle, maybe make this accessible in
32        # KSSingle?
33        wfs = c.wfs
34        gd = wfs.gd
35
36        kpt = None
37        for i in wfs.kpt_u:
38            if i.k == k and i.s == s:
39                kpt = i
40
41        pd = PairDensity(c)
42        pd.initialize(kpt, state1, state2)
43
44        # coarse grid contribution
45        # <i|r|j> is the negative of the dipole moment (because of negative
46        # e- charge)
47        me = -gd.calculate_dipole_moment(pd.get())
48
49        # augmentation contributions
50        ma = np.zeros(me.shape)
51        pos_av = c.get_atoms().get_positions() / Bohr
52        for a, P_ni in kpt.P_ani.items():
53            Ra = pos_av[a]
54            Pi_i = P_ni[state1]
55            Pj_i = P_ni[state2]
56            Delta_pL = wfs.setups[a].Delta_pL
57            ni = len(Pi_i)
58
59            ma0 = 0
60            ma1 = np.zeros(me.shape)
61            for i in range(ni):
62                for j in range(ni):
63                    pij = Pi_i[i] * Pj_i[j]
64                    ij = packed_index(i, j, ni)
65                    # L=0 term
66                    ma0 += Delta_pL[ij, 0] * pij
67                    # L=1 terms
68                    if wfs.setups[a].lmax >= 1:
69                        # see spherical_harmonics.py for
70                        # L=1:y L=2:z; L=3:x
71                        ma1 += np.array([
72                            Delta_pL[ij, 3], Delta_pL[ij, 1], Delta_pL[ij, 2]
73                        ]) * pij
74            ma += sqrt(4 * pi / 3) * ma1 + Ra * sqrt(4 * pi) * ma0
75        gd.comm.sum(ma)
76
77        me += ma
78
79        return me * Bohr
80
81    # Currently only works on a single processor
82    assert size == 1
83
84    maxfield = 0.01
85    nfs = 5  # number of field
86    nbands = 30  # number of bands
87    h = 0.20  # grid spacing
88
89    debug = not False
90
91    if debug:
92        txt = 'gpaw.out'
93    else:
94        txt = None
95
96    test1 = True
97    test2 = True
98    test3 = True
99
100    a0 = 6.0
101    a = Atoms('H', positions=[[a0 / 2, a0 / 2, a0 / 2]], cell=[a0, a0, a0])
102
103    alpha1 = None
104    alpha2 = None
105    alpha3 = None
106
107    # Test 1
108
109    if test1:
110        c = GPAW(h=h,
111                 nbands=nbands + 10,
112                 spinpol=True,
113                 hund=True,
114                 xc='LDA',
115                 eigensolver='cg',
116                 convergence={
117                     'bands': nbands,
118                     'eigenstates': 3.3e-4
119                 },
120                 maxiter=1000,
121                 txt=txt)
122        a.calc = c
123        a.get_potential_energy()
124
125        o1 = c.get_occupation_numbers(spin=0)
126
127        if o1[0] > 0.0:
128            spin = 0
129        else:
130            spin = 1
131
132        alpha = 0.0
133
134        ev = c.get_eigenvalues(0, spin)
135        for i in range(1, nbands):
136            mu_x, mu_y, mu_z = dipole_op(c, 0, i, k=0, s=spin)
137
138            alpha += mu_z**2 / (ev[i] - ev[0])
139
140        alpha *= 2
141
142        if rank == 0 and debug:
143            print('From perturbation theory:')
144            print('  alpha = ', alpha, ' A**2/eV')
145            print('  alpha = ', alpha * to_au, ' Bohr**3')
146
147        alpha1 = alpha
148
149    ###
150
151    c = GPAW(
152        h=h,
153        nbands=2,
154        spinpol=True,
155        hund=True,
156        xc='LDA',
157        # eigensolver  = 'cg',
158        convergence={
159            'bands': nbands,
160            'eigenstates': 3.3e-4
161        },
162        maxiter=1000,
163        txt=txt)
164    a.calc = c
165
166    # Test 2
167
168    if test2:
169        e = []
170        e1s = []
171        d = []
172        fields = np.linspace(-maxfield, maxfield, nfs)
173        for field in fields:
174            if rank == 0 and debug:
175                print(field)
176            c.set(external=ConstantElectricField(field))
177            etot = a.get_potential_energy()
178            e += [etot]
179            ev0 = c.get_eigenvalues(0)
180            ev1 = c.get_eigenvalues(1)
181            e1s += [min(ev0[0], ev1[0])]
182            dip = c.get_dipole_moment()
183            d += [dip[2]]
184
185        pol1, dummy = np.polyfit(fields, d, 1)
186        pol2, dummy1, dummy2 = np.polyfit(fields, e1s, 2)
187
188        if rank == 0 and debug:
189            print('From shift in 1s-state at constant electric field:')
190            print('  alpha = ', -pol2, ' A**2/eV')
191            print('  alpha = ', -pol2 * to_au, ' Bohr**3')
192
193            print('From dipole moment at constant electric field:')
194            print('  alpha = ', pol1, ' A**2/eV')
195            print('  alpha = ', pol1 * to_au, ' Bohr**3')
196
197            np.savetxt('ecf.out', np.transpose([fields, e, e1s, d]))
198
199        assert abs(pol1 + pol2) < 0.0001
200
201        alpha2 = (pol1 - pol2) / 2
202
203    # Test 3
204
205    if test3:
206        pcd = 1000.0  # distance of the two point charges
207        maxcharge = 100.0  # maximum charge on the point charge
208
209        e = []
210        e1s = []
211        d = []
212        charges = np.linspace(-maxcharge, maxcharge, nfs)
213        fields = []
214        for charge in charges:
215            ex = PointCharges(positions=[[a0 / 2, a0 / 2, -pcd / 2 + a0 / 2],
216                                         [a0 / 2, a0 / 2, pcd / 2 + a0 / 2]],
217                              charges=[charge, -charge])
218            c.set(external=ex)
219            etot = a.get_potential_energy()
220            e += [etot]
221            ev0 = c.get_eigenvalues(0)
222            ev1 = c.get_eigenvalues(1)
223            e1s += [min(ev0[0], ev1[0])]
224            dip = c.get_dipole_moment()
225            d += [dip[2]]
226            field = ex.get_taylor(position=a[0].position)[1][1]
227            if rank == 0 and debug:
228                print(field * to_eVA,
229                      2 * charge / ((pcd / 2)**2) * Hartree * Bohr)
230            fields += [field * to_eVA]
231
232        pol1, dummy = np.polyfit(fields, d, 1)
233        pol2, dummy1, dummy2 = np.polyfit(fields, e1s, 2)
234
235        if rank == 0 and debug:
236            print('From shift in 1s-state between two point charges:')
237            print('  alpha = ', -pol2, ' A**2/eV')
238            print('  alpha = ', -pol2 * to_au, ' Bohr**3')
239
240            # print 'From dipole moment between two point charges:'
241            # print '  alpha = ', pol1, ' A**2/eV'
242            # print '  alpha = ', pol1*to_au, ' Bohr**3'
243
244            np.savetxt('epc.out', np.transpose([fields, e, e1s, d]))
245
246        alpha3 = alpha
247
248    # This is a very, very rough test
249    assert abs(alpha1 - alpha2) < 0.01
250    assert abs(alpha3 - alpha2) < 0.01
251