Coverage for gpaw/response/coulomb_kernels.py: 89%
132 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1"""This module defines Coulomb and XC kernels for the response model.
2"""
3import warnings
4import numpy as np
5from ase.dft import monkhorst_pack
6from gpaw.response.qpd import SingleQPWDescriptor
7from gpaw.response.dyson import PWKernel
8from gpaw.kpt_descriptor import to1bz
11class NewCoulombKernel(PWKernel):
12 def __init__(self, Vbare_G):
13 self.Vbare_G = Vbare_G
15 @classmethod
16 def from_qpd(cls, qpd, **kwargs):
17 return cls(get_coulomb_kernel(qpd, **kwargs))
19 def get_number_of_plane_waves(self):
20 return len(self.Vbare_G)
22 def _add_to(self, x_GG):
23 x_GG.flat[::self.nG + 1] += self.Vbare_G
26class CoulombKernel:
27 def __init__(self, truncation, *, N_c, pbc_c, kd):
28 self.truncation = truncation
29 assert self.truncation in {None, '0D', '2D'}
30 self.N_c = N_c
31 self.pbc_c = pbc_c
32 self.kd = kd
34 @classmethod
35 def from_gs(cls, gs, *, truncation):
36 return cls(truncation,
37 N_c=gs.kd.N_c,
38 pbc_c=gs.atoms.get_pbc(),
39 kd=gs.kd)
41 def new(self, *, truncation):
42 return CoulombKernel(truncation,
43 N_c=self.N_c,
44 pbc_c=self.pbc_c,
45 kd=self.kd)
47 def description(self):
48 if self.truncation is None:
49 return 'No Coulomb truncation'
50 else:
51 return f'{self.truncation} Coulomb truncation'
53 def sqrtV(self, qpd, q_v=None):
54 return self.V(qpd, q_v=q_v)**0.5
56 def V(self, qpd, q_v=None):
57 assert isinstance(qpd, SingleQPWDescriptor)
58 return get_coulomb_kernel(qpd, self.N_c, pbc_c=self.pbc_c, q_v=q_v,
59 truncation=self.truncation)
61 def kernel(self, qpd, q_v=None):
62 return np.diag(self.V(qpd, q_v=q_v))
64 def integrated_kernel(self, qpd, reduced, tofirstbz=False, *, N):
65 return get_integrated_kernel(
66 qpd=qpd, N_c=self.N_c, pbc_c=self.pbc_c,
67 truncation=self.truncation, reduced=reduced,
68 tofirstbz=tofirstbz, N=N)
71def get_coulomb_kernel(qpd, N_c, q_v=None, truncation=None, *, pbc_c)\
72 -> np.ndarray:
73 """Factory function that calls the specified flavour
74 of the Coulomb interaction"""
76 if truncation is None:
77 qG_Gv = qpd.get_reciprocal_vectors(add_q=True)
78 if q_v is not None:
79 assert qpd.optical_limit
80 qG_Gv += q_v
81 if qpd.optical_limit and q_v is None:
82 v_G = np.zeros(len(qpd.G2_qG[0]))
83 v_G[0] = 4 * np.pi
84 v_G[1:] = 4 * np.pi / (qG_Gv[1:]**2).sum(axis=1)
85 else:
86 v_G = 4 * np.pi / (qG_Gv**2).sum(axis=1)
88 elif truncation == '2D':
89 v_G = calculate_2D_truncated_coulomb(qpd, pbc_c=pbc_c, q_v=q_v)
90 if qpd.optical_limit and q_v is None:
91 v_G[0] = 0.0
93 elif truncation == '0D':
94 from gpaw.hybrids.wstc import WignerSeitzTruncatedCoulomb
95 wstc = WignerSeitzTruncatedCoulomb(qpd.gd.cell_cv, np.ones(3, int))
96 v_G = wstc.get_potential(qpd)
98 elif truncation in {'1D'}:
99 raise feature_removed()
101 else:
102 raise ValueError('Truncation scheme %s not implemented' % truncation)
104 return v_G.astype(complex)
107def calculate_2D_truncated_coulomb(qpd, q_v=None, *, pbc_c):
108 """ Simple 2D truncation of Coulomb kernel PRB 73, 205119.
110 Note: q_v is only added to qG_Gv if qpd.optical_limit is True.
112 """
114 qG_Gv = qpd.get_reciprocal_vectors(add_q=True)
115 if qpd.optical_limit:
116 if q_v is not None:
117 qG_Gv += q_v
118 else: # only to avoid warning. Later set to zero in factory function
119 qG_Gv[0] = [1., 1., 1.]
120 else:
121 assert q_v is None
123 # The non-periodic direction is determined from pbc_c
124 Nn_c = np.where(~pbc_c)[0]
125 Np_c = np.where(pbc_c)[0]
126 assert len(Nn_c) == 1
127 assert len(Np_c) == 2
128 # Truncation length is half of cell vector in non-periodic direction
129 R = qpd.gd.cell_cv[Nn_c[0], Nn_c[0]] / 2.
131 qGp_G = ((qG_Gv[:, Np_c[0]])**2 + (qG_Gv[:, Np_c[1]]**2))**0.5
132 qGn_G = qG_Gv[:, Nn_c[0]]
134 v_G = 4 * np.pi / (qG_Gv**2).sum(axis=1)
135 if np.allclose(qGn_G[0], 0) or qpd.optical_limit:
136 """sin(qGn_G * R) = 0 when R = L/2 and q_n = 0.0"""
137 v_G *= 1.0 - np.exp(-qGp_G * R) * np.cos(qGn_G * R)
138 else:
139 """Normal component of q is not zero"""
140 a_G = qGn_G / qGp_G * np.sin(qGn_G * R) - np.cos(qGn_G * R)
141 v_G *= 1. + np.exp(-qGp_G * R) * a_G
143 return v_G.astype(complex)
146def get_integrated_kernel(qpd, N_c, truncation=None,
147 reduced=False, tofirstbz=False, *, pbc_c, N):
148 from scipy.special import j1, k0, j0, k1 # type: ignore
149 # ignore type hints for the above import
150 B_cv = 2 * np.pi * qpd.gd.icell_cv
151 Nf_c = np.array([N, N, N])
152 if reduced:
153 # Only integrate periodic directions if truncation is used
154 Nf_c[np.where(~pbc_c)[0]] = 1
156 q_qc = monkhorst_pack(Nf_c)
158 if tofirstbz:
159 # We make a 1st BZ shaped integration volume, by reducing the full
160 # Monkhorts-Pack grid to the 1st BZ.
161 q_qc = to1bz(q_qc, qpd.gd.cell_cv)
163 q_qc /= N_c
164 q_qc += qpd.q_c
166 if tofirstbz:
167 # Because we added the q_c q-point vector, the integration volume is
168 # no longer strictly inside the 1st BZ of the full cell. Thus, we
169 # reduce it again.
170 q_qc = to1bz(q_qc, qpd.gd.cell_cv)
172 q_qv = np.dot(q_qc, B_cv)
174 Nn_c = np.where(~pbc_c)[0]
175 Np_c = np.where(pbc_c)[0]
176 assert len(Nn_c) + len(Np_c) == 3
178 if truncation is None:
179 V_q = 4 * np.pi / np.sum(q_qv**2, axis=1)
180 if len(Np_c) < 3:
181 warnings.warn(f"You should be using truncation='{len(Np_c)}D'")
182 elif truncation == '2D':
183 assert len(Np_c) == 2
185 # Truncation length is half of cell vector in non-periodic direction
186 R = qpd.gd.cell_cv[Nn_c[0], Nn_c[0]] / 2.
188 qp_q = ((q_qv[:, Np_c[0]])**2 + (q_qv[:, Np_c[1]]**2))**0.5
189 qn_q = q_qv[:, Nn_c[0]]
191 V_q = 4 * np.pi / (q_qv**2).sum(axis=1)
192 a_q = qn_q / qp_q * np.sin(qn_q * R) - np.cos(qn_q * R)
193 V_q *= 1. + np.exp(-qp_q * R) * a_q
194 elif truncation == '1D':
195 assert len(Np_c) == 1
196 # The radius is determined from area of non-periodic part of cell
197 Acell_cv = qpd.gd.cell_cv[Nn_c, :][:, Nn_c]
198 R = abs(np.linalg.det(Acell_cv) / np.pi)**0.5
200 qnR_q = (q_qv[:, Nn_c[0]]**2 + q_qv[:, Nn_c[1]]**2)**0.5 * R
201 qpR_q = abs(q_qv[:, Np_c[0]]) * R
202 V_q = 4 * np.pi / (q_qv**2).sum(axis=1)
203 V_q *= (1.0 + qnR_q * j1(qnR_q) * k0(qpR_q)
204 - qpR_q * j0(qnR_q) * k1(qpR_q))
205 elif truncation == '0D':
206 assert len(Np_c) == 0
207 R = (3 * qpd.gd.volume / (4 * np.pi))**(1. / 3.)
208 q2_q = (q_qv**2).sum(axis=1)
209 V_q = 4 * np.pi / q2_q
210 V_q *= 1.0 - np.cos(q2_q**0.5 * R)
212 return np.sum(V_q) / len(V_q), np.sum(V_q**0.5) / len(V_q)
215def feature_removed():
216 return RuntimeError(
217 '0D and 1D truncation have been removed due to not being tested. '
218 'If you need them, please find them in '
219 'ec9e49e25613bb99cd69eec9d2613e38b9f6e6e1 '
220 'and make sure to add tests in order to have them re-added.')