Coverage for gpaw/lrtddft/apmb.py: 90%
286 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-12 00:18 +0000
1"""Omega matrix for functionals with Hartree-Fock exchange.
3"""
4from math import sqrt
6from ase.units import Hartree
7from ase.utils.timing import Timer
8from ase.utils import IOContext
9import numpy as np
10from numpy.linalg import inv
11from scipy.linalg import eigh
13from gpaw import debug
14import gpaw.mpi as mpi
15from gpaw.lrtddft.omega_matrix import OmegaMatrix
16from gpaw.pair_density import PairDensity
17from gpaw.helmholtz import HelmholtzSolver
18from gpaw.utilities.blas import mmm
21class ApmB(OmegaMatrix):
23 """Omega matrix for functionals with Hartree-Fock exchange.
25 """
27 def get_full(self):
29 hybrid = ((self.xc is not None) and
30 hasattr(self.xc, 'hybrid') and
31 (self.xc.hybrid > 0.0))
32 if self.fullkss.npspins < 2 and hybrid:
33 raise RuntimeError('Does not work spin-unpolarized ' +
34 'with hybrids (use nspins=2)')
36 if hasattr(self.xc, 'rsf') and (self.xc.rsf == 'Yukawa'):
37 self.screened_poissonsolver = HelmholtzSolver(
38 k2=-self.xc.omega**2, eps=1e-11, nn=3)
39 self.screened_poissonsolver.set_grid_descriptor(self.gd)
40 self.paw.timer.start('ApmB RPA')
41 self.ApB = self.Om
42 self.AmB = self.get_rpa()
43 self.paw.timer.stop()
45 if self.xc is not None:
46 self.paw.timer.start('ApmB XC')
47 self.get_xc() # inherited from OmegaMatrix
48 self.paw.timer.stop()
50 def get_rpa(self):
51 """Calculate RPA and Hartree-fock part of the A+-B matrices."""
53 # shorthands
54 kss = self.fullkss
55 finegrid = self.finegrid
56 yukawa = hasattr(self.xc, 'rsf') and (self.xc.rsf == 'Yukawa')
58 # calculate omega matrix
59 nij = len(kss)
60 self.log('RPAhyb', nij, 'transitions')
62 AmB = np.zeros((nij, nij))
63 ApB = self.ApB
65 # storage place for Coulomb integrals
66 integrals = {}
67 if yukawa:
68 rsf_integrals = {}
69 # setup things for IVOs
70 if (hasattr(self.xc, 'excitation') and
71 (self.xc.excitation is not None or self.xc.excited != 0)):
72 sin_tri_weight = 1
73 if self.xc.excitation is not None:
74 ex_type = self.xc.excitation.lower()
75 if ex_type == 'singlet':
76 sin_tri_weight = 2
77 elif ex_type == 'triplet':
78 sin_tri_weight = 0
79 homo = int(self.paw.get_number_of_electrons() // 2)
80 ivo_l = homo - self.xc.excited - 1
81 else:
82 ivo_l = None
84 for ij in range(nij):
85 self.log('RPAhyb kss[' + '%d' % ij + ']=', kss[ij])
87 timer = Timer()
88 timer.start('init')
89 timer2 = Timer()
91 # smooth density including compensation charges
92 timer2.start('with_compensation_charges 0')
93 rhot_p = kss[ij].with_compensation_charges(
94 finegrid != 0)
95 timer2.stop()
97 # integrate with 1/|r_1-r_2|
98 timer2.start('poisson')
99 phit_p = np.zeros(rhot_p.shape, rhot_p.dtype)
100 self.poisson.solve(phit_p, rhot_p, charge=None)
101 timer2.stop()
103 timer.stop()
104 t0 = timer.get_time('init')
105 timer.start(ij)
107 if finegrid == 1:
108 rhot = kss[ij].with_compensation_charges()
109 phit = self.gd.zeros()
110 self.restrict(phit_p, phit)
111 else:
112 phit = phit_p
113 rhot = rhot_p
115 for kq in range(ij, nij):
116 if kq != ij:
117 # smooth density including compensation charges
118 timer2.start('kq with_compensation_charges')
119 rhot = kss[kq].with_compensation_charges(
120 finegrid == 2)
121 timer2.stop()
122 pre = self.weight_Kijkq(ij, kq)
124 timer2.start('integrate')
125 I = self.Coulomb_integral_kss(kss[ij], kss[kq], phit, rhot)
126 if kss[ij].spin == kss[kq].spin:
127 name = self.Coulomb_integral_name(kss[ij].i, kss[ij].j,
128 kss[kq].i, kss[kq].j,
129 kss[ij].spin)
130 integrals[name] = I
131 ApB[ij, kq] = pre * I
132 timer2.stop()
134 if ij == kq:
135 epsij = kss[ij].get_energy() / kss[ij].get_weight()
136 AmB[ij, kq] += epsij
137 ApB[ij, kq] += epsij
139 timer.stop()
140# timer2.write()
141 if ij < (nij - 1):
142 self.log('RPAhyb estimated time left',
143 self.time_left(timer, t0, ij, nij))
145 # add HF parts and apply symmetry
146 if hasattr(self.xc, 'hybrid'):
147 weight = self.xc.hybrid
148 else:
149 weight = 0.0
150 for ij in range(nij):
151 self.log('HF kss[' + '%d' % ij + ']')
152 timer = Timer()
153 timer.start('init')
154 timer.stop()
155 t0 = timer.get_time('init')
156 timer.start(ij)
158 i = kss[ij].i
159 j = kss[ij].j
160 s = kss[ij].spin
161 for kq in range(ij, nij):
162 if kss[ij].pspin == kss[kq].pspin:
163 k = kss[kq].i
164 q = kss[kq].j
165 ikjq = self.Coulomb_integral_ijkq(i, k, j, q, s, integrals)
166 iqkj = self.Coulomb_integral_ijkq(i, q, k, j, s, integrals)
167 if yukawa: # Yukawa integrals might be caches
168 ikjq -= self.Coulomb_integral_ijkq(
169 i, k, j, q, s, rsf_integrals, yukawa)
170 iqkj -= self.Coulomb_integral_ijkq(
171 i, q, k, j, s, rsf_integrals, yukawa)
172 ApB[ij, kq] -= weight * (ikjq + iqkj)
173 AmB[ij, kq] -= weight * (ikjq - iqkj)
175 ApB[kq, ij] = ApB[ij, kq]
176 AmB[kq, ij] = AmB[ij, kq]
178 timer.stop()
179 if ij < (nij - 1):
180 self.log('HF estimated time left',
181 self.time_left(timer, t0, ij, nij))
183 if ivo_l is not None:
184 # IVO RPA after Berman, Kaldor, Chem. Phys. 43 (3) 1979
185 # doi: 10.1016/0301-0104(79)85205-2
186 l = ivo_l
187 for ij in range(nij):
188 i = kss[ij].i
189 j = kss[ij].j
190 s = kss[ij].spin
191 for kq in range(ij, nij):
192 if kss[kq].i == i and kss[ij].pspin == kss[kq].pspin:
193 k = kss[kq].i
194 q = kss[kq].j
195 jqll = self.Coulomb_integral_ijkq(j, q, l, l, s,
196 integrals)
197 jllq = self.Coulomb_integral_ijkq(j, l, l, q, s,
198 integrals)
199 if yukawa:
200 jqll -= self.Coulomb_integral_ijkq(j, q, l, l, s,
201 rsf_integrals,
202 yukawa)
203 jllq -= self.Coulomb_integral_ijkq(j, l, l, q, s,
204 rsf_integrals,
205 yukawa)
206 jllq *= sin_tri_weight
207 ApB[ij, kq] += weight * (jqll - jllq)
208 AmB[ij, kq] += weight * (jqll - jllq)
209 ApB[kq, ij] = ApB[ij, kq]
210 AmB[kq, ij] = AmB[ij, kq]
211 return AmB
213 def Coulomb_integral_name(self, i, j, k, l, spin):
214 """return a unique name considering the Coulomb integral
215 symmetry"""
216 def ij_name(i, j):
217 return str(max(i, j)) + ' ' + str(min(i, j))
219 # maximal gives the first
220 if max(i, j) >= max(k, l):
221 base = ij_name(i, j) + ' ' + ij_name(k, l)
222 else:
223 base = ij_name(k, l) + ' ' + ij_name(i, j)
224 return base + ' ' + str(spin)
226 def Coulomb_integral_ijkq(self, i, j, k, q, spin, integrals,
227 yukawa=False):
228 name = self.Coulomb_integral_name(i, j, k, q, spin)
229 if name in integrals:
230 return integrals[name]
232 # create the Kohn-Sham singles
233 kss_ij = PairDensity(self.paw)
234 kss_ij.initialize(self.paw.wfs.kpt_u[spin], i, j)
235 kss_kq = PairDensity(self.paw)
236 kss_kq.initialize(self.paw.wfs.kpt_u[spin], k, q)
238 rhot_p = kss_ij.with_compensation_charges(
239 self.finegrid != 0)
240 phit_p = np.zeros(rhot_p.shape, rhot_p.dtype)
241 if yukawa:
242 self.screened_poissonsolver.solve(phit_p, rhot_p, charge=None)
243 else:
244 self.poisson.solve(phit_p, rhot_p, charge=None)
246 if self.finegrid == 1:
247 phit = self.gd.zeros()
248 self.restrict(phit_p, phit)
249 else:
250 phit = phit_p
252 rhot = kss_kq.with_compensation_charges(
253 self.finegrid == 2)
255 integrals[name] = self.Coulomb_integral_kss(kss_ij, kss_kq,
256 phit, rhot,
257 yukawa=yukawa)
258 return integrals[name]
260 def timestring(self, t):
261 ti = int(t + .5)
262 td = int(ti // 86400)
263 st = ''
264 if td > 0:
265 st += '%d' % td + 'd'
266 ti -= td * 86400
267 th = int(ti // 3600)
268 if th > 0:
269 st += '%d' % th + 'h'
270 ti -= th * 3600
271 tm = int(ti // 60)
272 if tm > 0:
273 st += '%d' % tm + 'm'
274 ti -= tm * 60
275 st += '%d' % ti + 's'
276 return st
278 def mapAB(self, restrict={}):
279 """Map A+B, A-B matrices according to constraints."""
281 map, self.kss = self.get_map(restrict)
282 if map is None:
283 ApB = self.ApB.copy()
284 AmB = self.AmB.copy()
285 else:
286 nij = len(self.kss)
287 ApB = np.empty((nij, nij))
288 AmB = np.empty((nij, nij))
289 for ij in range(nij):
290 for kq in range(nij):
291 ApB[ij, kq] = self.ApB[map[ij], map[kq]]
292 AmB[ij, kq] = self.AmB[map[ij], map[kq]]
294 return ApB, AmB
296 def diagonalize(self, restrict={}, TDA=False):
297 """Evaluate Eigenvectors and Eigenvalues"""
299 ApB, AmB = self.mapAB(restrict)
300 nij = len(self.kss)
302 if TDA:
303 # Tamm-Dancoff approximation (B=0)
304 eigenvalues, evecs = eigh(0.5 * (ApB + AmB))
305 self.eigenvectors = evecs.T
306 self.eigenvalues = eigenvalues ** 2
307 else:
308 # the occupation matrix
309 C = np.empty((nij,))
310 for ij in range(nij):
311 C[ij] = 1. / self.kss[ij].fij
313 S = C * inv(AmB) * C
314 S = sqrt_matrix(inv(S).copy())
316 # get Omega matrix
317 M = np.zeros(ApB.shape)
318 mmm(1.0, S, 'N', ApB, 'N', 0.0, M)
319 self.eigenvectors = np.zeros(ApB.shape)
320 mmm(1.0, M, 'N', S, 'N', 0.0, self.eigenvectors)
322 self.eigenvalues, self.eigenvectors.T[:] = eigh(self.eigenvectors)
324 def read(self, filename=None, fh=None):
325 """Read myself from a file"""
326 if mpi.rank == 0:
327 with IOContext() as io:
328 if fh is None:
329 fd = io.openfile(filename, 'r', comm=self.paw.world)
330 else:
331 fd = fh
332 fd.readline()
333 nij = int(fd.readline())
334 ApB = np.zeros((nij, nij))
335 for ij in range(nij):
336 l = fd.readline().split()
337 for kq in range(ij, nij):
338 ApB[ij, kq] = float(l[kq - ij])
339 ApB[kq, ij] = ApB[ij, kq]
340 self.ApB = ApB
342 fd.readline()
343 nij = int(fd.readline())
344 AmB = np.zeros((nij, nij))
345 for ij in range(nij):
346 l = fd.readline().split()
347 for kq in range(ij, nij):
348 AmB[ij, kq] = float(l[kq - ij])
349 AmB[kq, ij] = AmB[ij, kq]
350 self.AmB = AmB
352 def weight_Kijkq(self, ij, kq):
353 """weight for the coupling matrix terms"""
354 return 2.
356 def write(self, filename=None, fh=None):
357 """Write current state to a file."""
358 if mpi.rank == 0:
359 with IOContext() as io:
360 if fh is None:
361 fd = io.openfile(filename, 'r', comm=self.paw.world)
362 else:
363 fd = fh
364 fd.write('# A+B\n')
365 nij = len(self.fullkss)
366 fd.write('%d\n' % nij)
367 for ij in range(nij):
368 for kq in range(ij, nij):
369 fd.write(' %g' % self.ApB[ij, kq])
370 fd.write('\n')
372 fd.write('# A-B\n')
373 nij = len(self.fullkss)
374 fd.write('%d\n' % nij)
375 for ij in range(nij):
376 for kq in range(ij, nij):
377 fd.write(' %g' % self.AmB[ij, kq])
378 fd.write('\n')
380 def __str__(self):
381 string = '<ApmB> '
382 if hasattr(self, 'eigenvalues'):
383 string += 'dimension ' + ('%d' % len(self.eigenvalues))
384 string += '\neigenvalues: '
385 for ev in self.eigenvalues:
386 string += ' ' + ('%f' % (sqrt(ev) * Hartree))
387 return string
390def sqrt_matrix(a, preserve=False):
391 """Get the sqrt of a symmetric matrix a (diagonalize is used).
392 The matrix is kept if preserve=True, a=sqrt(a) otherwise."""
393 n = len(a)
394 if debug:
395 assert a.flags.contiguous
396 assert a.dtype == float
397 assert a.shape == (n, n)
398 if preserve:
399 b = a.copy()
400 else:
401 b = a
403 # diagonalize to get the form b = Z * D * Z^T
404 # where D is diagonal
405 D = np.empty((n,))
406 D, b.T[:] = eigh(b, lower=True)
407 ZT = b.copy()
408 Z = np.transpose(b)
410 # c = Z * sqrt(D)
411 c = Z * np.sqrt(D)
413 # sqrt(b) = c * Z^T
414 mmm(1.0, np.ascontiguousarray(c), 'N', ZT, 'N', 0.0, b)
415 return b