Coverage for gpaw/coulomb.py: 50%
206 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1from math import pi, sqrt
3from ase.units import Hartree
4import numpy as np
5from numpy.fft import fftn
7from gpaw.lfc import LocalizedFunctionsCollection as LFC
8from gpaw.pair_density import PairDensity2 as PairDensity
9from gpaw.poisson import PoissonSolver, FFTPoissonSolver
10from gpaw.utilities import packed_index, unpack_hermitian, unpack_density
11from gpaw.utilities.ewald import madelung
12from gpaw.utilities.tools import construct_reciprocal, tri2full, symmetrize
13from gpaw.utilities.gauss import Gaussian
14from gpaw.utilities.blas import r2k
17def get_vxc(paw, spin=0, U=None):
18 """Calculate matrix elements of the xc-potential."""
19 assert not paw.hamiltonian.xc.xcfunc.orbital_dependent, "LDA/GGA's only"
20 assert paw.wfs.dtype == float, 'Complex waves not implemented'
22 if U is not None: # Rotate xc matrix
23 return np.dot(U.T.conj(), np.dot(get_vxc(paw, spin), U))
25 gd = paw.hamiltonian.gd
26 psit_nG = paw.wfs.kpt_u[spin].psit_nG[:]
27 if paw.density.nt_sg is None:
28 paw.density.interpolate_pseudo_density()
29 nt_g = paw.density.nt_sg[spin]
30 vxct_g = paw.density.finegd.zeros()
31 paw.hamiltonian.xc.get_energy_and_potential(nt_g, vxct_g)
32 vxct_G = gd.empty()
33 paw.hamiltonian.restrict(vxct_g, vxct_G)
34 Vxc_nn = np.zeros((paw.wfs.bd.nbands, paw.wfs.bd.nbands))
36 # Apply pseudo part
37 r2k(.5 * gd.dv, psit_nG, vxct_G * psit_nG, .0, Vxc_nn) # lower triangle
38 tri2full(Vxc_nn, 'L') # Fill in upper triangle from lower
39 gd.comm.sum(Vxc_nn)
41 # Add atomic PAW corrections
42 for a, P_ni in paw.wfs.kpt_u[spin].P_ani.items():
43 D_sp = paw.density.D_asp[a][:]
44 H_sp = np.zeros_like(D_sp)
45 paw.wfs.setups[a].xc_correction.calculate_energy_and_derivatives(
46 D_sp, H_sp)
47 H_ii = unpack_hermitian(H_sp[spin])
48 Vxc_nn += np.dot(P_ni, np.dot(H_ii, P_ni.T))
49 return Vxc_nn * Hartree
52class Coulomb:
53 """Class used to evaluate two index coulomb integrals."""
54 def __init__(self, gd, poisson=None):
55 """Class should be initialized with a grid_descriptor 'gd' from
56 the gpaw module.
57 """
58 self.gd = gd
59 self.poisson = poisson
61 def load(self, method):
62 """Make sure all necessary attributes have been initialized"""
64 assert method in ('real', 'recip_gauss', 'recip_ewald'), \
65 str(method) + ' is an invalid method name,\n' +\
66 'use either real, recip_gauss, or recip_ewald'
68 if method.startswith('recip'):
69 if self.gd.comm.size > 1:
70 raise RuntimeError("Cannot do parallel FFT, use method='real'")
71 if not hasattr(self, 'k2'):
72 self.k2, self.N3 = construct_reciprocal(self.gd)
73 if method.endswith('ewald') and not hasattr(self, 'ewald'):
74 # cutoff radius
75 assert self.gd.orthogonal
76 rc = 0.5 * np.average(self.gd.cell_cv.diagonal())
77 # ewald potential: 1 - cos(k rc)
78 self.ewald = (np.ones(self.gd.n_c) -
79 np.cos(np.sqrt(self.k2) * rc))
80 # lim k -> 0 ewald / k2
81 self.ewald[0, 0, 0] = 0.5 * rc**2
82 elif method.endswith('gauss') and not hasattr(self, 'ng'):
83 gauss = Gaussian(self.gd)
84 self.ng = gauss.get_gauss(0) / sqrt(4 * pi)
85 self.vg = gauss.get_gauss_pot(0) / sqrt(4 * pi)
86 else: # method == 'real'
87 if not hasattr(self, 'solve'):
88 if self.poisson is not None:
89 self.solve = self.poisson.solve
90 else:
91 solver = PoissonSolver('fd', nn=2, eps=1e-12)
92 solver.set_grid_descriptor(self.gd)
93 # solver.initialize()
94 self.solve = solver.solve
96 def coulomb(self, n1, n2=None, Z1=None, Z2=None, method='recip_gauss'):
97 """Evaluates the coulomb integral of n1 and n2
99 The coulomb integral is defined by::
101 *
102 / / n1(r) n2(r')
103 (n1 | n2) = | dr | dr' -------------,
104 / / |r - r'|
106 where n1 and n2 could be complex.
108 real:
109 Evaluate directly in real space using gaussians to neutralize
110 density n2, such that the potential can be generated by standard
111 procedures
113 recip_ewald:
114 Evaluate by Fourier transform.
115 Divergence at division by k^2 is avoided by utilizing the Ewald /
116 Tuckermann trick, which formally requires the densities to be
117 localized within half of the unit cell.
119 recip_gauss:
120 Evaluate by Fourier transform.
121 Divergence at division by k^2 is avoided by removing total charge
122 of n1 and n2 with gaussian density ng::
124 * * *
125 (n1|n2) = (n1 - Z1 ng|n2 - Z2 ng) + (Z2 n1 + Z1 n2 - Z1 Z2 ng | ng)
127 The evaluation of the integral (n1 - Z1 ng|n2 - Z2 ng) is done in
128 k-space using FFT techniques.
129 """
130 self.load(method)
131 # determine integrand using specified method
132 if method == 'real':
133 I = self.gd.zeros()
134 if n2 is None:
135 n2 = n1
136 Z2 = Z1
137 self.solve(I, n2, charge=Z2, zero_initial_phi=True)
138 I += madelung(self.gd.cell_cv) * self.gd.integrate(n2)
139 I *= n1.conj()
140 elif method == 'recip_ewald':
141 n1k = fftn(n1)
142 if n2 is None:
143 n2k = n1k
144 else:
145 n2k = fftn(n2)
146 I = n1k.conj() * n2k * self.ewald * 4 * pi / (self.k2 * self.N3)
147 else: # method == 'recip_gauss':
148 # Determine total charges
149 if Z1 is None:
150 Z1 = self.gd.integrate(n1)
151 if Z2 is None and n2 is not None:
152 Z2 = self.gd.integrate(n2)
154 # Determine the integrand of the neutral system
155 # (n1 - Z1 ng)* int dr' (n2 - Z2 ng) / |r - r'|
156 nk1 = fftn(n1 - Z1 * self.ng)
157 if n2 is None:
158 I = abs(nk1)**2 * 4 * pi / (self.k2 * self.N3)
159 else:
160 nk2 = fftn(n2 - Z2 * self.ng)
161 I = nk1.conj() * nk2 * 4 * pi / (self.k2 * self.N3)
163 # add the corrections to the integrand due to neutralization
164 if n2 is None:
165 I += (2 * np.real(np.conj(Z1) * n1) -
166 abs(Z1)**2 * self.ng) * self.vg
167 else:
168 I += (np.conj(Z1) * n2 + Z2 * n1.conj() -
169 np.conj(Z1) * Z2 * self.ng) * self.vg
170 if n1.dtype == float and (n2 is None or n2.dtype == float):
171 return np.real(self.gd.integrate(I))
172 else:
173 return self.gd.integrate(I)
176class CoulombNEW:
177 def __init__(self, gd, setups, spos_ac, fft=False):
178 assert gd.comm.size == 1
179 self.rhot1_G = gd.empty()
180 self.rhot2_G = gd.empty()
181 self.pot_G = gd.empty()
182 self.dv = gd.dv
183 if fft:
184 self.poisson = FFTPoissonSolver()
185 else:
186 self.poisson = PoissonSolver(name='fd', nn=3, eps=1e-12)
187 self.poisson.set_grid_descriptor(gd)
188 self.setups = setups
190 # Set coarse ghat
191 self.Ghat = LFC(gd, [setup.ghat_l for setup in setups],
192 integral=sqrt(4 * pi))
193 self.Ghat.set_positions(spos_ac)
195 def calculate(self, nt1_G, nt2_G, P1_ap, P2_ap):
196 I = 0.0
197 self.rhot1_G[:] = nt1_G
198 self.rhot2_G[:] = nt2_G
200 Q1_aL = {}
201 Q2_aL = {}
202 for a, P1_p in P1_ap.items():
203 P2_p = P2_ap[a]
204 setup = self.setups[a]
206 # Add atomic corrections to integral
207 I += 2 * np.dot(P1_p, np.dot(setup.M_pp, P2_p))
209 # Add compensation charges to pseudo densities
210 Q1_aL[a] = np.dot(P1_p, setup.Delta_pL)
211 Q2_aL[a] = np.dot(P2_p, setup.Delta_pL)
212 self.Ghat.add(self.rhot1_G, Q1_aL)
213 self.Ghat.add(self.rhot2_G, Q2_aL)
215 # Add coulomb energy of compensated pseudo densities to integral
216 self.poisson.solve(self.pot_G, self.rhot2_G, charge=None,
217 zero_initial_phi=True)
218 I += np.vdot(self.rhot1_G, self.pot_G) * self.dv
220 return I * Hartree
223class HF:
224 def __init__(self, paw):
225 paw.initialize_positions()
226 self.nspins = paw.wfs.nspins
227 self.nbands = paw.wfs.bd.nbands
228 self.restrict = paw.hamiltonian.restrict
229 self.pair_density = PairDensity(paw.density, paw.spos_ac,
230 finegrid=True)
231 self.dv = paw.wfs.gd.dv
232 self.dtype = paw.wfs.dtype
233 self.setups = paw.wfs.setups
235 # Allocate space for matrices
236 self.nt_G = paw.wfs.gd.empty()
237 self.rhot_g = paw.density.finegd.empty()
238 self.vt_G = paw.wfs.gd.empty()
239 self.vt_g = paw.density.finegd.empty()
240 self.poisson_solve = paw.hamiltonian.poisson.solve
242 def apply(self, paw, u=0):
243 H_nn = np.zeros((self.nbands, self.nbands), self.dtype)
244 self.soft_pseudo(paw, H_nn, u=u)
245 self.atomic_val_val(paw, H_nn, u=u)
246 self.atomic_val_core(paw, H_nn, u=u)
247 return H_nn * Hartree
249 def soft_pseudo(self, paw, H_nn, h_nn=None, u=0):
250 if h_nn is None:
251 h_nn = H_nn
252 kpt = paw.wfs.kpt_u[u]
253 pd = self.pair_density
254 deg = 2 / self.nspins
255 fmin = 1e-9
256 Htpsit_nG = np.zeros(kpt.psit_nG.shape, self.dtype)
258 for n1 in range(self.nbands):
259 psit1_G = kpt.psit_nG[n1]
260 f1 = kpt.f_n[n1] / deg
261 for n2 in range(n1, self.nbands):
262 psit2_G = kpt.psit_nG[n2]
263 f2 = kpt.f_n[n2] / deg
264 if f1 < fmin and f2 < fmin:
265 continue
267 pd.initialize(kpt, n1, n2)
268 pd.get_coarse(self.nt_G)
269 pd.add_compensation_charges(self.nt_G, self.rhot_g)
270 self.poisson_solve(self.vt_g, -self.rhot_g,
271 charge=-float(n1 == n2),
272 zero_initial_phi=True)
273 self.restrict(self.vt_g, self.vt_G)
274 Htpsit_nG[n1] += f2 * self.vt_G * psit2_G
275 if n1 != n2:
276 Htpsit_nG[n2] += f1 * self.vt_G * psit1_G
278 v_aL = paw.density.ghat.dict()
279 paw.density.ghat.integrate(self.vt_g, v_aL)
280 for a, v_L in v_aL.items():
281 v_ii = unpack_hermitian(
282 np.dot(paw.wfs.setups[a].Delta_pL, v_L))
283 P_ni = kpt.P_ani[a]
284 h_nn[:, n1] += f2 * np.dot(P_ni, np.dot(v_ii, P_ni[n2]))
285 if n1 != n2:
286 h_nn[:, n2] += f1 * np.dot(P_ni,
287 np.dot(v_ii, P_ni[n1]))
289 symmetrize(h_nn) # Grrrr why!!! XXX
291 # Fill in lower triangle
292 r2k(0.5 * self.dv, kpt.psit_nG[:], Htpsit_nG, 1.0, H_nn)
294 # Fill in upper triangle from lower
295 tri2full(H_nn, 'L')
297 def atomic_val_val(self, paw, H_nn, u=0):
298 deg = 2 / self.nspins
299 kpt = paw.wfs.kpt_u[u]
300 for a, P_ni in kpt.P_ani.items():
301 # Add atomic corrections to the valence-valence exchange energy
302 # --
303 # > D C D
304 # -- ii iiii ii
305 setup = paw.wfs.setups[a]
306 D_p = paw.density.D_asp[a][kpt.s]
307 H_p = np.zeros_like(D_p)
308 D_ii = unpack_density(D_p)
309 ni = len(D_ii)
310 for i1 in range(ni):
311 for i2 in range(ni):
312 A = 0.0
313 for i3 in range(ni):
314 p13 = packed_index(i1, i3, ni)
315 for i4 in range(ni):
316 p24 = packed_index(i2, i4, ni)
317 A += setup.M_pp[p13, p24] * D_ii[i3, i4]
318 p12 = packed_index(i1, i2, ni)
319 H_p[p12] -= 2 / deg * A / ((i1 != i2) + 1)
320 H_nn += np.dot(P_ni, np.inner(unpack_hermitian(H_p), P_ni.conj()))
322 def atomic_val_core(self, paw, H_nn, u=0):
323 kpt = paw.wfs.kpt_u[u]
324 for a, P_ni in kpt.P_ani.items():
325 dH_ii = unpack_hermitian(-paw.wfs.setups[a].X_p)
326 H_nn += np.dot(P_ni, np.inner(dH_ii, P_ni.conj()))