Coverage for gpaw/response/q0_correction.py: 100%
74 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1import numpy as np
4class Q0Correction:
5 def __init__(self, cell_cv, bzk_kc, N_c):
6 self.cell_cv = cell_cv
7 self.bzk_kc = bzk_kc
8 self.N_c = N_c
10 # Check that basic assumptions of cell and k-grid
11 # for Q0Correction are fulfilled
12 assert N_c[2] == 1 # z-axis is non periodic direction
13 eps = 1e-14
14 assert (abs(cell_cv[:2, 2]).max() < eps and
15 abs(cell_cv[2, :2]).max() < eps and
16 cell_cv[2, 2] > 0)
18 # Hardcoded crap?
19 x0density = 0.1 # ? 0.01
21 # Generate numerical q-point grid
22 self.rcell_cv = 2 * np.pi * np.linalg.inv(cell_cv).T
24 # Prepare stuff
25 self.q0cell_cv = np.array([1, 1, 1])**0.5 * self.rcell_cv / self.N_c
26 L = cell_cv[2, 2]
27 q0density = 2. / L * x0density
28 npts_c = np.ceil(np.sum(self.q0cell_cv**2, axis=1)**0.5 /
29 q0density).astype(int)
30 npts_c[2] = 1
31 npts_c += (npts_c + 1) % 2
32 self.npts_c = npts_c
34 def add_q0_correction(self, qpd, W_GG, einv_GG,
35 chi0_xvG, chi0_vv, sqrtV_G):
36 from ase.dft import monkhorst_pack
37 qpts_qc = self.bzk_kc
38 pi = np.pi
39 L = self.cell_cv[2, 2]
41 vc_G0 = sqrtV_G[1:]**2
43 B_GG = einv_GG[1:, 1:]
44 u_v0G = vc_G0[np.newaxis, :]**0.5 * chi0_xvG[0, :, 1:]
45 u_vG0 = vc_G0[np.newaxis, :]**0.5 * chi0_xvG[1, :, 1:]
46 U_vv = -chi0_vv
47 a_v0G = -np.dot(u_v0G, B_GG)
48 a_vG0 = -np.dot(u_vG0, B_GG.T)
49 A_vv = U_vv - np.dot(a_v0G, u_vG0.T)
50 S_v0G = a_v0G
51 S_vG0 = a_vG0
52 L_vv = A_vv
54 # Get necessary G vectors.
55 G_Gv = qpd.get_reciprocal_vectors(add_q=False)[1:]
56 G_Gv += np.array([1e-14, 1e-14, 0])
57 G2_G = np.sum(G_Gv**2, axis=1)
58 Gpar_G = np.sum(G_Gv[:, 0:2]**2, axis=1)**0.5
60 # There is still a lot of stuff here,
61 # which could go to the constructor! XXX
62 iq = np.argmin(np.sum(qpts_qc**2, axis=1))
63 assert np.allclose(qpts_qc[iq], 0)
64 q0vol = abs(np.linalg.det(self.q0cell_cv))
66 qpts_qc = monkhorst_pack(self.npts_c)
67 qgamma = np.argmin(np.sum(qpts_qc**2, axis=1))
69 qpts_qv = np.dot(qpts_qc, self.q0cell_cv)
70 qpts_q = np.sum(qpts_qv**2, axis=1)**0.5
71 qpts_q[qgamma] = 1e-14
72 qdir_qv = qpts_qv / qpts_q[:, np.newaxis]
73 qdir_qvv = qdir_qv[:, :, np.newaxis] * qdir_qv[:, np.newaxis, :]
74 nq = len(qpts_qc)
75 q0area = q0vol / self.q0cell_cv[2, 2]
76 dq0 = q0area / nq
77 dq0rad = (dq0 / pi)**0.5
78 R = L / 2.
79 x0area = q0area * R**2
80 dx0rad = dq0rad * R
82 exp_q = 4 * pi * (1 - np.exp(-qpts_q * R))
83 dv_G = ((pi * L * G2_G * np.exp(-Gpar_G * R) * np.cos(G_Gv[:, 2] * R) -
84 4 * pi * Gpar_G * (1 - np.exp(-Gpar_G * R) *
85 np.cos(G_Gv[:, 2] * R))) /
86 (G2_G**1.5 * Gpar_G *
87 (4 * pi * (1 - np.exp(-Gpar_G * R) *
88 np.cos(G_Gv[:, 2] * R)))**0.5))
90 dv_Gv = dv_G[:, np.newaxis] * G_Gv
92 # Add corrections
93 W_GG[:, 0] = 0.0
94 W_GG[0, :] = 0.0
96 A_q = np.sum(qdir_qv * np.dot(qdir_qv, L_vv), axis=1)
97 frac_q = 1. / (1 + exp_q * A_q)
99 # HEAD:
100 w00_q = -(exp_q / qpts_q)**2 * A_q * frac_q
101 w00_q[qgamma] = 0.0
102 W_GG[0, 0] = w00_q.sum() / nq
103 Axy = 0.5 * (L_vv[0, 0] + L_vv[1, 1]) # in-plane average
104 a0 = 4 * pi * Axy + 1
106 W_GG[0, 0] += -((a0 * dx0rad - np.log(a0 * dx0rad + 1)) /
107 a0**2 / x0area * 2 * np.pi * (2 * pi * L)**2 * Axy)
109 # WINGS:
110 u_q = -exp_q / qpts_q * frac_q
111 W_GG[1:, 0] = 1. / nq * np.dot(
112 np.sum(qdir_qv * u_q[:, np.newaxis], axis=0),
113 S_vG0 * sqrtV_G[np.newaxis, 1:])
115 W_GG[0, 1:] = 1. / nq * np.dot(
116 np.sum(qdir_qv * u_q[:, np.newaxis], axis=0),
117 S_v0G * sqrtV_G[np.newaxis, 1:])
119 # BODY:
120 # Constant corrections:
121 W_GG[1:, 1:] += 1. / nq * sqrtV_G[1:, None] * sqrtV_G[None, 1:] * \
122 np.tensordot(S_v0G, np.dot(S_vG0.T,
123 np.sum(-qdir_qvv *
124 exp_q[:, None, None] *
125 frac_q[:, None, None],
126 axis=0)), axes=(0, 1))
127 u_vvv = np.tensordot(u_q[:, None] * qpts_qv, qdir_qvv, axes=(0, 0))
128 # Gradient corrections:
129 W_GG[1:, 1:] += 1. / nq * np.sum(
130 dv_Gv[:, :, None] * np.tensordot(
131 S_v0G, np.tensordot(u_vvv, S_vG0 * sqrtV_G[None, 1:],
132 axes=(2, 0)), axes=(0, 1)), axis=1)