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 if units == 'eVA': 

113 D = D * eV_to_au 

114 K = K * eV_to_au 

115 mu = mu * eA_to_au 

116 elif units != 'au': 

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

118 

119 self._D_n = D 

120 self._K_nn = K 

121 self._mu_nv = mu 

122 self._f_n = occ 

123 

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

125 # use lazy evaluation 

126 # See also issue #34 

127 

128 # compose and diagonalize Casida matrix 

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

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

131 eigval_I, F_nI = np.linalg.eig(C_nn) 

132 

133 # Calculate transition dipole moments 

134 omega_I = np.sqrt(eigval_I) 

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

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

137 optimize=True) 

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

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

140 self._eigvec_nI = F_nI 

141 

142 def __str__(self) -> str: 

143 fmt = ClassFormatter(self, pad=15) 

144 fmt.append_class_name() 

145 

146 fmt.append_attr('name') 

147 fmt.append_attr('n_states') 

148 

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

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

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

152 fmt.append_attr('broadening') 

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

154 fmt.append_attr('atoms', formula) 

155 

156 return fmt.to_string() 

157 

158 @property 

159 def n_states(self) -> int: 

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

161 return len(self._D_n) 

162 

163 @property 

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

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

166 return self._D_n * au_to_eV # type: ignore 

167 

168 @property 

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

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

171 return self._K_nn * au_to_eV # type: ignore 

172 

173 @property 

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

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

176 return self._mu_nv * au_to_eA # type: ignore 

177 

178 @property 

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

180 """Unitless eigenvectors of Casida matrix.""" 

181 return self._eigvec_nI # type: ignore 

182 

183 @property 

184 def excitations(self) -> Excitations: 

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

186 return self._excitations 

187 

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

189 return self._excitations._get_dynamic_polarizability(frequencies) 

190 

191 def _get_dynamic_polarizability_imaginary_frequency( 

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

193 return self._excitations._get_dynamic_polarizability_imaginary_frequency( 

194 frequencies)