Coverage for gpaw/cdft/cdft_coupling.py: 56%
532 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1'''A class for computing the
2parameters for Marcus theory
3from two constrained DFT
4wave functions
6Computes:
7-the coupling matrix Hab
8 Hab = <Psi_a|H|Psi_b>
9-reorganization energy lambda
10lambda = E_a(Rb)-E_a(Ra)
11'''
13import numpy as np
14from gpaw.cdft.cdft import (WeightFunc, get_ks_energy_wo_external,
15 get_all_weight_functions)
16from ase.units import kB as kb
17from gpaw.utilities.ps2ae import PS2AE, interpolate_weight
18from ase.units import Bohr
19from gpaw.mpi import rank
20import warnings
23spin_state_error = ('The cDFT wave functions have\n' +
24 'different spin states! Similar\n' +
25 'spin states are required for coupling constant\n' +
26 'calculation!')
28ae_ibz_error = ('The all electron calculation is unreliable with kpts\n' +
29 'Please set AE to false')
31nab_missing_error = ('The pair density matrix n_ab must\n' +
32 'be provided for weight matrix calculation')
33migliore_warning = ('WARNING! Migliore coupling might be unreliable!:\n' +
34 '<F|GS> =<I|GS> and dE_if=0')
37class CouplingParameters:
38 def __init__(self,
39 cdft_a=None,
40 cdft_b=None,
41 h=0.04,
42 calc_a=None,
43 calc_b=None,
44 calc_gs=None,
45 wfs_a='initial.gpw',
46 wfs_b='final.gpw',
47 AE=False,
48 charge_difference=False,
49 charge_regions_A=None,
50 spin_regions_A=None,
51 charge_regions_B=None,
52 spin_regions_B=None,
53 FA=None,
54 FB=None,
55 Va=None,
56 Vb=None,
57 NA=None,
58 NB=None,
59 n_ab=None,
60 VW_ab=None,
61 VW_ba=None,
62 w_k=None,
63 Rc={},
64 mu={},
65 specific_bands_A=None,
66 specific_bands_B=None,
67 kpoint=None,
68 spin=None,
69 S_matrix='S_matrix',
70 VW_matrix='VW_matrix',
71 S=None,
72 A=None,
73 B=None,
74 save_matrix=False,
75 freq=None,
76 temp=None,
77 E_tst=None,
78 reorg=None,
79 use_adiabatic_correction=False,
80 reaction_energy=0.,
81 band_occupation_cutoff=0.5,
82 energy_gap=None):
83 '''cdft_a cdft_b: cdft calculators
84 h = grid spacing in S_AB and W_AB calculations in AE mode
85 AE = Use all electron wave functions
88 if cdft calculators are not provided,
89 calc_a, calc_b, wfs_a, wfs_b,gd,
90 FA,FB,Va, Vb, NA, NB,
91 charge_regions_A/B must be provided
93 calc_a/b = GPAW calculators for a/b states
94 wfs_a/b = wfs of the a/b states
95 FA/FB = pseudo cdft free energies of A/B states in eV
96 as from cdft.free_energy() or cdft.Ecdft
98 Va/b = list of lagrange multipliers in eV
99 NA/NB = charge + spin constraints in A/B states
100 !NOTE! Total number of electrons as from cdft.get_constraints()
102 weight_A/B = the weight functions for A/B
103 Array
105 N_charge_regions_A/B = number of constrained charges
106 int
108 E_KS_A/B = KS energies with out the external potential
109 (Ekin+Epot+Ebar+Exc+S) in eV
111 n_ab = S_AB_ij(k) = <A_i(k)|B_j(k)>
113 VW_ab,VW_ba = <A(k)|sum_c Vc_B*Wc_B|B(k)>, <B(k)|sum_c Vc_A*Wc_A|B(k)>
115 WAB, WBA = sum_k sum_ij sum_c <A_i(k)|Wc_B|B_j(k)>
117 w_k = combined kpt weight, (w_kA*w_kB) in eq 48 and 49
118 can be provided to compute overlap and
119 weight matrices
121 specific_bands, kpoint, spin = [[],[]], int, int
122 allows to extract density matrix elements
123 and weight-density matrix elements for
124 specified orbitals Eq. 48 and 49
125 [[],[]] = indices of first for alpha,
126 second for beta spin bands
128 If one these is activated, matrix of
129 int[psi_sk^i*psi_sk^j] is returned instead
130 of S_ab or W_ab
132 S_matrix, W_matrix = str
133 place to print density and
134 weight-density matrix
135 Turned on automatically if bands,
136 spin, or kpoint are specified
138 S, A, B = np.array, stored migliore matrices
140 for calculation of electron transfer rates
141 using the Marcus theory the following can be
142 specified:
144 freq = effective frequency
145 temp = temperature
146 E_tst = adiabatic reaction barrier
147 use_adiabatic_correction = False/True
148 add the adiabadicity correction to
149 Marcus barrier (can some times give very
150 too small barriers for an adiabatic reaction
151 with large coupling...)
152 reaction_energy: E_B(RB)-E_A(RA)
153 band_occupation_cutoff = minimum filling of band
154 to be included in the wfs
155 Rc, mu: Gaussian parameters. If not provided either standard or
156 the ones from a cdft will be used
157 '''
159 self.charge_difference = charge_difference
160 self.AE = AE
162 if cdft_a is not None and cdft_b is not None:
163 self.cdft_A = cdft_a
164 self.cdft_B = cdft_b
165 self.Rc = self.cdft_A.Rc
166 self.mu = self.cdft_A.mu
168 self.calc_A = self.cdft_A.calc
169 self.calc_B = self.cdft_B.calc
171 if self.AE:
172 self.gd = self.cdft_A.get_grid()
173 else:
174 self.gd = self.cdft_A.calc.density.gd
176 # cDFT free energies
177 self.FA = self.cdft_A.cdft_free_energy()
178 self.FB = self.cdft_B.cdft_free_energy()
180 # cDFT energies
181 self.EA = self.cdft_A.dft_energy()
182 self.EB = self.cdft_B.dft_energy()
184 # lagrange multipliers
185 self.Va = self.cdft_A.get_lagrangians()
186 self.Vb = self.cdft_B.get_lagrangians()
188 # constraint values
189 self.NA = self.cdft_A.get_constraints()
190 self.NB = self.cdft_B.get_constraints()
192 # number of charge regions
193 self.n_charge_regionsA = self.cdft_A.n_charge_regions
194 self.n_charge_regionsB = self.cdft_B.n_charge_regions
195 self.regionsA = self.cdft_A.regions
196 self.regionsB = self.cdft_B.regions
197 self.atoms = self.calc_A.get_atoms()
199 if self.AE:
200 self.fineweightA = self.cdft_A.get_weight(save=False, pad=True)
201 self.fineweightB = self.cdft_B.get_weight(save=False, pad=True)
202 else:
203 self.coarseweightA = get_all_weight_functions(
204 self.calc_A.atoms, self.gd, self.regionsA,
205 self.charge_difference, self.Rc, self.mu)
206 self.coarseweightB = get_all_weight_functions(
207 self.calc_B.atoms, self.gd, self.regionsB,
208 self.charge_difference, self.Rc, self.mu)
210 else:
211 self.calc_A = calc_a
212 self.calc_B = calc_b
214 if self.AE:
215 self.gd = self.calc_A.density.finegd
216 else:
217 self.gd = self.calc_A.density.gd
219 # cDFT free energies
220 self.FA = FA
221 self.FB = FB
223 # Weights, i.e., lagrange multipliers
224 self.Va = Va
225 self.Vb = Vb
226 regionsA = []
227 if charge_regions_A is not None:
228 regionsA += charge_regions_A
229 if spin_regions_A is not None:
230 regionsA += spin_regions_A
232 regionsB = []
233 if charge_regions_B is not None:
234 regionsB += charge_regions_B
235 if spin_regions_B is not None:
236 regionsB += spin_regions_B
238 self.regionsA = regionsA
239 self.regionsB = regionsB
240 self.NA = NA
241 self.NB = NB
243 # see if a correct set of Rc and mu are provided
244 if (set(Rc) == set(self.calc_A.atoms.get_chemical_symbols()) and
245 set(mu) == set(self.calc_A.atoms.get_chemical_symbols())):
246 self.Rc = Rc
247 self.mu = mu
248 else:
249 # get mu and Rc the long way
250 w_temp = WeightFunc(gd=self.gd,
251 atoms=self.calc_A.atoms,
252 indices=self.regionsA,
253 Rc=Rc,
254 mu=mu)
255 self.Rc, self.mu = w_temp.get_Rc_and_mu()
256 del w_temp
258 if self.AE:
259 self.fineweightA = get_all_weight_functions(
260 self.calc_A.atoms, self.gd, self.regionsA,
261 self.charge_difference, self.Rc, self.mu)
263 self.fineweightB = get_all_weight_functions(
264 self.calc_B.atoms, self.gd, self.regionsB,
265 self.charge_difference, self.Rc, self.mu)
267 else:
268 self.coarseweightA = get_all_weight_functions(
269 self.calc_A.atoms, self.gd, self.regionsA,
270 self.charge_difference, self.Rc, self.mu)
271 self.coarseweightB = get_all_weight_functions(
272 self.calc_B.atoms, self.gd, self.regionsB,
273 self.charge_difference, self.Rc, self.mu)
275 self.n_charge_regionsA = len(charge_regions_A)
276 self.n_charge_regionsB = len(charge_regions_B)
278 # Do not include external potential
279 self.EA = get_ks_energy_wo_external(self.calc_A)
280 self.EB = get_ks_energy_wo_external(self.calc_B)
282 self.finegd = self.calc_A.density.finegd
283 self.energy_gap = energy_gap
284 self.Va = np.asarray(self.Va)
285 self.Vb = np.asarray(self.Vb)
287 # wave functions
288 if wfs_a is not None and wfs_b is not None:
289 self.wfs_A = wfs_a
290 self.wfs_B = wfs_b
291 else:
292 self.wfs_A = self.calc_A.wfs
293 self.wfs_B = self.calc_B.wfs
295 # ground state calculator for migliore
297 if calc_gs is not None:
298 self.calc_gs = calc_gs
299 if S is not None:
300 self.S = np.load(S)
301 if A is not None:
302 self.A = np.load(A)
303 if B is not None:
304 self.B = np.load(B)
305 if w_k is not None:
306 self.w_AB = self.w_AC = self.w_BC = np.load(w_k)
308 # use precalculated pair density/weight matrices
309 if w_k:
310 self.w_k = np.load(w_k)
311 if n_ab:
312 self.n_ab = np.load(n_ab)
314 if VW_ab:
315 self.VW_AB = np.load(VW_ab)
316 self.VW_BA = np.load(VW_ba)
318 # specify bands to be treated
319 self.specific_bands_A = specific_bands_A
320 self.specific_bands_B = specific_bands_B
322 if kpoint:
323 self.kpoints = kpoint
324 else:
325 self.kpoints = self.calc_A.wfs.kd.nibzkpts
326 self.spin = spin
328 # where to save matrices
329 self.S_matrix = S_matrix
330 self.VW_matrix = VW_matrix
332 self.save_matrix = save_matrix
333 # Only part of bands treated --> matrices saved
334 if self.specific_bands_A or self.specific_bands_B or self.spin:
335 self.save_matrix = True
337 self.H_AA = self.EA
338 self.H_BB = self.EB
340 self.density = self.calc_A.density.finegd.zeros()
341 self.atoms = self.calc_A.atoms
342 self.natoms = len(self.atoms)
343 self.h = h # all-electron interpolation grid spacing
345 # controls for rate calculation
347 self.E_tst = E_tst # adiabatic barrier
349 # effective frequency
350 if freq:
351 self.freq = freq
352 else:
353 self.freq = 1.e12
354 # temperature
355 if temp:
356 self.temp = temp
357 else:
358 self.temp = 298
359 self.reorg = reorg
361 self.reaction_energy = reaction_energy
362 self.use_adiabatic_correction = use_adiabatic_correction
363 self.band_occupation_cutoff = band_occupation_cutoff
365 def get_coupling_term(self):
366 # coupling by Migliore
367 # dx.doi.org/10.1021/ct200192d, eq 6
368 # check that needed terms exist
370 if (hasattr(self, 'H') and hasattr(self, 'S')):
371 if rank == 0:
372 H = self.H.copy()
373 S = self.S.copy()
375 else:
376 self.make_hamiltonian_matrix()
377 if rank == 0:
378 H = self.H.copy()
379 S = self.S.copy()
380 if rank == 0:
381 H_av = 0.5 * (H[0][1] + H[1][0])
382 S_av = 0.5 * (S[0][1] + S[1][0])
384 self.ct = (1. / (1 - S_av**2)) * np.abs(H_av - S_av *
385 (H[0][0] + H[1][1]) / 2.)
386 return self.ct
388 def get_coupling_term_from_lowdin(self):
389 self.make_hamiltonian_matrix()
390 if rank == 0:
391 H_orth = self.lowdin_orthogonalize_cdft_hamiltonian()
392 self.ct = 1. / 2. * np.real(H_orth[0][1] + H_orth[1][0])
393 return self.ct
395 def lowdin_orthogonalize_cdft_hamiltonian(self):
396 # H_orth = S^(-1/2)*H_cdft*S^(-1/2)
398 U, D, V = np.linalg.svd(self.S, full_matrices=True) # V=U(+)
399 s = np.diag(D**(-0.5))
400 S_square = np.dot(U, np.dot(s, V)) # S^(-1/2)
402 H_orth = np.dot(S_square, np.dot(self.H, S_square))
403 return H_orth
405 def get_migliore_coupling(self):
406 # the entire migliore method
407 # from dx.doi.org/10.1021/ct200192d,
408 # requires also ground state adiabatic wf
409 if not hasattr(self, 'S') or not hasattr(self, 'A') or not hasattr(
410 self, 'B'):
411 self.S, self.A, self.B, self.w_AB, self.w_AC, self.w_BC = (
412 self.get_migliore_wf_overlaps(
413 self.calc_A, self.calc_B, self.calc_gs))
414 if rank == 0:
415 SIF = 0.
416 SFI = 0.
417 A = 0.
418 B = 0.
420 for k in range(len(self.w_AB)):
421 SIF += self.w_AB[k] * np.linalg.det(self.S[k])
422 SFI += self.w_AB[k] * np.linalg.det(
423 np.transpose(np.conjugate(self.S[k])))
424 for k in range(len(self.w_AC)):
425 A += self.w_AC[k] * np.linalg.det(self.A[k])
426 for k in range(len(self.w_BC)):
427 B += self.w_BC[k] * np.linalg.det(self.B[k])
429 if self.energy_gap is None:
430 self.get_energy_gap()
431 dE_if = self.energy_gap
432 if np.abs(dE_if) < 0.001 and np.abs(A**2 - B**2) < 0.001:
433 warnings.warn(migliore_warning)
434 self.ct = np.abs((A * B / (A**2 - B**2)) * dE_if *
435 (1 - SIF * (A**2 + B**2) / (2 * A * B)) *
436 1 / (1 - SIF**2))
437 return self.ct
439 def make_hamiltonian_matrix(self):
440 # this returns a 2x2 cDFT Hamiltonian
441 # in a non-orthogonal basis -->
442 # made orthogonal in get_coupling_term
443 # self.get_diagonal_H_elements()
444 H_AA = self.H_AA # diabat A with H^KS_A
445 H_BB = self.H_BB # diabat B with H^KS_B
447 if hasattr(self, 'S'):
448 S = self.S
449 else:
450 S = self.get_overlap_matrix()
451 self.S = S
453 if hasattr(self, 'VW'):
454 VW = self.VW
455 else:
456 VW = self.get_weight_matrix()
457 if rank == 0:
458 S_AB = S[0][1]
459 S_BA = S[1][0]
461 h_AB = self.FB * S_AB - VW[0][1]
462 h_BA = self.FA * S_BA - VW[1][0]
464 # Ensure that H is hermitian
465 H_AB = 1. / 2. * (h_AB + h_BA)
466 H_BA = 1. / 2. * (h_AB + h_BA).conj()
468 self.H = np.array([[H_AA, H_BA], [H_AB, H_BB]])
469 return self.H
471 def get_ae_pair_weight_matrix(self):
472 if self.calc_A.wfs.kd.nibzkpts != 1:
473 raise ValueError(ae_ibz_error)
474 if hasattr(self, 'n_ab'):
475 pass
476 else:
477 raise ValueError(nab_missing_error)
479 # pseudo wfs to all-electron wfs
480 psi_A = PS2AE(self.calc_A, grid_spacing=self.h)
481 psi_B = PS2AE(self.calc_B, grid_apspacing=self.h)
483 ns = self.calc_A.wfs.nspins
485 # weight functions sum_c VcWc
486 wa = []
487 wb = []
489 for a in range(len(self.Va)):
490 wa.append(
491 interpolate_weight(self.calc_A, self.fineweightA[a], h=self.h))
492 for b in range(len(self.Vb)):
493 wb.append(
494 interpolate_weight(self.calc_B, self.fineweightB[b], h=self.h))
495 wa = np.asarray(wa)
496 wb = np.asarray(wb)
498 # check number of occupied and total number of bands
499 n_occup_A = self.get_n_occupied_bands(self.calc_A)
500 n_occup_B = self.get_n_occupied_bands(self.calc_B)
502 # place to store <i_A(k)|w|j_B(k)>
503 w_kij_AB = []
504 w_kij_BA = []
506 # k-point weights
507 w_k = np.zeros(self.calc_A.wfs.kd.nibzkpts, dtype=np.float32)
509 # get weight matrix at for each ij band at kpt and spin
510 # the resulting matrix is organized in alpha and beta blocks
511 # | | |
512 # | a | 0 | a:<psi_a|Vb*wb|psi_a> != 0
513 # VW_AB =|_______|____| , <psi_a|w|psi_b> = 0
514 # | 0 | b | b:<psi_b|Vb*wb|psi_b> != 0
515 # | | |
516 #
517 # a = nAa x nAa, b = nAb x nAb
519 for spin in range(ns):
520 for k in range(self.calc_A.wfs.kd.nibzkpts):
522 # k-dependent overlap/pair density matrices
523 inv_S = np.linalg.inv(self.n_ab[k])
524 det_S = np.linalg.det(self.n_ab[k])
525 I = np.identity(inv_S.shape[0])
526 C_ab = np.transpose(np.dot(inv_S, (det_S * I)))
528 nAa, nAb, nBa, nBb = self.check_bands(n_occup_A, n_occup_B, k)
529 nas, n_occup, n_occup_s = self.check_spin_and_occupations(
530 nAa, nAb, nBa, nBb)
532 # check that a and b cDFT states have similar spin state
533 if np.sign(nAa - nAb) != np.sign(nBa - nBb):
534 warning = UserWarning(spin_state_error)
535 warnings.warn(warning)
536 # form overlap matrices of correct size for each kpt
538 if spin == 0:
539 w_kij_AB.append(
540 np.zeros((n_occup, n_occup), dtype=complex))
541 w_kij_BA.append(
542 np.zeros((n_occup, n_occup), dtype=complex))
544 # store k-point weights
545 kd = self.calc_A.wfs.kd
546 w_kA = kd.weight_k[k]
547 kd = self.calc_B.wfs.kd
548 w_kB = kd.weight_k[k]
550 for i in range(n_occup_s[spin]):
551 for j in range(n_occup_s[spin]):
552 I = spin * nas + i
553 J = spin * nas + j
555 psi_kA = psi_A.get_wave_function(n=i,
556 k=k,
557 s=spin,
558 ae=True)
559 psi_kB = psi_B.get_wave_function(n=j,
560 k=k,
561 s=spin,
562 ae=True)
563 w_ij_AB = []
564 w_ji_BA = []
566 for b in range(len(self.Vb)):
567 integral = psi_B.gd.integrate(
568 psi_kA.conj() * wb[b] * psi_kB,
569 global_integral=True) * C_ab[I][J]
571 if b >= self.n_charge_regionsB and spin == 1:
572 # for charge constraint w > 0
573 integral *= -1.
574 w_ij_AB.append(-integral)
576 for a in range(len(self.Va)):
577 integral = psi_A.gd.integrate(
578 psi_kB.conj() * wa[a] * psi_kA,
579 global_integral=True) * C_ab[J][I]
580 if a >= self.n_charge_regionsA and spin == 1:
581 integral *= -1.
582 w_ji_BA.append(-integral)
584 w_ij_AB = np.asarray(w_ij_AB) * Bohr**3
585 w_ji_BA = np.asarray(w_ji_BA) * Bohr**3
587 # collect kpt weight, only once per kpt
588 if spin == 0 and i == 0 and j == 0:
589 w_k[k] = (w_kA + w_kB) / 2.
591 w_kij_AB[k][I][J] += np.dot(self.Vb, w_ij_AB).sum()
592 w_kij_BA[k][J][I] += np.dot(self.Va, w_ji_BA).sum()
594 self.w_k = w_k
595 self.VW_AB = w_kij_AB
596 self.VW_BA = w_kij_BA
598 if self.save_matrix:
599 np.save(self.VW_matrix + 'final_AB', self.VW_AB)
600 np.save(self.VW_matrix + 'final_BA', self.VW_BA)
601 return self.VW_AB, self.VW_BA, self.w_k
603 def get_pair_weight_matrix(self):
604 # <Psi_A|Psi_B> using pseudo wave functions and atomic corrections
606 if not hasattr(self, 'n_ab'):
607 raise ValueError(nab_missing_error)
609 # check number of occupied and total number of bands
610 n_occup_A = self.get_n_occupied_bands(
611 self.calc_A) # total of filled a and b bands
612 n_occup_B = self.get_n_occupied_bands(self.calc_B)
614 # place to store <i_A(k)|w|j_B(k)>
615 w_kij_AB = []
616 w_kij_BA = []
617 # k-point weights
618 w_k = np.zeros(self.calc_A.wfs.kd.nibzkpts)
620 # get weight matrix at for each ij band at kpt and spin
621 # the resulting matrix is organized in alpha and beta blocks
622 # | | |
623 # | a | 0 | a:<psi_a|Vb*wb|psi_a> != 0
624 # VW_AB =|_______|____| , <psi_a|w|psi_b> = 0
625 # | 0 | b | b:<psi_b|Vb*wb|psi_b> != 0
626 # | | |
627 #
628 # a = nAa x nAa, b = nAb x nAb
630 for kpt_a, kpt_b in zip(self.calc_A.wfs.kpt_u, self.calc_B.wfs.kpt_u):
631 k = kpt_a.k
632 spin = kpt_a.s
634 # k-dependent overlap/pair density matrices
635 inv_S = np.linalg.inv(self.n_ab[k])
636 det_S = np.linalg.det(self.n_ab[k])
637 I = np.identity(inv_S.shape[0])
638 C_ab = np.transpose(np.dot(inv_S, (det_S * I)))
640 inv_S = np.linalg.inv(np.transpose(self.n_ab[k]).conj())
641 det_S = np.linalg.det(np.transpose(self.n_ab[k]))
642 I = np.identity(inv_S.shape[0])
643 C_ba = np.transpose(np.dot(inv_S, (det_S * I)))
645 nAa, nAb, nBa, nBb = self.check_bands(n_occup_A, n_occup_B, k)
646 nas, n_occup, n_occup_s = self.check_spin_and_occupations(
647 nAa, nAb, nBa, nBb)
649 # check that a and b cDFT states have similar spin state
650 if np.sign(nAa - nAb) != np.sign(nBa - nBb):
651 warning = UserWarning(spin_state_error)
652 warnings.warn(warning)
653 # form overlap matrices of correct size for each kpt
655 if spin == 0:
656 w_kij_AB.append(np.zeros((n_occup, n_occup), dtype=complex))
657 w_kij_BA.append(np.zeros((n_occup, n_occup), dtype=complex))
659 # store k-point weights
660 kd = self.calc_A.wfs.kd
661 w_kA = kd.weight_k[k]
662 kd = self.calc_B.wfs.kd
663 w_kB = kd.weight_k[k]
665 for b in range(len(self.Vb)):
666 self.get_matrix_element(kpt_a.psit_nG,
667 kpt_a.P_ani,
668 kpt_b.psit_nG,
669 kpt_b.P_ani,
670 n_occup_s,
671 n_occup_s,
672 w_kij_AB,
673 self.regionsB[b],
674 k,
675 spin,
676 nas,
677 V=self.Vb[b],
678 W=self.coarseweightB[b],
679 C12=C_ab)
681 if b >= self.n_charge_regionsB and spin == 1:
682 # change sign of b spin constraint
683 w_kij_AB[k][nas:, nas:] *= -1.
685 for a in range(len(self.Va)):
686 self.get_matrix_element(kpt_b.psit_nG,
687 kpt_b.P_ani,
688 kpt_a.psit_nG,
689 kpt_a.P_ani,
690 n_occup_s,
691 n_occup_s,
692 w_kij_BA,
693 self.regionsA[a],
694 k,
695 spin,
696 nas,
697 V=self.Va[a],
698 W=self.coarseweightA[a],
699 C12=C_ba)
700 if a >= self.n_charge_regionsA and spin == 1:
701 w_kij_BA[k][nas:, nas:] *= -1.
703 if spin == 0:
704 w_k[k] = (w_kA + w_kB) / 2.
706 self.VW_BA = np.asarray(w_kij_BA)
707 self.VW_AB = np.asarray(w_kij_AB)
708 self.w_k = w_k
710 return self.VW_AB, self.VW_BA, self.w_k
712 def get_weight_matrix(self):
713 # from the pair density matrix
714 if not (hasattr(self, 'VW_AB')):
715 if self.AE:
716 self.get_ae_pair_weight_matrix()
717 else:
718 self.get_pair_weight_matrix()
719 if not (hasattr(self, 'w_k')):
720 self.get_ae_pair_density_matrix()
722 if rank == 0:
723 # diagonal of V*weight matrix
724 self.VW = np.zeros((2, 2))
725 # fill diagonal
727 self.VW[0][0] += np.sum(self.NA)
728 self.VW[1][1] += np.sum(self.NB)
730 W_k_AB = np.zeros(len(self.w_k))
731 W_k_BA = np.zeros(len(self.w_k))
733 for k in range(len(self.w_k)):
734 # sum_k (sum_ij <i|sum_c Vc*wc| j> * C_ij
735 W_k_AB[k] = self.VW_AB[k].real.sum()
736 W_k_AB[k] *= self.w_k[k]
737 W_k_BA[k] = self.VW_BA[k].real.sum()
738 W_k_BA[k] *= self.w_k[k]
740 self.VW[0][1] = W_k_AB.sum()
741 self.VW[1][0] = W_k_BA.sum()
743 return self.VW
745 def get_ae_pair_density_matrix(self, calc_A, calc_B, matrix_name=None):
746 if calc_A.wfs.kd.nibzkpts != 1:
747 raise ValueError(ae_ibz_error)
748 # <Psi_A|Psi_B> using the all-electron pair density
749 psi_A = PS2AE(calc_A, grid_spacing=self.h)
750 psi_B = PS2AE(calc_B, grid_spacing=self.h)
752 ns = calc_A.wfs.nspins
754 # total of filled a and b bands for each spin and kpt
755 n_occup_A = self.get_n_occupied_bands(calc_A)
756 n_occup_B = self.get_n_occupied_bands(calc_B)
758 # list to store k-dependent pair density
759 n_AB = []
760 w_k = np.zeros(calc_A.wfs.kd.nibzkpts) # store kpt weights
762 # get overlap at for each ij band at kpt and spin
763 # the resulting matrix is organized in alpha and beta blocks
764 # | | |
765 # | a | 0 | a:<psi_a|psi_a> != 0
766 # S =|_______|____| , <psi_a|psi_b> = 0
767 # | 0 | b | b:<psi_b|psi_b> != 0
768 # | | |
769 #
770 # a = naa x naa, b = nab x nab
772 for spin in range(ns):
773 for k in range(calc_A.wfs.kd.nibzkpts):
775 nAa, nAb, nBa, nBb = self.check_bands(n_occup_A, n_occup_B, k)
776 nas, n_occup, n_occup_s = self.check_spin_and_occupations(
777 nAa, nAb, nBa, nBb)
778 kd = calc_A.wfs.kd
779 w_kA = kd.weight_k[k]
780 kd = calc_B.wfs.kd
781 w_kB = kd.weight_k[k]
783 # check that a and b cDFT states have similar spin state
784 if np.sign(nAa - nAb) != np.sign(nBa - nBb):
785 warning = UserWarning(
786 'The cDFT wave functions have\n'
787 'different spin states! Similar\n'
788 'spin states are required for coupling constant\n'
789 'calculation!')
790 warnings.warn(warning)
791 # form overlap matrices of correct size for each kpt
792 if spin == 0:
793 n_AB.append(np.zeros((n_occup, n_occup), dtype=complex))
795 for i in range(n_occup_s[spin]):
796 for j in range(n_occup_s[spin]):
797 # take only the bands which contain electrons
798 # in spin-orbital
799 psi_kA = psi_A.get_wave_function(n=i,
800 k=k,
801 s=spin,
802 ae=True)
803 psi_kB = psi_B.get_wave_function(n=j,
804 k=k,
805 s=spin,
806 ae=True)
807 n_ij = psi_A.gd.integrate(psi_kA.conj(),
808 psi_kB,
809 global_integral=True)
811 if spin == 0 and i == 0 and j == 0:
812 w_k[k] = (w_kA + w_kB) / 2.
814 I = spin * nas + i
815 J = spin * nas + j
816 n_AB[k][I][J] = n_ij * Bohr**3
818 n_AB = np.asarray(n_AB)
819 self.w_k = w_k
820 self.n_ab = n_AB
822 if self.save_matrix:
823 if matrix_name is None:
824 np.save(self.S_matrix + 'final', self.n_ab)
825 else:
826 np.save('%s_final' % matrix_name, self.n_ab)
827 return self.n_ab, self.w_k
829 def get_pair_density_matrix(self, calc_A, calc_B, matrix_name=None):
830 # <Psi_A|Psi_B> using pseudo wave functions and atomic corrections
831 assert calc_A.wfs.kd.nibzkpts == calc_B.wfs.kd.nibzkpts
833 # total of filled a and b bands for each spin and kpt
834 n_occup_A = self.get_n_occupied_bands(calc_A)
835 n_occup_B = self.get_n_occupied_bands(calc_B)
837 # list to store k-dependent pair density
838 n_AB = []
840 w_k = np.zeros(calc_A.wfs.kd.nibzkpts) # store kpt weights
842 # get overlap at for each ij band at kpt and spin
843 # the resulting matrix is organized in alpha and beta blocks
844 # | | |
845 # | a | 0 | a:<psi_a|psi_a> != 0
846 # S =|_______|____| , <psi_a|psi_b> = 0
847 # | 0 | b | b:<psi_b|psi_b> != 0
848 # | | |
849 #
850 # a = naa x naa, b = nab x nab
852 for kpt_a, kpt_b in zip(calc_A.wfs.kpt_u, calc_B.wfs.kpt_u):
854 k = kpt_a.k
855 spin = kpt_a.s
856 nAa, nAb, nBa, nBb = self.check_bands(n_occup_A, n_occup_B, k)
857 nas, n_occup, n_occup_s = self.check_spin_and_occupations(
858 nAa, nAb, nBa, nBb)
860 # check that a and b cDFT states have similar spin state
861 if np.sign(nAa - nAb) != np.sign(nBa - nBb):
862 warning = UserWarning(
863 'The cDFT wave functions have\n'
864 'different spin states! Similar\n'
865 'spin states are required for coupling constant\n'
866 'calculation!')
867 warnings.warn(warning)
868 # form overlap matrices of correct size for each kpt
869 if spin == 0:
870 n_AB.append(np.zeros((n_occup, n_occup), dtype=complex))
872 kd = calc_A.wfs.kd
873 w_kA = kd.weight_k[k]
874 kd = calc_B.wfs.kd
875 w_kB = kd.weight_k[k]
877 self.get_matrix_element(kpt_b.psit_nG, kpt_b.P_ani, kpt_a.psit_nG,
878 kpt_a.P_ani, n_occup_s, n_occup_s, n_AB,
879 None, k, spin, nas)
881 if spin == 0:
882 w_k[k] = (w_kA + w_kB) / 2.
883 n_AB = np.asarray(n_AB)
884 self.w_k = w_k
885 self.n_ab = n_AB
887 if self.save_matrix:
888 if matrix_name is None:
889 np.save(self.S_matrix + 'final', self.n_ab)
890 else:
891 np.save('%s_final' % matrix_name, self.n_ab)
893 return self.n_ab, self.w_k
895 def get_overlap_matrix(self, save=False, name='wf_overlap'):
897 # from the pair density matrix
898 if not (hasattr(self, 'n_ab')):
899 if self.AE:
900 self.n_ab, self.w_k = self.get_ae_pair_density_matrix(
901 self.calc_A, self.calc_B)
902 else:
903 self.n_ab, self.w_k = self.get_pair_density_matrix(
904 self.calc_A, self.calc_B)
906 if not (hasattr(self, 'w_k')):
907 self.get_ae_pair_density_matrix(self.calc_A, self.calc_B)
908 if rank == 0:
909 self.S = np.identity(2)
911 S_k_AB = np.zeros(len(self.w_k))
912 S_k_BA = np.zeros(len(self.w_k))
914 for k in range(len(self.w_k)):
915 S_k_AB[k] = self.w_k[k] * np.linalg.det(self.n_ab[k]).real
916 # determinant of the complex conjugate
917 S_k_BA[k] = self.w_k[k] * np.linalg.det(
918 np.transpose(self.n_ab[k]).conj()).real
920 S_AB = S_k_AB.sum()
921 S_BA = S_k_BA.sum()
923 # fill 2x2 overlap matrix
924 self.S[0][1] = S_AB
925 self.S[1][0] = S_BA
927 if save:
928 np.save('%s' % name, self.S)
929 return self.S
931 def get_migliore_wf_overlaps(self, calc_A, calc_B, calc_gs):
932 '''A = <I|GS>, B= <F|GS>, S = <I|F> '''
933 if self.AE:
934 S, w_AB = self.get_ae_pair_density_matrix(calc_A,
935 calc_B,
936 matrix_name='Sif')
937 A, w_AC = self.get_ae_pair_density_matrix(calc_A,
938 calc_gs,
939 matrix_name='A')
940 B, w_BC = self.get_ae_pair_density_matrix(calc_B,
941 calc_gs,
942 matrix_name='B')
943 else:
944 S, w_AB = self.get_pair_density_matrix(calc_A,
945 calc_B,
946 matrix_name='Sif')
947 A, w_AC = self.get_pair_density_matrix(calc_A,
948 calc_gs,
949 matrix_name='A')
950 B, w_BC = self.get_pair_density_matrix(calc_B,
951 calc_gs,
952 matrix_name='B')
954 return S, A, B, w_AB, w_AC, w_BC
956 def get_matrix_element(self,
957 psit1_nG,
958 P1_ani,
959 psit2_nG,
960 P2_ani,
961 n_occup_1,
962 n_occup_2,
963 vw,
964 region,
965 k,
966 s,
967 nas,
968 V=1.,
969 W=1.,
970 C12=None):
971 # weight: V*W acting on psitb_nG: VW_ij = <psita_nG_i|VW|psitb_nG_j> or
972 # overlap <psita_nG_i|psitb_nG_j>
973 # includes occupation dependency and only states with filled
974 # orbitals are considered
975 # Cab = cofactor matrix
976 # k and s --> kpt and spin
978 VW_ij = np.zeros((len(psit1_nG), len(psit2_nG)), complex)
979 VW = V * W
981 for n1 in range(len(psit1_nG)):
982 for n2 in range(len(psit2_nG)):
983 nijt_G = np.multiply(psit1_nG[n1].conj(), psit2_nG[n2])
984 VW_ij[n1][n2] = self.gd.integrate(VW * nijt_G,
985 global_integral=True)
987 P_array = np.zeros(VW_ij.shape, complex)
988 for a, P1_ni in P1_ani.items():
989 P2_ni = P2_ani[a]
990 if region is None or a in region:
991 # the atomic correction is zero outside
992 # the augmentation regions --> w=0 if
993 # a is not in the w.
994 # region = None --> overlap and all terms are used
995 inner = np.dot(P1_ni, self.calc_A.wfs.setups[a].dO_ii)
996 outer = (np.dot(P2_ni, V * np.conjugate(inner.T))).T
997 P_array += outer
999 self.calc_A.density.gd.comm.sum(P_array)
1000 VW_ij += P_array
1002 for i, j in enumerate(range(n_occup_1[s])):
1003 for x, y in enumerate(range(n_occup_2[s])):
1004 I = s * nas + i
1005 X = s * nas + x
1006 if C12 is None:
1007 vw[k][I][X] += VW_ij[j][y]
1008 else:
1009 vw[k][I][X] += VW_ij[j][y] * C12[I][X]
1011 def get_reorganization_energy(self):
1012 # get Ea (Rb) - Ea(Ra) -->
1013 # cdft energy at geometry Rb
1014 # with charge constraint A
1015 geometry = self.calc_B.get_atoms()
1017 cdft = self.cdft_A
1018 # set cdft_a on geometry of B
1019 geometry.calc = cdft
1020 self.reorg = geometry.get_potential_energy()
1021 self.reorg -= self.EA
1022 return self.reorg
1024 def get_energy_gap(self):
1025 # get Eb (Ra) - Ea(Ra) -->
1026 # cdft energy at geometry Ra
1027 # with charge constraint A or B
1028 geometry = self.calc_A.get_atoms()
1030 cdft = self.cdft_B
1031 # set cdft_a on geometry of B
1032 geometry.calc = cdft
1033 self.energy_gap = geometry.get_potential_energy()
1034 self.energy_gap -= self.EA
1035 return self.energy_gap
1037 def get_landau_zener(self):
1038 # computes the Landau-Zener factor
1040 planck = 4.135667e-15 # eV s
1042 if self.reorg is not None:
1043 Lambda = self.reorg
1044 else:
1045 Lambda = self.get_reorganization_energy()
1047 if hasattr(self, 'ct'):
1048 Hab = self.ct
1049 else:
1050 Hab = self.get_coupling_term()
1052 self.P_lz = 1 - np.exp(-np.pi**(3 / 2) * (np.abs(Hab))**2 /
1053 (planck * self.freq *
1054 np.sqrt(self.temp * Lambda * kb)))
1056 return self.P_lz
1058 def get_marcus_rate(self):
1059 # adiabatic transition state energy
1060 if self.E_tst:
1061 E_tst = self.E_tst
1062 else:
1063 # compute the activation energy
1064 # from the classical marcus
1065 # parabolas
1066 E_tst = self.get_marcus_barrier()
1068 # electron transmission coeffiecient
1069 # is the reaction diabatic or adiabatic?
1070 P_lz = self.get_landau_zener()
1071 dE = self.EA - self.EB
1072 if dE >= -self.reorg:
1073 # normal
1074 kappa = 2. * P_lz / (1 + P_lz)
1075 else:
1076 kappa = 2 * P_lz * (1 - P_lz)
1077 rate = kappa * self.freq * np.exp(-E_tst / (kb * self.temp))
1079 return rate
1081 def get_marcus_barrier(self):
1082 # approximate barrier from
1083 # two marcus parabolas
1084 # and an adiabatic correction
1086 # reaction energy
1087 dE = self.reaction_energy
1089 # crossing of the parabolas
1090 barrier = 1 / (4 * self.reorg) * (self.reorg + dE)**2
1091 if self.use_adiabatic_correction:
1092 # adiabatic correction
1093 correction = (
1094 np.abs(self.ct) +
1095 (self.reorg + dE) / 2. -
1096 np.sqrt(((self.reorg + self.dE)**2) / 4. +
1097 (np.abs(self.ct))**2))
1099 return barrier - correction
1100 else:
1101 return barrier
1103 def get_n_occupied_bands(self, calc):
1104 ''' how many occupied bands?'''
1106 ns = calc.wfs.nspins
1107 occup_ks = np.zeros((len(calc.wfs.kd.weight_k), ns), dtype=int)
1108 w = calc.wfs.kd.weight_k
1109 for k in range(calc.wfs.kd.nibzkpts):
1110 for s in range(2):
1111 f_n = calc.get_occupation_numbers(kpt=k, spin=s)
1112 f_n *= 1. / w[k]
1113 f_N = f_n > self.band_occupation_cutoff
1114 occup_ks[k][s] += f_N.sum()
1115 return occup_ks
1117 def check_bands(self, n_occup_A, n_occup_B, k):
1118 if not self.specific_bands_A and not self.specific_bands_B:
1119 nAa, nAb = n_occup_A[k][0], n_occup_A[k][1]
1120 nBa, nBb = n_occup_B[k][0], n_occup_B[k][1]
1121 else: # choose subset of orbitals for coupling
1122 if self.specific_bands_A:
1123 nAa, nAb = self.specific_bands_A[0], self.specific_bands_A[1]
1124 if self.specific_bands_B:
1125 nBa, nBb = self.specific_bands_B[0], self.specific_bands_B[1]
1127 return nAa, nAb, nBa, nBb
1129 def check_spin_and_occupations(self, nAa, nAb, nBa, nBb):
1130 nas = np.max((nAa, nBa))
1131 nbs = np.max((nAb, nBb))
1132 n_occup_s = [nas, nbs]
1133 n_occup = sum(n_occup_s)
1135 return nas, n_occup, n_occup_s