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

1from typing import Optional 

2import logging 

3 

4import numpy as np 

5 

6from .base import BaseResponse 

7from ..types import Array 

8from ..units import au_to_eV, eV_to_au, au_to_eA, eA_to_au 

9from ..utilities import ClassFormatter 

10from .excitations import Excitations 

11from .utilities import Broadening, NoArtificialBroadening 

12 

13 

14logger = logging.getLogger(__name__) 

15 

16 

17class CasidaResponse(BaseResponse): 

18 """Objects of this class hold different representations of the response 

19 functions for a Casida system. 

20 

21 Parameters 

22 ---------- 

23 D 

24 Diagonal part of decomposition of Casida matrix 

25 (initial-final energy difference of uncoupled single particle levels). 

26 

27 Units are eV by default, optionally atomic units (see :attr:`units`). 

28 K 

29 Coupling matrix of decomposition of Casida matrix 

30 (electron-hole coupling). 

31 

32 Units are eV by default, optionally atomic units (see :attr:`units`). 

33 mu 

34 Transition dipoles. 

35 

36 Units are eÅ by default, optionally atomic units (see :attr:`units`). 

37 occ 

38 Occupation number differences, defaults to one, unitless. 

39 broadening 

40 Artificial broadening used for the continuous response; 

41 defaults to no broadening. 

42 units 

43 `eVA` to specify D, K in eV and mu in eÅ or `au` to specify inputs 

44 in atomic units. 

45 

46 This parameter determines whether conversion should be performed during 

47 initialization and has no effect on instance methods and variables. 

48 name 

49 Name of response functioN. 

50 

51 Examples 

52 -------- 

53 

54 The following snippet illustrates how a Casida response object can be constructed. 

55 Here, artificial Casida matrix components are used. In practice the Casida matrix 

56 would usually be imported from, e.g., NWChem calculations: 

57 

58 >>> import numpy as np 

59 >>> from strongcoca.response import CasidaResponse 

60 >>> 

61 >>> D = np.array([0.1, 0.3]) 

62 >>> K = np.zeros((2, 2)) 

63 >>> mu = np.array([[0.1, 0.0, 0.0], [0.0, 0.1, 0.0]]) 

64 >>> response = CasidaResponse(D, K, mu, name='Panda') 

65 >>> print(response) 

66 ---------------------- CasidaResponse ---------------------- 

67 name : Panda 

68 n_states : 2 

69 D (eV) : [0.1 0.3] 

70 K (eV) : [[0. 0.] 

71 [0. 0.]] 

72 mu (eÅ) : [[0.1 0. 0. ] 

73 [0. 0.1 0. ]] 

74 broadening : NoArtificialBroadening 

75 atoms : None 

76 """ 

77 

78 def __init__(self, 

79 D: np.ndarray, 

80 K: np.ndarray, 

81 mu: np.ndarray, 

82 broadening: Broadening = NoArtificialBroadening(), 

83 occ: Optional[np.ndarray] = None, 

84 units: str = 'eVA', 

85 name: str = 'Casida') -> None: 

86 logger.debug(f'Entering {self.__class__.__name__}.__init__') 

87 super().__init__(broadening=broadening, pbc=False, name=name) 

88 

89 # TODO: explore way to generate these automatically based on type hints 

90 if not isinstance(D, np.ndarray): 

91 raise TypeError(f'D matrix must be provided as numpy array: {D}') 

92 if not isinstance(K, np.ndarray): 

93 raise TypeError(f'K matrix must be provided as numpy array: {K}') 

94 if not isinstance(mu, np.ndarray): 

95 raise TypeError(f'mu matrix must be provided as numpy array: {mu}') 

96 if len(D.shape) != 1: 

97 raise ValueError(f'D matrix has the wrong shape: {D.shape}') 

98 n = len(D) 

99 if K.shape != (n, n): 

100 raise ValueError(f'K matrix must be {n}x{n}; actual shape: {K.shape}') 

101 if mu.shape != (n, 3): 

102 raise ValueError(f'mu matrix must be {n}x{3}; actual shape: {mu.shape}') 

103 if occ is None: 

104 occ = np.ones_like(D) 

105 else: 

106 # TODO: check that all code supports occupation numbers 

107 raise NotImplementedError('Occupation numbers not implemented') 

108 evals = np.linalg.eigvalsh(K) 

109 if not np.all(evals > -1e-8): 

110 raise ValueError(f'K matrix must be positive definite; eigenvalues: {evals}') 

111 

112 D = np.array(D, dtype=float) 

113 K = np.array(K, dtype=float) 

114 mu = np.array(mu, dtype=float) 

115 

116 if units == 'eVA': 

117 D *= eV_to_au 

118 K *= eV_to_au 

119 mu *= eA_to_au 

120 elif units != 'au': 

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

122 

123 self._D_n = D 

124 self._K_nn = K 

125 self._mu_nv = mu 

126 self._f_n = occ 

127 

128 # TODO: do not diagonalize and evaluate the following in init but 

129 # use lazy evaluation 

130 # See also issue #34 

131 

132 # compose and diagonalize Casida matrix 

133 sq_fD_n = np.sqrt(self._f_n * self._D_n) 

134 C_nn = np.diag(self._D_n)**2 + 2 * np.outer(sq_fD_n, sq_fD_n) * self._K_nn 

135 eigval_I, F_nI = np.linalg.eigh(C_nn) 

136 

137 # Calculate transition dipole moments 

138 omega_I = np.sqrt(eigval_I) 

139 mu_Iv = np.einsum('nv,n,nI,I->Iv', self._mu_nv, sq_fD_n, 

140 F_nI, 1. / np.sqrt(omega_I), 

141 optimize=True) 

142 self._excitations = Excitations(omega_I, mu_Iv, broadening=self._broadening, 

143 name=self._name, units='au') 

144 self._eigvec_nI = F_nI 

145 

146 def __str__(self) -> str: 

147 fmt = ClassFormatter(self, pad=15) 

148 fmt.append_class_name() 

149 

150 fmt.append_attr('name') 

151 fmt.append_attr('n_states') 

152 

153 fmt.append_attr('D', unit='eV') 

154 fmt.append_attr('K', unit='eV') 

155 fmt.append_attr('mu', unit='eÅ') 

156 fmt.append_attr('broadening') 

157 formula = self.atoms.get_chemical_formula() if self.atoms is not None else None 

158 fmt.append_attr('atoms', formula) 

159 

160 return fmt.to_string() 

161 

162 @property 

163 def n_states(self) -> int: 

164 """Number of states in Casida system.""" 

165 return len(self._D_n) 

166 

167 @property 

168 def D(self) -> np.ndarray: 

169 """D matrix of Casida system in units of eV.""" 

170 return self._D_n * au_to_eV 

171 

172 @property 

173 def K(self) -> np.ndarray: 

174 """K matrix of Casida system in units of eV.""" 

175 return self._K_nn * au_to_eV 

176 

177 @property 

178 def mu(self) -> np.ndarray: 

179 """mu matrix of Casida system in units of eÅ.""" 

180 return self._mu_nv * au_to_eA 

181 

182 @property 

183 def U(self) -> np.ndarray: 

184 """Unitless eigenvectors of Casida matrix.""" 

185 return self._eigvec_nI 

186 

187 @property 

188 def excitations(self) -> Excitations: 

189 """Excitations as :class:`~strongcoca.response.excitations.Excitations`.""" 

190 return self._excitations 

191 

192 def _get_dynamic_polarizability(self, frequencies: Array) -> np.ndarray: 

193 return self._excitations._get_dynamic_polarizability(frequencies) 

194 

195 def _get_dynamic_polarizability_imaginary_frequency( 

196 self, frequencies: Array) -> np.ndarray: 

197 return self._excitations._get_dynamic_polarizability_imaginary_frequency( 

198 frequencies)