Coverage for gpaw/lrtddft2/__init__.py: 71%
181 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1"""Module for linear response TDDFT class with indexed K-matrix storage."""
3import os
4import datetime
5import glob
7import numpy as np
9from ase.units import Hartree
10from ase.utils import IOContext
12from gpaw.xc import XC
14# a KS determinant with a single occ-uncc excitation
15# from gpaw.lrtddft2.ks_singles import KohnShamSingleExcitation
17# a list of KS determinants with single occ-uncc excitations
18from gpaw.lrtddft2.ks_singles import KohnShamSingles
20# Matrix Kip,jq <ia|f_Hxc|jq>
21from gpaw.lrtddft2.k_matrix import Kmatrix
23# a set of linear combinations of KS single determinants
24from gpaw.lrtddft2.lr_transitions import LrtddftTransitions
26# a linear combination of KS single determinants
27# for a CW laser with Lorentzian width (in energy)
28from gpaw.lrtddft2.lr_response import LrResponse
30# communicators
31from gpaw.lrtddft2.lr_communicators import LrCommunicators
34class LrTDDFT2:
35 """Linear response TDDFT (Casida) class with indexed K-matrix storage."""
37 def __init__(self,
38 basefilename,
39 gs_calc,
40 fxc=None,
41 min_occ=None,
42 max_occ=None,
43 min_unocc=None,
44 max_unocc=None,
45 max_energy_diff=None,
46 recalculate=None,
47 lr_communicators=None,
48 txt='-'):
49 """Initialize linear response TDDFT without calculating anything.
51 Note
52 ----
53 Does NOT support spin polarized calculations yet.
55 Tip
56 ---
57 If K_matrix file is too large and you keep running out of
58 memory when trying to calculate spectrum or response wavefunction,
59 you can try
60 ``split -l 100000
61 xxx.K_matrix.ddddddofDDDDDD xxx.K_matrix.ddddddofDDDDDD``.
64 Parameters
65 ----------
66 basefilename
67 All files associated with this calculation are stored as
68 *<basefilename>.<extension>*
69 gs_calc
70 Ground state calculator (if you are using eh_communicator,
71 you need to take care that calc has suitable dd_communicator.)
72 fxc
73 Name of the exchange-correlation kernel (fxc) used in calculation.
74 (optional)
75 min_occ
76 Index of the first occupied state to be included in the calculation.
77 (optional)
78 max_occ
79 Index of the last occupied state (inclusive) to be included in the
80 calculation. (optional)
81 min_unocc
82 Index of the first unoccupied state to be included in the
83 calculation. (optional)
84 max_unocc
85 Index of the last unoccupied state (inclusive) to be included in the
86 calculation. (optional)
87 max_energy_diff
88 Noninteracting Kohn-Sham excitations above this value are not
89 included in the calculation. Units: eV (optional)
90 recalculate
91 | Force recalculation.
92 | 'eigen' : recalculate only eigensystem (useful for on-the-fly
93 | spectrum calculations and convergence checking)
94 | 'matrix' : recalculate matrix without solving the eigensystem
95 | 'all' : recalculate everything
96 | None : do not recalculate anything if not needed (default)
97 lr_communicators
98 Communicators for parallelizing over electron-hole pairs (i.e.,
99 rows of K-matrix) and domain. Note that ground state calculator
100 must have a matching (domain decomposition) communicator, which
101 can be assured by using lr_communicators
102 to create both communicators.
103 txt
104 Filename for text output
105 """
107 # Save input params
108 self.basefilename = basefilename
109 self.fxc_name = fxc
110 self.xc = XC(self.fxc_name)
111 self.min_occ = min_occ
112 self.max_occ = max_occ
113 self.min_unocc = min_unocc
114 self.max_unocc = max_unocc
115 if max_energy_diff is not None:
116 self.max_energy_diff = max_energy_diff / Hartree
117 else:
118 self.max_energy_diff = None
119 self.recalculate = recalculate
120 # Don't init calculator yet if it's not needed (to save memory)
121 self.calc = gs_calc
122 self.calc_ready = False
124 # FIXME: SUPPORT ALSO SPIN POLARIZED
125 self.kpt_ind = 0
127 # Input paramers?
128 self.deriv_scale = 1e-5 # fxc finite difference step
129 # ignore transition if population difference is below this value:
130 self.min_pop_diff = 1e-3
132 # set up communicators
133 self.lr_comms = lr_communicators
135 if self.lr_comms is None:
136 self.lr_comms = LrCommunicators(None, None)
137 self.lr_comms.initialize(gs_calc)
139 # Init text output
140 self.iocontext = IOContext()
141 self.txt = self.iocontext.openfile(txt, self.lr_comms.parent_comm)
143 # Check and set unset params
144 kpt = self.calc.wfs.kpt_u[self.kpt_ind]
145 nbands = len(kpt.f_n)
147 # If min/max_occ/unocc were not given, but max_energy_diff was,
148 # check that calc has enough states for max_energy_diff
149 # (i.e., KS eigenvalue difference between HOMO and highest
150 # state is below max_energy_diff)
151 if ((self.min_occ is None or self.min_unocc is None
152 or self.max_occ is None or self.max_unocc is None)
153 and self.max_energy_diff is not None):
154 n_homo = np.sum(kpt.f_n > self.min_pop_diff) - 1
155 n_highest = nbands - 1 # XXX use highest converged state instead
156 eps_n = kpt.eps_n
157 eps_diff = eps_n[n_highest] - eps_n[n_homo]
158 if eps_diff <= self.max_energy_diff:
159 msg = ('Error in LrTDDFT2: not enough states in '
160 'the calculator for the requested max_energy_diff='
161 f'{self.max_energy_diff * Hartree:.4f} eV. '
162 f'Max eigenvalue difference from HOMO (n={n_homo}) is '
163 f'{eps_diff * Hartree:.4f} eV.')
164 raise RuntimeError(msg)
166 # If min/max_occ/unocc were not given, initialized them to include
167 # everything: min_occ/unocc => 0, max_occ/unocc to nubmer of wfs,
168 # energy diff to numerical infinity
169 if self.min_occ is None:
170 self.min_occ = 0
171 if self.min_unocc is None:
172 self.min_unocc = self.min_occ
173 if self.max_occ is None:
174 self.max_occ = nbands - 1
175 if self.max_unocc is None:
176 self.max_unocc = self.max_occ
177 if self.max_energy_diff is None:
178 self.max_energy_diff = np.inf
180 self.min_occ = max(self.min_occ, 0)
181 self.min_unocc = max(self.min_unocc, 0)
183 if self.max_occ >= nbands:
184 raise RuntimeError('Error in LrTDDFT2: max_occ >= nbands')
185 if self.max_unocc >= nbands:
186 raise RuntimeError('Error in LrTDDFT2: max_unocc >= nbands')
188 # Only spin unpolarized calculations are supported atm
189 # > FIXME
190 assert len(self.calc.wfs.kpt_u) == 1, \
191 'LrTDDFT2 does not support more than one k-point/spin.'
192 self.kpt_ind = 0
193 # <
195 # Internal classes
197 # list of singly excited Kohn-Sham Slater determinants
198 # (ascending KS energy difference)
199 self.ks_singles = KohnShamSingles(self.basefilename, self.calc,
200 self.kpt_ind, self.min_occ,
201 self.max_occ, self.min_unocc,
202 self.max_unocc, self.max_energy_diff,
203 self.min_pop_diff, self.lr_comms,
204 self.txt)
206 # Response kernel matrix K = <ip|f_Hxc|jq>
207 # Note: this is not the Casida matrix
208 self.K_matrix = Kmatrix(self.ks_singles, self.xc, self.deriv_scale)
210 self.sl_lrtddft = self.calc.parallel['sl_lrtddft']
212 # LR-TDDFT transitions
213 self.lr_transitions = LrtddftTransitions(self.ks_singles,
214 self.K_matrix,
215 self.sl_lrtddft)
217 # Response wavefunction
218 self.lr_response_wf = None
220 # Pair density
221 self.pair_density = None # pair density class
223 # Timer
224 self.timer = self.calc.timer
225 self.timer.start('LrTDDFT')
227 # If a previous calculation exist, read info file
228 # self.read(self.basefilename)
230 # Write info file
231 # self.parent_comm.barrier()
232 # if self.parent_comm.rank == 0:
233 # self.write_info(self.basefilename)
235 # self.custom_axes = None
236 # self.K_matrix_scaling_factor = 1.0
237 self.K_matrix_values_ready = False
239 def get_transitions(self,
240 filename=None,
241 min_energy=0.0,
242 max_energy=30.0,
243 units='eVcgs'):
244 """Get transitions: energy, dipole strength and rotatory strength.
246 Returns transitions as (w,S,R, Sx,Sy,Sz) where
247 w is an array of frequencies,
248 S is an array of corresponding dipole strengths,
249 and R is an array of corresponding rotatory strengths.
251 Parameters
252 ----------
253 min_energy
254 Minimum energy
255 min_energy
256 Maximum energy
257 units
258 Units for spectrum: 'au' or 'eVcgs'
259 """
261 self.calculate()
262 self.txt.write('Calculating transitions (%s).\n' %
263 str(datetime.datetime.now()))
264 trans = self.lr_transitions.get_transitions(filename, min_energy,
265 max_energy, units)
266 self.txt.write('Transitions calculated (%s).\n' %
267 str(datetime.datetime.now()))
268 return trans
270 #################################################################
271 def get_spectrum(self,
272 filename=None,
273 min_energy=0.0,
274 max_energy=30.0,
275 energy_step=0.01,
276 width=0.1,
277 units='eVcgs'):
278 """Get spectrum for dipole and rotatory strength.
280 Returns folded spectrum as (w,S,R) where w is an array of frequencies,
281 S is an array of corresponding dipole strengths, and R is an array of
282 corresponding rotatory strengths.
284 Parameters
285 ----------
286 min_energy
287 Minimum energy
288 min_energy
289 Maximum energy
290 energy_step
291 Spacing between calculated energies
292 width
293 Width of the Gaussian
294 units
295 Units for spectrum: 'au' or 'eVcgs'
296 """
297 self.calculate()
298 self.txt.write('Calculating spectrum (%s).\n' %
299 str(datetime.datetime.now()))
300 spec = self.lr_transitions.get_spectrum(filename, min_energy,
301 max_energy, energy_step, width,
302 units)
303 self.txt.write('Spectrum calculated (%s).\n' %
304 str(datetime.datetime.now()))
305 return spec
307 #################################################################
308 def get_transition_contributions(self, index_of_transition):
309 """Get contributions of Kohn-Sham singles to a given transition
310 as number of electrons contributing.
312 Includes population difference.
314 This method is meant to be used for small number of transitions.
315 It is not suitable for analysing densely packed transitions of
316 large systems. Use transition contribution map (TCM) or similar
317 approach for this.
319 Parameters
320 ----------
321 index_of_transition:
322 index of transition starting from zero
323 """
324 self.calculate()
325 return self.lr_transitions.get_transition_contributions(
326 index_of_transition)
328 def calculate_response(self,
329 excitation_energy,
330 excitation_direction,
331 lorentzian_width,
332 units='eVang'):
333 """Calculates and returns response using TD-DFPT.
335 Parameters
336 ----------
337 excitation_energy
338 Energy of the laser in given units
339 excitation_direction
340 Vector for direction (will be normalized)
341 lorentzian_width
342 Life time or width parameter. Larger width results in wider
343 energy envelope around excitation energy.
344 """
346 # S_z(omega) = 2 * omega sum_ip n_ip C^(im)_ip(omega) * mu^(z)_ip
348 self.calculate()
350 omega_au = excitation_energy
351 width_au = lorentzian_width
352 # always unit field in au !!!
353 direction_au = np.array(excitation_direction)
354 direction_au = direction_au / np.sqrt(
355 np.vdot(direction_au, direction_au))
357 if units == 'au':
358 pass
359 elif units == 'eVang':
360 omega_au /= Hartree
361 width_au /= Hartree
362 else:
363 raise RuntimeError(
364 'Error in calculate_response_wavefunction: Invalid units.')
366 lr_response = LrResponse(self, omega_au, direction_au, width_au,
367 self.sl_lrtddft)
368 lr_response.solve()
370 return lr_response
372 def calculate(self):
373 """Calculates linear response matrix and properties of KS
374 electron-hole pairs.
376 This is called implicitly by get_spectrum, get_transitions, etc.
377 but there is no harm for calling this explicitly.
378 """
379 if not self.calc_ready:
380 # Initialize wfs, paw corrections and xc, if not done yet
381 # FIXME: CHECK THIS STUFF, DOES GLLB WORK???
382 if not self.calc_ready:
383 mode = self.calc.wfs.mode
384 C_nM = self.calc.wfs.kpt_u[self.kpt_ind].C_nM
385 psit_nG = self.calc.wfs.kpt_u[self.kpt_ind].psit_nG
386 if ((mode == 'lcao' and C_nM is None)
387 or (mode == 'fd' and psit_nG is None)):
388 raise RuntimeError('Use ground state calculator ' +
389 'containing the wave functions.')
391 if not self.calc.scf.converged:
392 # Do not allow new scf cycles, it could change
393 # the wave function signs after every restart.
394 raise RuntimeError('Converge wave functions first.')
396 spos_ac = self.calc.initialize_positions()
397 self.calc.wfs.calculate_occupation_numbers()
398 self.calc.wfs.initialize(self.calc.density,
399 self.calc.hamiltonian, spos_ac)
400 self.xc.initialize(self.calc.density, self.calc.hamiltonian,
401 self.calc.wfs)
402 if mode == 'lcao':
403 self.calc.wfs.initialize_wave_functions_from_lcao()
404 self.calc_ready = True
406 # Singles logic
407 if self.recalculate == 'all' or self.recalculate == 'matrix':
408 self.ks_singles.update_list()
409 self.ks_singles.calculate()
410 elif self.recalculate == 'eigen':
411 self.ks_singles.read()
412 self.ks_singles.kss_list_ready = True
413 self.ks_singles.kss_prop_ready = True
414 elif ((not self.ks_singles.kss_list_ready)
415 or (not self.ks_singles.kss_prop_ready)):
416 self.ks_singles.read()
417 self.ks_singles.update_list()
418 self.ks_singles.calculate()
419 else:
420 pass
422 # K-matrix logic
423 if self.recalculate == 'all' or self.recalculate == 'matrix':
424 # delete files
425 if self.lr_comms.parent_comm.rank == 0:
426 for ready_file in glob.glob(self.basefilename +
427 '.ready_rows.*'):
428 os.remove(ready_file)
429 for K_file in glob.glob(self.basefilename + '.K_matrix.*'):
430 os.remove(K_file)
431 self.K_matrix.calculate()
432 elif self.recalculate == 'eigen':
433 self.K_matrix.read_indices()
434 self.K_matrix.K_matrix_ready = True
435 elif not self.K_matrix.K_matrix_ready:
436 self.K_matrix.read_indices()
437 self.K_matrix.calculate()
438 else:
439 pass
441 # Wait... we don't want to read incomplete files
442 self.lr_comms.parent_comm.barrier()
444 if not self.K_matrix_values_ready:
445 self.K_matrix.read_values()
447 # lr_transitions logic
448 if not self.lr_transitions.trans_prop_ready:
449 trans_file = self.basefilename + '.transitions'
450 if os.path.exists(trans_file) and os.path.isfile(trans_file):
451 os.remove(trans_file)
452 self.lr_transitions.calculate()
454 # recalculate only once
455 self.recalculate = None
457 def read(self, basename):
458 """Does not do much at the moment."""
459 info_file = open(basename + '.lr_info')
460 for line in info_file:
461 if line[0] == '#':
462 continue
463 if len(line.split()) <= 2:
464 continue
465 # key = line.split('=')[0]
466 # value = line.split('=')[1]
467 # .....
468 # FIXME: do something, like warn if changed
469 # ...
470 info_file.close()
472 def write_info(self, basename):
473 """Writes used parameters to a file."""
474 f = open(basename + '.lr_info', 'a+')
475 f.write('# LrTDDFTindexed\n')
476 f.write('%20s = %s\n' % ('xc_name', self.xc_name))
477 f.write('%20s = %d\n' % ('min_occ', self.min_occ))
478 f.write('%20s = %d\n' % ('min_unocc', self.min_unocc))
479 f.write('%20s = %d\n' % ('max_occ', self.max_occ))
480 f.write('%20s = %d\n' % ('max_unocc', self.max_unocc))
481 f.write('%20s = %18.12lf\n' %
482 ('max_energy_diff', self.max_energy_diff))
483 f.write('%20s = %18.12lf\n' % ('deriv_scale', self.deriv_scale))
484 f.write('%20s = %18.12lf\n' % ('min_pop_diff', self.min_pop_diff))
485 f.close()
487 def __del__(self):
488 try:
489 self.iocontext.close()
490 self.timer.stop('LrTDDFT')
491 except AttributeError:
492 pass