Coverage for gpaw/lrtddft2/ks_singles.py: 94%
211 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import os
3import io
5import numpy as np
7import ase.units
9import gpaw.mpi
11from gpaw.utilities.tools import coordinates
12from gpaw.utilities.tools import pick
13from gpaw.utilities import pack_density
15from gpaw.fd_operators import Gradient
18class KohnShamSingleExcitation:
19 def __init__(self, gs_calc, kpt_ind, occ_ind, unocc_ind):
20 self.calc = gs_calc
21 self.kpt_ind = kpt_ind
22 self.spin_ind = self.calc.wfs.kpt_u[self.kpt_ind].s
23 self.occ_ind = occ_ind
24 self.unocc_ind = unocc_ind
25 self.energy_diff = None
26 self.pop_diff = None
27 self.dip_mom_r = None
28 self.magn_mom = None
29 self.pair_density = None
31 # Calculates pair density of (i,p) transition (and given kpt)
32 # dnt_Gp: pair density without compensation charges on coarse grid
33 # dnt_gp: pair density without compensation charges on fine grid (XC)
34 # drhot_gp: pair density with compensation charges on fine grid (poisson)
35 def calculate_pair_density(self, dnt_Gip, dnt_gip, drhot_gip):
36 if self.pair_density is None:
37 # FIXME: USE EXISTING PairDensity method
38 self.pair_density = LRiPairDensity(self.calc.density)
39 self.pair_density.initialize(self.calc.wfs.kpt_u[self.kpt_ind],
40 self.occ_ind, self.unocc_ind)
41 self.pair_density.get(dnt_Gip, dnt_gip, drhot_gip)
43 def __str__(self):
44 if self.dip_mom_r is not None and self.dip_mom_v is not None\
45 and self.magn_mom is not None:
46 str = "# KS single excitation from state %05d to state %05d:"\
47 " dE_pi = %18.12lf, f_pi = %18.12lf,"\
48 " dmr_ip = (%18.12lf, %18.12lf, %18.12lf),"\
49 " dmv_ip = (%18.12lf, %18.12lf, %18.12lf),"\
50 " dmm_ip = %18.12lf" \
51 % (self.occ_ind,
52 self.unocc_ind,
53 self.energy_diff,
54 self.pop_diff,
55 self.dip_mom_r[0], self.dip_mom_r[1], self.dip_mom_r[2],
56 self.dip_mom_v[0], self.dip_mom_v[1], self.dip_mom_v[2],
57 self.magn_mom)
58 elif self.energy_diff is not None and self.pop_diff is not None:
59 str = "# KS single excitation from state %05d to state %05d:"\
60 " dE_pi = %18.12lf, f_pi = %18.12lf" \
61 % (self.occ_ind,
62 self.unocc_ind,
63 self.energy_diff,
64 self.pop_diff)
65 elif self.occ_ind is not None and self.unocc_ind is not None:
66 str = "# KS single excitation from state %05d to state %05d" \
67 % (self.occ_ind, self.unocc_ind)
68 else:
69 raise RuntimeError("Uninitialized KSSingle")
70 return str
73class KohnShamSingles:
74 def __init__(self, basefilename, gs_calc, kpt_index, min_occ, max_occ,
75 min_unocc, max_unocc, max_energy_diff, min_pop_diff, lr_comms,
76 txt):
77 self.basefilename = basefilename
78 self.calc = gs_calc
79 self.min_occ = min_occ
80 self.max_occ = max_occ
81 self.min_unocc = min_unocc
82 self.max_unocc = max_unocc
83 self.max_energy_diff = max_energy_diff
84 self.min_pop_diff = min_pop_diff
85 self.kpt_ind = kpt_index
86 self.lr_comms = lr_comms
87 self.txt = txt
89 self.kss_list = None
90 self.kss_list_ready = False
92 self.kss_prop = None
93 self.kss_prop_ready = False
95 def update_list(self):
96 # shorthands
97 eps_n = self.calc.wfs.kpt_u[self.kpt_ind].eps_n # eigen energies
98 f_n = self.calc.wfs.kpt_u[self.kpt_ind].f_n # occupations
100 # Create Kohn-Sham single excitation list with energy filter
101 old_kss_list = self.kss_list # save old list for later
102 self.kss_list = [] # create a completely new list
103 # Occupied loop
104 for i in range(self.min_occ, self.max_occ + 1):
105 # Unoccupied loop
106 for p in range(self.min_unocc, self.max_unocc + 1):
107 deps_pi = eps_n[p] - eps_n[i] # energy diff
108 df_ip = f_n[i] - f_n[p] # population diff
109 # Filter it
110 if (np.abs(deps_pi) <= self.max_energy_diff
111 and df_ip > self.min_pop_diff):
112 # i, p, deps, df, mur, muv, magn
113 kss = KohnShamSingleExcitation(self.calc, self.kpt_ind, i,
114 p)
115 kss.energy_diff = deps_pi
116 kss.pop_diff = df_ip
117 self.kss_list.append(kss)
119 # Sort by energy diff:
120 self.kss_list = sorted(self.kss_list, key=lambda x: x.energy_diff)
122 # Remove old transitions and add new
123 # BUT only add to the end of the list (otherwise lower triangle
124 # matrix is not filled completely)
125 if old_kss_list is not None:
126 new_kss_list = self.kss_list # required list
127 self.kss_list = [] # final list with correct order
128 # Old list first
129 # If in new and old lists
130 for kss_o in old_kss_list:
131 found = False
132 for kss_n in new_kss_list:
133 if (kss_o.occ_ind == kss_n.occ_ind
134 and kss_o.unocc_ind == kss_n.unocc_ind):
135 found = True
136 break
137 if found:
138 self.kss_list.append(kss_o) # Found, add to final list
139 else:
140 pass # else drop
142 # Now old transitions which are not in new list where dropped
144 # If only in new list
145 app_list = []
146 for kss_n in new_kss_list:
147 found = False
148 for kss in self.kss_list:
149 if (kss.occ_ind == kss_n.occ_ind
150 and kss.unocc_ind == kss_n.unocc_ind):
151 found = True
152 break
153 if not found:
154 app_list.append(kss_n) # Not found, add to final list
155 else:
156 pass # else skip to avoid duplicates
158 # Create the final list
159 self.kss_list += app_list
161 self.txt.write('Number of electron-hole pairs: %d\n' %
162 len(self.kss_list))
163 self.txt.write('Maximum energy difference: %6.3lf\n' %
164 (self.kss_list[-1].energy_diff * ase.units.Hartree))
166 # Prevent repeated work
167 self.kss_list_ready = True
169 def read(self):
170 # Read KS_singles file if exists
171 # occ_index | unocc index | energy diff | population diff |
172 # dmx | dmy | dmz | magnx | magny | magnz
173 kss_file = self.basefilename + '.KS_singles'
174 if os.path.exists(kss_file) and os.path.isfile(kss_file):
175 self.kss_list = []
177 # root reads and broadcasts
178 # a bit ugly and slow but works
179 data = None
180 if self.lr_comms.parent_comm.rank == 0:
181 with open(kss_file, 'r', 1024 * 1024) as fd:
182 data = fd.read()
183 data = gpaw.mpi.broadcast_string(data,
184 root=0,
185 comm=self.lr_comms.parent_comm)
187 # just parse
188 for line in io.StringIO(data):
189 line = line.split()
190 i, p = int(line[0]), int(line[1])
191 ediff, fdiff = float(line[2]), float(line[3])
192 dm = np.array([float(line[4]), float(line[5]), float(line[6])])
193 mm = np.array([float(line[7]), float(line[8]), float(line[9])])
194 kss = KohnShamSingleExcitation(self.calc, self.kpt_ind, i, p)
195 kss.energy_diff = ediff
196 kss.pop_diff = fdiff
197 kss.dip_mom_r = dm
198 kss.magn_mom = mm
199 # assert self.index_of_kss(i,p) is None, 'KS transition'\
200 # '%d->%d found twice in KS_singles files.' % (i,p)
201 self.kss_list.append(kss)
203 # if none read
204 if len(self.kss_list) <= 0:
205 self.kss_list = None
207 # Calculate dipole moment and magnetic moment for noninteracting
208 # Kohn-Sham transitions
209 # FIXME: Do it in parallel, currently doing repeated work
210 # Should distribute kss_list over eh_comm
211 def calculate(self):
212 # Check that kss_list is up-to-date
213 if not self.kss_list_ready:
214 self.update_list()
215 self.kss_prop_ready = False
217 # Check if we already have properties for all singles
218 if self.kss_prop_ready:
219 return
220 self.kss_prop_ready = True
221 for kss_ip in self.kss_list:
222 if kss_ip.dip_mom_r is None or kss_ip.magn_mom is None:
223 self.kss_prop_ready = False
224 break
225 # If had, done
226 if self.kss_prop_ready:
227 return
229 # self.timer.start('Calculate KS properties')
231 # Init pair densities
232 dnt_Gip = self.calc.wfs.gd.empty()
233 dnt_gip = self.calc.density.finegd.empty()
234 drhot_gip = self.calc.density.finegd.empty()
236 # Init gradients of wfs
237 grad_psit2_G = [
238 self.calc.wfs.gd.empty(),
239 self.calc.wfs.gd.empty(),
240 self.calc.wfs.gd.empty()
241 ]
243 # Init gradient operators
244 grad = []
245 dtype = pick(self.calc.wfs.kpt_u[self.kpt_ind].psit_nG, 0).dtype
246 for c in range(3):
247 grad.append(Gradient(self.calc.wfs.gd, c, dtype=dtype, n=2))
249 # Coordinate vector r
250 R0 = 0.5 * np.diag(self.calc.wfs.gd.cell_cv)
251 # R0 = np.array([0.0,0.0,0.0]) # old_version
252 # print 'R0', R0
253 r_cG, r2_G = coordinates(self.calc.wfs.gd, origin=R0)
254 r_cg, r2_g = coordinates(self.calc.density.finegd, origin=R0)
256 # Loop over all KS single excitations
257 #
258 # Transition dipole moment, mu_ip = <p| (-e r) |i>
259 # Magnetic transition dipole, m_ip = -(1/2c) <i|L|p>
260 # For total m_0I = -m_I0 = -(m_0I)^*, but not for m_ip (right?)
261 # R_0I = Im[mu_0I * m_0I]
262 # mu_ip^0I = omega_0I^(-1/2) D_ip S^(-1/2) F_0I
263 # m_ip^0I = -omega_0I^(+1/2) M_ip C^-1 S^(+1/2) F_0I
264 #
265 # S_ip,ip = - (eps_p - eps_i) / (n_p - n_i) (note: n_p < n_i)
266 # C_ip,ip = 1 / (n_p - n_i)
267 #
268 # See:
269 # WIREs Comput Mol Sci 2012, 2: 150-166 doi: 10.1002/wcms.55
270 # J. Chem. Phys., Vol. 116, 6930 (2002)
271 for kss_ip in self.kss_list:
272 # If have dipole moment and magnetic moment, already done and skip
273 if (kss_ip.dip_mom_r is not None and kss_ip.magn_mom is not None):
274 continue
276 # Transition dipole moment, mu_ip = <p| (-r) |i>
277 kss_ip.calculate_pair_density(dnt_Gip, dnt_gip, drhot_gip)
278 # kss_ip.dip_mom_r
279 # = self.calc.density.finegd.calculate_dipole_moment(drhot_gip)
280 # kss_ip.dip_mom_r
281 # = self.calc.density.finegd.calculate_dipole_moment(drhot_gip)
282 kss_ip.dip_mom_r = np.zeros(3)
283 kss_ip.dip_mom_r[0] = -self.calc.density.finegd.integrate(
284 r_cg[0] * drhot_gip)
285 kss_ip.dip_mom_r[1] = -self.calc.density.finegd.integrate(
286 r_cg[1] * drhot_gip)
287 kss_ip.dip_mom_r[2] = -self.calc.density.finegd.integrate(
288 r_cg[2] * drhot_gip)
290 # Magnetic transition dipole,
291 # m_ip = -(1/2c) <i|L|p> = i/2c <i|r x p|p>
292 # see Autschbach et al., J. Chem. Phys., 116, 6930 (2002)
294 # Gradients
295 for c in range(3):
296 grad[c].apply(kss_ip.pair_density.psit2_G, grad_psit2_G[c],
297 self.calc.wfs.kpt_u[self.kpt_ind].phase_cd)
299 # <psi1|r x grad|psi2>
300 # i j k
301 # x y z = (y pz - z py)i + (z px - x pz)j + (x py - y px)
302 # px py pz
303 rxnabla_g = np.zeros(3)
304 rxnabla_g[0] = self.calc.wfs.gd.integrate(
305 kss_ip.pair_density.psit1_G *
306 (r_cG[1] * grad_psit2_G[2] - r_cG[2] * grad_psit2_G[1]))
307 rxnabla_g[1] = self.calc.wfs.gd.integrate(
308 kss_ip.pair_density.psit1_G *
309 (r_cG[2] * grad_psit2_G[0] - r_cG[0] * grad_psit2_G[2]))
310 rxnabla_g[2] = self.calc.wfs.gd.integrate(
311 kss_ip.pair_density.psit1_G *
312 (r_cG[0] * grad_psit2_G[1] - r_cG[1] * grad_psit2_G[0]))
314 # augmentation contributions to magnetic moment
315 # <psi1| r x nabla |psi2> = <psi1| (r-Ra+Ra) x nabla |psi2>
316 # = <psi1| (r-Ra) x nabla |psi2> + Ra x <psi1| nabla |psi2>
317 rxnabla_a = np.zeros(3)
318 # <psi1| (r-Ra) x nabla |psi2>
319 for a, P_ni in self.calc.wfs.kpt_u[self.kpt_ind].P_ani.items():
320 Pi_i = P_ni[kss_ip.occ_ind]
321 Pp_i = P_ni[kss_ip.unocc_ind]
322 rxnabla_iiv = self.calc.wfs.setups[a].rxnabla_iiv
323 for c in range(3):
324 for i1, Pi in enumerate(Pi_i):
325 for i2, Pp in enumerate(Pp_i):
326 rxnabla_a[c] += Pi * Pp * rxnabla_iiv[i1, i2, c]
328 self.calc.wfs.gd.comm.sum(rxnabla_a) # sum up from different procs
330 # Ra x <psi1| nabla |psi2>
331 Rxnabla_a = np.zeros(3)
332 for a, P_ni in self.calc.wfs.kpt_u[self.kpt_ind].P_ani.items():
333 Pi_i = P_ni[kss_ip.occ_ind]
334 Pp_i = P_ni[kss_ip.unocc_ind]
335 nabla_iiv = self.calc.wfs.setups[a].nabla_iiv
336 Ra = (self.calc.atoms[a].position / ase.units.Bohr) - R0
337 for i1, Pi in enumerate(Pi_i):
338 for i2, Pp in enumerate(Pp_i):
339 # (y pz - z py)i + (z px - x pz)j + (x py - y px)k
340 Rxnabla_a[0] += Pi * Pp * (
341 Ra[1] * nabla_iiv[i1, i2, 2] -
342 Ra[2] * nabla_iiv[i1, i2, 1])
343 Rxnabla_a[1] += Pi * Pp * (
344 Ra[2] * nabla_iiv[i1, i2, 0] -
345 Ra[0] * nabla_iiv[i1, i2, 2])
346 Rxnabla_a[2] += Pi * Pp * (
347 Ra[0] * nabla_iiv[i1, i2, 1] -
348 Ra[1] * nabla_iiv[i1, i2, 0])
350 self.calc.wfs.gd.comm.sum(Rxnabla_a) # sum up from different procs
352 # print (kss_ip.occ_ind, kss_ip.unocc_ind),
353 # kss_ip.dip_mom_r, rxnabla_g, rxnabla_a, Rxnabla_a
355 # m_ip = -1/2c <i|r x p|p> = i/2c <i|r x nabla|p>
356 # just imaginary part!!!
357 kss_ip.magn_mom = ase.units.alpha / 2. * (rxnabla_g + rxnabla_a +
358 Rxnabla_a)
360 # Wait... to avoid io problems, and write KS_singles file
361 self.lr_comms.parent_comm.barrier()
362 if self.lr_comms.parent_comm.rank == 0:
363 self.kss_file = open(self.basefilename + '.KS_singles', 'w')
364 for kss_ip in self.kss_list:
365 format = '%08d %08d %18.12lf %18.12lf '
366 format += '%18.12lf %18.12lf %18.12lf '
367 format += '%18.12lf %18.12lf %18.12lf\n'
368 self.kss_file.write(
369 format %
370 (kss_ip.occ_ind, kss_ip.unocc_ind, kss_ip.energy_diff,
371 kss_ip.pop_diff, kss_ip.dip_mom_r[0], kss_ip.dip_mom_r[1],
372 kss_ip.dip_mom_r[2], kss_ip.magn_mom[0],
373 kss_ip.magn_mom[1], kss_ip.magn_mom[2]))
374 self.kss_file.close()
375 self.lr_comms.parent_comm.barrier()
377 self.kss_prop_ready = True # avoid repeated work
378 # self.timer.stop('Calculate KS properties')
381# Small utility class
383# FIXME: USE EXISTING METHOD
384class LRiPairDensity:
385 """Pair density calculator class."""
387 def __init__(self, density):
388 """Initialization needs density instance"""
389 self.density = density
391 def initialize(self, kpt, n1, n2):
392 """Set wave function indices."""
393 self.n1 = n1
394 self.n2 = n2
395 self.spin = kpt.s
396 self.P_ani = kpt.P_ani
397 self.psit1_G = pick(kpt.psit_nG, n1)
398 self.psit2_G = pick(kpt.psit_nG, n2)
400 def get(self, nt_G, nt_g, rhot_g):
401 """Get pair densities.
403 nt_G
404 Pair density without compensation charges on the coarse grid
406 nt_g
407 Pair density without compensation charges on the fine grid
409 rhot_g
410 Pair density with compensation charges on the fine grid
411 """
412 # Coarse grid product
413 np.multiply(self.psit1_G.conj(), self.psit2_G, nt_G)
414 # Interpolate to fine grid
415 self.density.interpolator.apply(nt_G, nt_g)
417 # Determine the compensation charges for each nucleus
418 Q_aL = {}
419 for a, P_ni in self.P_ani.items():
420 assert P_ni.dtype == float
421 # Generate density matrix
422 P1_i = P_ni[self.n1]
423 P2_i = P_ni[self.n2]
424 D_ii = np.outer(P1_i.conj(), P2_i)
425 # Pack with pack_density as it is used in the scalar product with
426 # the symmetric array Delta_pL
427 D_p = pack_density(D_ii)
428 # FIXME: CHECK THIS D_p = pack_density(D_ii, tolerance=1e30)
430 # Determine compensation charge coefficients:
431 Q_aL[a] = np.dot(D_p, self.density.setups[a].Delta_pL)
433 # Add compensation charges
434 rhot_g[:] = nt_g[:]
435 self.density.ghat.add(rhot_g, Q_aL)
436 # print 'dens', self.density.finegd.integrate(rhot_g)