1import pytest
2import numpy as np
3from ase.units import Bohr
4from ase.build import bulk
5from gpaw import GPAW, FermiDirac
6from gpaw.response.df import DielectricFunction
7from gpaw.test import equal, findpeak
8
9
10@pytest.mark.libxc
11def test_response_diamond_absorption(in_tmp_dir):
12    a = 6.75 * Bohr
13    atoms = bulk('C', 'diamond', a=a)
14
15    calc = GPAW(mode='pw',
16                kpts=(3, 3, 3),
17                eigensolver='rmm-diis',
18                occupations=FermiDirac(0.001))
19
20    atoms.calc = calc
21    atoms.get_potential_energy()
22    calc.write('C.gpw', 'all')
23
24    eM1_ = 9.727
25    eM2_ = 9.548
26    w0_ = 10.7782
27    I0_ = 5.47
28    w_ = 10.7532
29    I_ = 5.98
30
31    # Macroscopic dielectric constant calculation
32    df = DielectricFunction('C.gpw', frequencies=(0.,), eta=0.001, ecut=50,
33                            hilbert=False)
34    eM1, eM2 = df.get_macroscopic_dielectric_constant()
35    equal(eM1, eM1_, 0.01)
36    equal(eM2, eM2_, 0.01)
37
38    # Absorption spectrum calculation RPA
39    df = DielectricFunction('C.gpw',
40                            eta=0.25,
41                            ecut=50,
42                            frequencies=np.linspace(0, 24., 241),
43                            hilbert=False)
44    a0, a = df.get_dielectric_function(filename=None)
45    df.check_sum_rule(a.imag)
46
47    equal(a0[0].real, eM1_, 0.01)
48    equal(a[0].real, eM2_, 0.01)
49    w, I = findpeak(np.linspace(0, 24., 241), a0.imag)
50    equal(w, w0_, 0.01)
51    equal(I / (4 * np.pi), I0_, 0.1)
52    w, I = findpeak(np.linspace(0, 24., 241), a.imag)
53    equal(w, w_, 0.01)
54    equal(I / (4 * np.pi), I_, 0.1)
55
56    a0, a = df.get_polarizability(filename=None)
57
58    w, I = findpeak(np.linspace(0, 24., 241), a0.imag)
59    equal(w, w0_, 0.01)
60    equal(I, I0_, 0.01)
61    w, I = findpeak(np.linspace(0, 24., 241), a.imag)
62    equal(w, w_, 0.01)
63    equal(I, I_, 0.01)
64
65    # Absorption spectrum calculation ALDA
66    w0_ = 10.7931
67    I0_ = 5.36
68    w_ = 10.7562
69    I_ = 5.8803
70
71    a0, a = df.get_polarizability(filename=None, xc='ALDA')
72
73    w, I = findpeak(np.linspace(0, 24., 241), a0.imag)
74    equal(w, w0_, 0.01)
75    equal(I, I0_, 0.1)
76    w, I = findpeak(np.linspace(0, 24., 241), a.imag)
77    equal(w, w_, 0.01)
78    equal(I, I_, 0.1)
79
80    # Absorption spectrum calculation long-range kernel
81    w0_ = 10.2189
82    I0_ = 5.14
83    w_ = 10.2906
84    I_ = 5.6955
85
86    a0, a = df.get_polarizability(filename=None, xc='LR0.25')
87
88    w, I = findpeak(np.linspace(0, 24., 241), a0.imag)
89    equal(w, w0_, 0.01)
90    equal(I, I0_, 0.1)
91    w, I = findpeak(np.linspace(0, 24., 241), a.imag)
92    equal(w, w_, 0.01)
93    equal(I, I_, 0.1)
94
95    # Absorption spectrum calculation Bootstrap
96    w0_ = 10.37
97    I0_ = 5.27
98    w_ = 10.4600
99    I_ = 6.0263
100
101    a0, a = df.get_polarizability(filename=None, xc='Bootstrap')
102
103    w, I = findpeak(np.linspace(0, 24., 241), a0.imag)
104    equal(w, w0_, 0.02)
105    equal(I, I0_, 0.2)
106    w, I = findpeak(np.linspace(0, 24., 241), a.imag)
107    equal(w, w_, 0.02)
108    equal(I, I_, 0.2)
109
110    # Absorption spectrum calculation RPA - Wigner-Seitz truncation
111    w0_ = 10.7780
112    I0_ = 5.5472
113    w_ = 10.7532
114    I_ = 6.0750
115
116    df_ws = DielectricFunction('C.gpw',
117                               eta=0.25,
118                               ecut=50,
119                               frequencies=np.linspace(0, 24., 241),
120                               hilbert=False,
121                               truncation='wigner-seitz')
122
123    a0_ws, a_ws = df_ws.get_polarizability(filename=None)
124
125    w, I = findpeak(np.linspace(0, 24., 241), a0_ws.imag)
126    equal(w, w0_, 0.02)
127    equal(I, I0_, 0.2)
128    w, I = findpeak(np.linspace(0, 24., 241), a_ws.imag)
129    equal(w, w_, 0.1)
130    equal(I, I_, 0.2)
131    # The Wigner-Seitz truncation does not give exactly the
132    # same for kpts=(3,3,3)
133