Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : !> \brief Matrix algebra for analytic k derivatives through fifth order.
9 : !> \note Each array stores normalized derivatives X_α=(∂_k^α X)/α!,
10 : !> with α=(α_x,α_y,α_z), α!=α_x!α_y!α_z! and |α|<=5.
11 : !> Examples: X_(1,1,0)=∂_kx ∂_ky X and X_(2,0,0)=(∂_kx² X)/2.
12 : !> Products are formed by the multivariable product rule; no finite differences
13 : !> are used. The overlap inverse square root X=S^(-1/2) satisfies X S X=I.
14 : !> Its derivative of each order is obtained from lower orders by a Sylvester solve.
15 : ! **************************************************************************************************
16 : MODULE nlo_derivatives
17 : USE ieee_arithmetic, ONLY: ieee_is_finite
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: gaussi
20 : USE mathlib, ONLY: diag_complex
21 : #include "../base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nlo_derivatives'
28 :
29 : INTEGER, PARAMETER, PUBLIC :: nlo_max_derivatives_order = 5, &
30 : nlo_n_derivatives = 56
31 : PUBLIC :: nlo_derivative_indices, nlo_derivative_product, nlo_lowdin_derivatives, &
32 : nlo_orthonormal_operator_derivatives, &
33 : nlo_equilibrium_density_derivatives, &
34 : nlo_derivative_count, derivatives_index
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief Number of Cartesian Taylor coefficients through the given total grade.
40 : !> \param grade Highest total derivative grade, from zero through five.
41 : !> \return Number of Cartesian coefficients through the requested grade.
42 : ! **************************************************************************************************
43 1049 : INTEGER FUNCTION nlo_derivative_count(grade) RESULT(count)
44 : INTEGER, INTENT(IN) :: grade
45 :
46 1049 : count = 0
47 1049 : IF (grade < 0 .OR. grade > nlo_max_derivatives_order) THEN
48 : RETURN
49 : END IF
50 1049 : count = (grade + 1)*(grade + 2)*(grade + 3)/6
51 :
52 1049 : END FUNCTION nlo_derivative_count
53 :
54 : ! **************************************************************************************************
55 : !> \brief Enumerate the 56 three-dimensional multi-indices through total order five.
56 : !> \param alpha Cartesian derivative counts (α_x, α_y, α_z) for each Taylor coefficient.
57 : ! **************************************************************************************************
58 215 : SUBROUTINE nlo_derivative_indices(alpha)
59 : INTEGER, INTENT(OUT) :: alpha(3, nlo_n_derivatives)
60 :
61 : INTEGER :: ia, ib, ij, iorder
62 :
63 215 : ij = 0
64 1505 : DO iorder = 0, nlo_max_derivatives_order
65 6020 : DO ia = iorder, 0, -1
66 17845 : DO ib = iorder - ia, 0, -1
67 12040 : ij = ij + 1
68 52675 : alpha(:, ij) = [ia, ib, iorder - ia - ib]
69 : END DO
70 : END DO
71 : END DO
72 :
73 215 : END SUBROUTINE nlo_derivative_indices
74 :
75 : ! **************************************************************************************************
76 : !> \brief Apply the product rule to two arrays of normalized matrix k derivatives.
77 : !> \param alpha Cartesian derivative counts (α_x, α_y, α_z) for each Taylor coefficient.
78 : !> \param left Normalized Taylor coefficients of the left matrix factor.
79 : !> \param right Normalized Taylor coefficients of the right matrix factor.
80 : !> \param product Normalized Taylor coefficients of the matrix product.
81 : !> \param active_grade Highest total derivative grade to compute; omit for all five grades.
82 : ! **************************************************************************************************
83 168 : SUBROUTINE nlo_derivative_product(alpha, left, right, product, active_grade)
84 : INTEGER, INTENT(IN) :: alpha(3, nlo_n_derivatives)
85 : COMPLEX(dp), INTENT(IN) :: left(:, :, :), right(:, :, :)
86 : COMPLEX(dp), INTENT(OUT) :: product(:, :, :)
87 : INTEGER, INTENT(IN), OPTIONAL :: active_grade
88 :
89 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_derivative_product'
90 :
91 168 : COMPLEX(dp), ALLOCATABLE :: term(:, :)
92 : INTEGER :: handle, ib, ij, ik, nj
93 :
94 168 : CALL timeset(routineN, handle)
95 :
96 168 : nj = nlo_n_derivatives
97 168 : IF (PRESENT(active_grade)) nj = nlo_derivative_count(active_grade)
98 1787046 : product = 0.0_dp
99 672 : ALLOCATE (term(SIZE(left, 1), SIZE(right, 2)))
100 2894 : DO ij = 1, nj
101 43320 : DO ib = 1, ij
102 161704 : ik = derivatives_index(alpha, alpha(:, ij) - alpha(:, ib))
103 40426 : IF (ik == 0) CYCLE
104 176558732 : term(:, :) = MATMUL(left(:, :, ib), right(:, :, ik))
105 7106578 : PRODUCT(:, :, ij) = PRODUCT(:, :, ij) + term
106 : END DO
107 : END DO
108 :
109 168 : CALL timestop(handle)
110 :
111 336 : END SUBROUTINE nlo_derivative_product
112 :
113 : ! **************************************************************************************************
114 : !> \brief Solve X(k) S(k) X(k)=I analytically through fifth Cartesian order.
115 : !> \param overlap Normalized Taylor coefficients of S; leading coefficient Hermitian positive
116 : !> \param inverse_root Normalized Taylor coefficients of S**(-1/2)
117 : !> \param active_grade Highest total derivative grade to compute; omit for all five grades.
118 : ! **************************************************************************************************
119 12 : SUBROUTINE nlo_lowdin_derivatives(overlap, inverse_root, active_grade)
120 : COMPLEX(dp), INTENT(IN) :: overlap(:, :, :)
121 : COMPLEX(dp), INTENT(OUT) :: inverse_root(:, :, :)
122 : INTEGER, INTENT(IN), OPTIONAL :: active_grade
123 :
124 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_lowdin_derivatives'
125 :
126 12 : COMPLEX(dp), ALLOCATABLE :: matrix_work(:, :), metric_u(:, :), &
127 : rhs(:, :)
128 12 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: matrix_product
129 : INTEGER :: handle, ib, ic, id, ij, n, nj
130 : INTEGER, DIMENSION(3, nlo_n_derivatives) :: alpha
131 : REAL(KIND=dp) :: scale
132 12 : REAL(KIND=dp), ALLOCATABLE :: metric_e(:)
133 :
134 12 : CALL timeset(routineN, handle)
135 :
136 12 : n = SIZE(overlap, 1)
137 157299 : inverse_root = 0.0_dp
138 12 : nj = nlo_n_derivatives
139 12 : IF (PRESENT(active_grade)) nj = nlo_derivative_count(active_grade)
140 48 : CPASSERT(.NOT. (n < 1 .OR. ANY(SHAPE(overlap) /= [n, n, nj])))
141 48 : CPASSERT(.NOT. (ANY(SHAPE(inverse_root) /= [n, n, nj])))
142 12 : CPASSERT(nj >= 1)
143 325 : DO ij = 1, nj
144 157287 : scale = MAX(1.0_dp, MAXVAL(ABS(overlap(:, :, ij))))
145 : IF (.NOT. ALL(ieee_is_finite(REAL(overlap(:, :, ij), KIND=dp))) .OR. &
146 471861 : .NOT. ALL(ieee_is_finite(AIMAG(overlap(:, :, ij)))) .OR. &
147 : MAXVAL(ABS(overlap(:, :, ij) - &
148 12 : CONJG(TRANSPOSE(overlap(:, :, ij))))) > 1.e-10_dp*scale) THEN
149 0 : CPABORT('NONLINEAR_OPTICS Löwdin derivative check failed.')
150 : END IF
151 : END DO
152 : ALLOCATE (metric_u(n, n), rhs(n, n), matrix_work(n, n), &
153 144 : matrix_product(n, n), metric_e(n))
154 :
155 12 : CALL diag_complex(overlap(:, :, 1), metric_u, metric_e)
156 :
157 508 : IF (MAXVAL(metric_e) <= 0.0_dp .OR. MINVAL(metric_e) <= 1.e-10_dp*MAXVAL(metric_e)) THEN
158 0 : CPABORT('NONLINEAR_OPTICS Löwdin derivative check failed.')
159 : END IF
160 12 : rhs = 0.0_dp
161 260 : DO ib = 1, n
162 260 : rhs(ib, ib) = 1.0_dp/SQRT(metric_e(ib))
163 : END DO
164 295284 : matrix_work(:, :) = MATMUL(rhs, CONJG(TRANSPOSE(metric_u)))
165 301340 : inverse_root(:, :, 1) = MATMUL(metric_u, matrix_work)
166 12 : CALL nlo_derivative_indices(alpha)
167 :
168 313 : DO ij = 2, nj
169 301 : rhs = 0.0_dp
170 6612 : DO ib = 1, ij
171 212393 : DO ic = 1, ij
172 823124 : id = derivatives_index(alpha, alpha(:, ij) - alpha(:, ib) - alpha(:, ic))
173 205781 : IF (id == 0) CYCLE
174 7644 : IF (ib == ij .OR. id == ij) CYCLE
175 89423166 : matrix_work(:, :) = MATMUL(overlap(:, :, ic), inverse_root(:, :, id))
176 89423166 : matrix_product(:, :) = MATMUL(inverse_root(:, :, ib), matrix_work)
177 3800488 : rhs(:, :) = rhs - matrix_product
178 : END DO
179 : END DO
180 7350467 : matrix_work(:, :) = MATMUL(rhs, metric_u)
181 7350467 : rhs(:, :) = MATMUL(CONJG(TRANSPOSE(metric_u)), matrix_work)
182 6495 : DO ib = 1, n
183 151219 : DO ic = 1, n
184 150918 : rhs(ib, ic) = rhs(ib, ic)/(SQRT(metric_e(ib)) + SQRT(metric_e(ic)))
185 : END DO
186 : END DO
187 7350467 : matrix_work(:, :) = MATMUL(rhs, CONJG(TRANSPOSE(metric_u)))
188 7501385 : inverse_root(:, :, ij) = MATMUL(metric_u, matrix_work)
189 151219 : scale = MAX(1.0_dp, MAXVAL(ABS(inverse_root(:, :, ij))))
190 151219 : IF (MAXVAL(ABS(inverse_root(:, :, ij) - &
191 12 : CONJG(TRANSPOSE(inverse_root(:, :, ij))))) > 1.e-9_dp*scale) THEN
192 0 : CPABORT('NONLINEAR_OPTICS Löwdin derivative check failed.')
193 : END IF
194 : END DO
195 :
196 12 : CALL timestop(handle)
197 :
198 24 : END SUBROUTINE nlo_lowdin_derivatives
199 :
200 : ! **************************************************************************************************
201 : !> \brief Express H and the AO position matrix in an orthonormal AO basis.
202 : !> \param overlap Normalized Taylor coefficients of the AO overlap S(k).
203 : !> \param hamiltonian Normalized Taylor coefficients of the AO Hamiltonian H(k).
204 : !> \param moment Normalized Taylor coefficients of the AO position matrices d_a(k).
205 : !> \param orth_h Normalized Taylor coefficients of the orthogonal Hamiltonian.
206 : !> \param connection Normalized Taylor coefficients of the orthogonal position connection.
207 : !> \param active_grade Highest total derivative grade to compute; omit for all five grades.
208 : !> \note AO basis functions overlap, so S(k) is not the identity. The matrix
209 : !> X(k)=S(k)^(-1/2) makes their overlap equal to the identity.
210 : !> This routine differentiates H^orth=X H X and the position matrix
211 : !> A_a^orth=X d_a X-i (dX/dk_a) S X, where d_a is the AO position matrix.
212 : !> It does not select occupied bands or build the equilibrium density rho0.
213 : !> A_a^orth is needed through order four; its order-five slots remain zero.
214 : ! **************************************************************************************************
215 24 : SUBROUTINE nlo_orthonormal_operator_derivatives(overlap, hamiltonian, moment, orth_h, &
216 12 : connection, active_grade)
217 : COMPLEX(dp), DIMENSION(:, :, :), INTENT(IN) :: overlap, hamiltonian
218 : COMPLEX(dp), INTENT(IN) :: moment(:, :, :, :)
219 : COMPLEX(dp), INTENT(OUT) :: orth_h(:, :, :), connection(:, :, :, :)
220 : INTEGER, INTENT(IN), OPTIONAL :: active_grade
221 :
222 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_orthonormal_operator_derivatives'
223 :
224 : COMPLEX(dp), ALLOCATABLE :: derivative_x(:, :, :), &
225 12 : inverse_root(:, :, :), &
226 12 : matrix_derivatives(:, :, :), &
227 : moment_derivatives(:, :, :)
228 : INTEGER :: alpha(3, nlo_n_derivatives), grade, &
229 : handle, ia, ij, ik, n, na, nh, &
230 : sought(3)
231 : REAL(KIND=dp) :: scale
232 :
233 12 : CALL timeset(routineN, handle)
234 :
235 : ! AO overlap and Hamiltonian.
236 12 : n = SIZE(overlap, 1)
237 157299 : orth_h = 0.0_dp
238 264478 : connection = 0.0_dp
239 12 : grade = nlo_max_derivatives_order
240 12 : IF (PRESENT(active_grade)) grade = active_grade
241 12 : nh = nlo_derivative_count(grade)
242 12 : na = nlo_derivative_count(grade - 1)
243 12 : CPASSERT(grade >= 2 .AND. nh >= 1)
244 48 : CPASSERT(ALL(SHAPE(overlap) == [n, n, nh]))
245 48 : CPASSERT(ALL(SHAPE(hamiltonian) == [n, n, nh]))
246 60 : CPASSERT(ALL(SHAPE(moment) == [n, n, 3, na]))
247 48 : CPASSERT(ALL(SHAPE(orth_h) == [n, n, nh]))
248 60 : CPASSERT(ALL(SHAPE(connection) == [n, n, 3, na]))
249 : IF (.NOT. ALL(ieee_is_finite(REAL(hamiltonian, KIND=dp))) .OR. &
250 : .NOT. ALL(ieee_is_finite(AIMAG(hamiltonian))) .OR. &
251 843554 : .NOT. ALL(ieee_is_finite(REAL(moment, KIND=dp))) .OR. &
252 : .NOT. ALL(ieee_is_finite(AIMAG(moment)))) THEN
253 0 : CPABORT('NONLINEAR_OPTICS AO operator derivatives are inconsistent.')
254 : END IF
255 12 : CALL nlo_derivative_indices(alpha)
256 325 : DO ij = 1, nh
257 157287 : scale = MAX(1.0_dp, MAXVAL(ABS(hamiltonian(:, :, ij))))
258 157287 : IF (MAXVAL(ABS(hamiltonian(:, :, ij) - &
259 : CONJG(TRANSPOSE(hamiltonian(:, :, ij))))) > 1.e-9_dp*scale) THEN
260 0 : CPABORT('NONLINEAR_OPTICS AO operator derivatives are inconsistent.')
261 : END IF
262 1252 : IF (SUM(alpha(:, ij)) >= grade) CYCLE
263 712 : DO ia = 1, 3
264 2100 : sought = alpha(:, ij)
265 525 : sought(ia) = sought(ia) + 1
266 525 : ik = derivatives_index(alpha, sought)
267 : scale = MAX(1.0_dp, MAXVAL(ABS(moment(:, :, ia, ij))), &
268 528057 : MAXVAL(ABS(overlap(:, :, ik))))
269 264291 : IF (MAXVAL(ABS(moment(:, :, ia, ij) - &
270 : CONJG(TRANSPOSE(moment(:, :, ia, ij))) + &
271 : gaussi*REAL(sought(ia), KIND=dp)*overlap(:, :, ik))) > &
272 313 : 1.e-9_dp*scale) THEN
273 0 : CPABORT('NONLINEAR_OPTICS AO operator derivatives are inconsistent.')
274 : END IF
275 : END DO
276 : END DO
277 : ALLOCATE (inverse_root(n, n, nh), derivative_x(n, n, na), &
278 180 : matrix_derivatives(n, n, nh), moment_derivatives(n, n, na))
279 :
280 : ! X = S^(-1/2), with X S X = 1 at every derivative order.
281 12 : CALL nlo_lowdin_derivatives(overlap, inverse_root, active_grade=grade)
282 :
283 : ! H^o = X H X; coefficient products use the multivariable product rule.
284 : CALL nlo_derivative_product(alpha, inverse_root, hamiltonian, matrix_derivatives, &
285 12 : active_grade=grade)
286 : CALL nlo_derivative_product(alpha, matrix_derivatives, inverse_root, orth_h, &
287 12 : active_grade=grade)
288 :
289 48 : DO ia = 1, 3
290 36 : derivative_x = 0.0_dp
291 561 : DO ij = 1, na
292 2100 : IF (SUM(alpha(:, ij)) >= grade) CYCLE
293 2100 : sought = alpha(:, ij)
294 525 : sought(ia) = sought(ia) + 1
295 525 : ik = derivatives_index(alpha, sought)
296 264327 : derivative_x(:, :, ij) = REAL(alpha(ia, ij) + 1, KIND=dp)*inverse_root(:, :, ik)
297 : END DO
298 : CALL nlo_derivative_product(alpha, derivative_x, overlap, matrix_derivatives, &
299 36 : active_grade=grade - 1)
300 : CALL nlo_derivative_product(alpha, matrix_derivatives, inverse_root, &
301 36 : moment_derivatives, active_grade=grade - 1)
302 264327 : connection(:, :, ia, :) = -gaussi*moment_derivatives
303 : CALL nlo_derivative_product(alpha, inverse_root, moment(:, :, ia, :), &
304 36 : matrix_derivatives, active_grade=grade - 1)
305 : CALL nlo_derivative_product(alpha, matrix_derivatives, inverse_root, &
306 36 : moment_derivatives, active_grade=grade - 1)
307 264327 : connection(:, :, ia, :) = connection(:, :, ia, :) + moment_derivatives
308 573 : DO ij = 1, na
309 264291 : scale = MAX(1.0_dp, MAXVAL(ABS(connection(:, :, ia, ij))))
310 264291 : IF (MAXVAL(ABS(connection(:, :, ia, ij) - &
311 36 : CONJG(TRANSPOSE(connection(:, :, ia, ij))))) > 1.e-9_dp*scale) THEN
312 0 : CPABORT('NONLINEAR_OPTICS AO operator derivatives are inconsistent.')
313 : END IF
314 : END DO
315 : END DO
316 :
317 12 : CALL timestop(handle)
318 :
319 24 : END SUBROUTINE nlo_orthonormal_operator_derivatives
320 :
321 : ! **************************************************************************************************
322 : !> \brief Derivatives of the equilibrium density matrix for filled Kohn-Sham bands.
323 : !> \param energies Band energies ε_m at the central k point.
324 : !> \param occupied Unit occupation flags for the bands at the central k point.
325 : !> \param hamiltonian Normalized Taylor coefficients of the AO Hamiltonian H(k).
326 : !> \param rho0 Normalized Taylor coefficients of the equilibrium density matrix.
327 : !> \param active_grade Highest total derivative grade to compute; omit for all five grades.
328 : !> \note At the central k point, rho0_mn = occupied_m delta_mn with occupied_m=0 or 1.
329 : !> At nearby k, rho0(k) = sum_(m occupied) |u_m(k)><u_m(k)|.
330 : !> This is a spectral projector because rho0**2 = rho0; it is also the
331 : !> equilibrium density matrix for unit band occupation. The actual spin
332 : !> occupation factor is applied later during Brillouin-zone integration.
333 : !> The routine obtains analytic k derivatives from H rho0 = rho0 H
334 : !> and rho0**2 = rho0, without differentiating individual eigenvectors.
335 : !> It does not calculate the field-induced rho^(1), rho^(2), etc.
336 : ! **************************************************************************************************
337 12 : SUBROUTINE nlo_equilibrium_density_derivatives(energies, occupied, hamiltonian, rho0, &
338 : active_grade)
339 : REAL(KIND=dp), INTENT(IN) :: energies(:)
340 : INTEGER, INTENT(IN) :: occupied(:)
341 : COMPLEX(dp), INTENT(IN) :: hamiltonian(:, :, :)
342 : COMPLEX(dp), INTENT(OUT) :: rho0(:, :, :)
343 : INTEGER, INTENT(IN), OPTIONAL :: active_grade
344 :
345 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_equilibrium_density_derivatives'
346 :
347 12 : COMPLEX(dp), ALLOCATABLE :: block(:, :), forcing(:, :), &
348 12 : matrix_term(:, :)
349 : INTEGER :: handle, ib, ic, ij, ik, n, nj
350 : INTEGER, DIMENSION(3, nlo_n_derivatives) :: alpha
351 : REAL(KIND=dp) :: gap, scale, symmetry_error
352 :
353 12 : CALL timeset(routineN, handle)
354 :
355 12 : n = SIZE(energies)
356 157299 : rho0 = 0.0_dp
357 12 : CPASSERT(.NOT. (n < 2 .OR. SIZE(occupied) /= n))
358 12 : nj = nlo_n_derivatives
359 12 : IF (PRESENT(active_grade)) nj = nlo_derivative_count(active_grade)
360 48 : CPASSERT(ALL(SHAPE(hamiltonian) == [n, n, nj]))
361 48 : CPASSERT(ALL(SHAPE(rho0) == [n, n, nj]))
362 520 : CPASSERT(.NOT. (ANY(occupied < 0) .OR. ANY(occupied > 1)))
363 68 : CPASSERT(.NOT. (ALL(occupied == 0) .OR. ALL(occupied == 1)))
364 : IF (.NOT. ALL(ieee_is_finite(energies)) .OR. &
365 314858 : .NOT. ALL(ieee_is_finite(REAL(hamiltonian, KIND=dp))) .OR. &
366 : .NOT. ALL(ieee_is_finite(AIMAG(hamiltonian)))) THEN
367 0 : CPABORT('NONLINEAR_OPTICS equilibrium density input is not finite.')
368 : END IF
369 260 : scale = MAX(1.0_dp, MAXVAL(ABS(energies)))
370 : ! The occupied and empty energies must remain separated at this k point.
371 508 : gap = MINVAL(energies, MASK=occupied == 0) - MAXVAL(energies, MASK=occupied == 1)
372 12 : IF (gap <= 1.e-10_dp*scale) THEN
373 0 : CPABORT('NONLINEAR_OPTICS requires a positive occupied-empty gap.')
374 : END IF
375 96 : ALLOCATE (block(n, n), forcing(n, n), matrix_term(n, n))
376 6068 : matrix_term(:, :) = hamiltonian(:, :, 1)
377 260 : DO ib = 1, n
378 260 : matrix_term(ib, ib) = matrix_term(ib, ib) - energies(ib)
379 : END DO
380 6068 : IF (MAXVAL(ABS(matrix_term)) > 1.e-9_dp*scale) THEN
381 0 : CPABORT('NONLINEAR_OPTICS Hamiltonian is not diagonal in the band basis.')
382 : END IF
383 12 : CALL nlo_derivative_indices(alpha)
384 : ! Zeroth Taylor coefficient: the 0/1 equilibrium density in the band basis.
385 260 : DO ib = 1, n
386 260 : rho0(ib, ib, 1) = occupied(ib)
387 : END DO
388 : ! Differentiate [H,rho0]=0 and rho0**2=rho0 order by order in k.
389 313 : DO ij = 2, nj
390 301 : block = 0.0_dp
391 301 : forcing = 0.0_dp
392 6311 : DO ik = 2, ij
393 24040 : ic = derivatives_index(alpha, alpha(:, ij) - alpha(:, ik))
394 6010 : IF (ic == 0) CYCLE
395 21140589 : matrix_term(:, :) = MATMUL(hamiltonian(:, :, ik), rho0(:, :, ic))
396 850741 : forcing(:, :) = forcing + matrix_term
397 21140589 : matrix_term(:, :) = MATMUL(rho0(:, :, ic), hamiltonian(:, :, ik))
398 850741 : forcing(:, :) = forcing - matrix_term
399 1675 : IF (ic == 1) CYCLE
400 17389746 : matrix_term(:, :) = MATMUL(rho0(:, :, ik), rho0(:, :, ic))
401 704459 : block(:, :) = block + matrix_term
402 : END DO
403 6507 : DO ic = 1, n
404 151219 : DO ib = 1, n
405 : ! Cross-gap blocks use the energy difference; same-group blocks
406 : ! follow from rho0**2 = rho0.
407 150918 : IF (occupied(ib) /= occupied(ic)) THEN
408 43626 : rho0(ib, ic, ij) = -forcing(ib, ic)/(energies(ib) - energies(ic))
409 101098 : ELSE IF (occupied(ib) == 1) THEN
410 5077 : rho0(ib, ic, ij) = -block(ib, ic)
411 : ELSE
412 96021 : rho0(ib, ic, ij) = block(ib, ic)
413 : END IF
414 : END DO
415 : END DO
416 : END DO
417 : symmetry_error = 0.0_dp
418 325 : DO ij = 1, nj
419 : symmetry_error = MAX(symmetry_error, &
420 : MAXVAL(ABS(rho0(:, :, ij) - CONJG(TRANSPOSE(rho0(:, :, ij)))))/ &
421 314586 : MAX(1.0_dp, MAXVAL(ABS(rho0(:, :, ij)))))
422 : END DO
423 : IF (.NOT. ALL(ieee_is_finite(REAL(rho0, KIND=dp))) .OR. &
424 314598 : .NOT. ALL(ieee_is_finite(AIMAG(rho0))) .OR. symmetry_error > 1.e-8_dp) THEN
425 0 : CPABORT('NONLINEAR_OPTICS equilibrium density derivatives are invalid.')
426 : END IF
427 :
428 12 : CALL timestop(handle)
429 :
430 24 : END SUBROUTINE nlo_equilibrium_density_derivatives
431 :
432 : ! **************************************************************************************************
433 : !> \brief Find the Taylor coefficient for the requested (x,y,z) derivative counts.
434 : !> \param alpha Cartesian derivative counts (α_x, α_y, α_z) for each Taylor coefficient.
435 : !> \param sought Requested Cartesian derivative counts; negative counts return zero.
436 : !> \return Its index in alpha, or zero when any count is negative or absent.
437 : ! **************************************************************************************************
438 262315 : INTEGER FUNCTION derivatives_index(alpha, sought) RESULT(index)
439 : INTEGER, INTENT(IN) :: alpha(3, nlo_n_derivatives), sought(3)
440 :
441 : INTEGER :: ij
442 :
443 262315 : index = 0
444 454678 : IF (ANY(sought < 0)) RETURN
445 222190 : DO ij = 1, nlo_n_derivatives
446 384436 : IF (ALL(alpha(:, ij) == sought)) THEN
447 262315 : index = ij
448 : RETURN
449 : END IF
450 : END DO
451 :
452 : END FUNCTION derivatives_index
453 :
454 : END MODULE nlo_derivatives
|
