Coverage for strongcoca/response/time_domain_casida.py: 100%

1import logging 

2from typing import Optional, Tuple 

3import numpy as np 

4from . import CasidaResponse, TimeDomainResponse 

5from ..units import fs_to_au 

6 

7 

8logger = logging.getLogger(__name__) 

9 

10 

11def _check_propagation_parameters(eigenvals: np.ndarray, 

12 dt: Optional[float] = None, 

13 n_steps: Optional[int] = None) -> Tuple[float, int]: 

14 """Checks propagation parameters based on eigenvalues. 

15 

16 Parameters 

17 ---------- 

18 dt 

19 Time step in atomic units. 

20 

21 The time step is used to propagate the Casida system and 

22 generate the dipole tensor in the time domain; by default the time 

23 step is set based on the highest frequency of the eigenvalue spectrum. 

24 n_steps 

25 Number of time steps to propagate the Casida system and generate the 

26 dipole tensor in the time domain; by default the value is set 

27 based on the lowest frequency of the eigenvalue spectrum. 

28 

29 Returns 

30 ------- 

31 dt 

32 Time step in atomic units (the same as given or recommended). 

33 n_steps 

34 Number of time steps (the same as given or recommended). 

35 """ 

36 

37 # determine suitable integration parameters and 

38 # check sensibility of user-provided values (if applicable) 

39 # - time step 

40 dt_recommended = 1 / sorted(eigenvals)[-1] 

41 logger.debug(f'dt_recommended: {dt_recommended} a.u.') 

42 if dt is None: 

43 dt = dt_recommended 

44 else: 

45 if dt > 2 * dt_recommended or dt < 0.5 * dt_recommended: 

46 logger.warning(f'dt ({dt} a.u.) deviates strongly from the' 

47 f' recommended value ({dt_recommended} a.u.).') 

48 logger.debug(f'dt_actual: {dt} a.u.') 

49 # - maximum time and number of time steps 

50 t_max_recommended = 100 / sorted(eigenvals)[0] 

51 logger.debug(f't_max_recommended: {t_max_recommended} a.u.') 

52 if n_steps is None: 

53 n_steps = int(t_max_recommended / dt) 

54 else: 

55 t_max = dt * (n_steps - 1) 

56 if t_max < 0.5 * t_max_recommended or t_max > 2 * t_max_recommended: 

57 logger.warning(f't_max ({t_max} a.u.) deviates strongly from the' 

58 f' recommended value ({t_max_recommended} a.u.).') 

59 logger.debug(f'n_steps: {n_steps}') 

60 

61 return dt, n_steps 

62 

63 

64def build_time_domain_response_from_casida( 

65 casida_response: CasidaResponse, 

66 dt: Optional[float] = None, 

67 n_steps: Optional[int] = None, 

68 units: str = 'fs', 

69 name: str = 'PropagatedCasida') -> TimeDomainResponse: 

70 """Builds a :class:`~strongcoca.response.TimeDomainResponse` instance 

71 from an existing Casida representation via time propagation. 

72 

73 Parameters 

74 ---------- 

75 casida_response 

76 Response in Casida representation. 

77 dt 

78 Time step; units are fs by default, optionally atomic units (see :attr:`units`). 

79 n_steps 

80 Number of time steps for which to propagate the equations of motion. 

81 units 

82 `fs` to specify dt in fs or `au` to specify `dt` in atomic units. 

83 """ 

84 

85 if units == 'fs': 

86 if dt is not None: 

87 dt = dt * fs_to_au 

88 elif units != 'au': 

89 raise ValueError(f"units has to be 'fs' or 'au', not '{units}'") 

90 

91 # TODO: The argument below should probably be just `energies`. Test this. See #41. 

92 dt, n_steps = _check_propagation_parameters( 

93 casida_response.excitations._energies**2, dt, n_steps) 

94 

95 # Matrices 

96 D_n = casida_response._D_n 

97 K_nn = casida_response._K_nn 

98 mu_nv = casida_response._mu_nv 

99 

100 # Propagate 

101 time_t = np.linspace(0, dt * (n_steps - 1), n_steps) 

102 dm_tvv = np.zeros((n_steps, 3, 3)) 

103 

104 # Integrate equations of motion 

105 Re_rho_nv = 0 * mu_nv 

106 Im_rho_nv = -mu_nv 

107 for t in range(0, n_steps): 

108 # Evaluate dipole moment 

109 dm_tvv[t] = 2 * np.dot(Re_rho_nv.T, mu_nv) 

110 

111 # Symplectic propagator for real-time TDDFT 

112 Im_rho_nv += dt * D_n[:, np.newaxis] * Re_rho_nv + 2 * dt * np.dot(K_nn, Re_rho_nv) 

113 Re_rho_nv -= dt * D_n[:, np.newaxis] * Im_rho_nv 

114 

115 # Broadening 

116 broadening = casida_response.broadening 

117 return TimeDomainResponse(time_t, dm_tvv, broadening=broadening, 

118 name=name, units='au')