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 Local full-tensor independent-particle current response through fifth field order.
9 : !> \note This is a Bloch-point kernel. It does not integrate the Brillouin zone,
10 : !> convert current to susceptibility, or build AO derivatives from CP2K data.
11 : ! **************************************************************************************************
12 : MODULE nlo_current
13 : USE ieee_arithmetic, ONLY: ieee_is_finite
14 : USE kinds, ONLY: dp
15 : USE mathconstants, ONLY: gaussi
16 : USE nlo_T1_T2_band_transition, ONLY: nlo_band_transition_correction,&
17 : nlo_band_transition_prepare,&
18 : nlo_band_transition_rate,&
19 : nlo_band_transition_type
20 : USE nlo_derivatives, ONLY: derivatives_index,&
21 : nlo_derivative_count,&
22 : nlo_derivative_indices,&
23 : nlo_n_derivatives
24 : #include "../base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nlo_current'
31 :
32 : PUBLIC :: nlo_current_tensor_derivatives
33 :
34 : CONTAINS
35 :
36 : ! **************************************************************************************************
37 : !> \brief Compute intrinsically permutation-symmetric χ^(n) current at one k point.
38 : !> \param order Number of incident fields, 2:5
39 : !> \param omega Incident photon energies in Hartree
40 : !> \param energies Central band energies in Hartree
41 : !> \param rate1 Population damping energy in Hartree
42 : !> \param rate2 Coherence damping energy in Hartree
43 : !> \param transition_width Band-energy width of their smooth transition in Hartree
44 : !> \param hamiltonian Orthogonal Bloch Hamiltonian normalized k derivatives through grade n
45 : !> \param connection Orthogonal Bloch position connection derivatives through grade n-1
46 : !> \param curvature Three position-curvature matrices for this k point and spin
47 : !> \param rho0 Equilibrium density derivatives through grade n
48 : !> \param velocity Physical central-band velocity in bohr/atomic time
49 : !> \param current Complex current tensor, 3 by 3**order
50 : !> \param elements Optional output/input selection, same shape and indexing as current.
51 : !> Only required subset responses are evaluated. Unselected outputs remain zero.
52 : !> \note Label the n incident fields by j=1,...,n, each with frequency ω_j and
53 : !> Cartesian direction a_j. For every nonempty labeled subset S, form
54 : !> Ω_S=Σ_{j∈S}ω_j and the normalized density response
55 : !> ρ_S = |S|⁻¹ Σ_(j in S) R(Ω_S)L_(a_j)ρ_(S without j),
56 : !> where R L is the analytic normalized k-derivative step of nlo_density_step_derivatives.
57 : !> This average equals the sum of all ordered responses divided by n!.
58 : !> Only k-derivative grades n-|S| are retained for each subset. States
59 : !> with the same labeled subset and field directions are reused across
60 : !> all full-tensor components, retaining only two adjacent field orders
61 : !> in memory. The final current is -Tr[v_c ρ_{1,...,n}]
62 : !> minus n⁻¹ Σ_j Tr[F_(c,a_j) ρ_(S without j)], in atomic units.
63 : !> This follows from v_c(E)=i[H+E_a r_a,r_c]=v_c(0)+E_a F_ca,
64 : !> with [r_c,r_a]=iF_ca. The curvature term is always included. The second
65 : !> array index stores a_1,...,a_n as 1+Σ_j(a_j-1)3^(j-1).
66 : ! **************************************************************************************************
67 13 : SUBROUTINE nlo_current_tensor_derivatives(order, omega, energies, rate1, rate2, &
68 13 : transition_width, hamiltonian, &
69 26 : connection, curvature, rho0, velocity, current, &
70 13 : elements)
71 : INTEGER, INTENT(IN) :: order
72 : REAL(KIND=dp), INTENT(IN) :: omega(:), energies(:), rate1, rate2, &
73 : transition_width
74 : COMPLEX(dp), INTENT(IN) :: hamiltonian(:, :, :), &
75 : connection(:, :, :, :), &
76 : curvature(:, :, :), rho0(:, :, :), &
77 : velocity(:, :, :)
78 : COMPLEX(dp), INTENT(OUT) :: current(:, :)
79 : LOGICAL, INTENT(IN), OPTIONAL :: elements(:, :)
80 :
81 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_current_tensor_derivatives'
82 :
83 : COMPLEX(dp) :: curvature_term, extra(3, 243)
84 13 : COMPLEX(dp), ALLOCATABLE :: next_density(:, :, :), &
85 13 : previous_density(:, :, :)
86 13 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :, :) :: density
87 : INTEGER :: code, current_start, direction, field, full_mask, handle, i_code, i_field, &
88 : i_mask, i_radix, i_subset_code, index, iother, iother_direction, iother_position, &
89 : ioutput, mask, multiplicity, n, na, ncodes(0:31), ncoeff(0:31), nh, &
90 : offset_current(0:242, 0:31), offset_previous(0:242, 0:31), order_subset, position, &
91 : previous_code, previous_end, previous_mask, previous_start, radix, &
92 : representative_code(0:242, 0:31), total_coefficients
93 : LOGICAL :: evaluate(0:242, 0:31), &
94 : required(0:242, 0:31), selected(3, 243)
95 : REAL(KIND=dp) :: subset_omega, weight
96 :
97 13 : CALL timeset(routineN, handle)
98 :
99 3793 : current = 0.0_dp
100 13 : n = SIZE(energies)
101 13 : CPASSERT(.NOT. (order < 2 .OR. order > 5 .OR. n < 2))
102 13 : CPASSERT(SIZE(omega) == order)
103 13 : nh = nlo_derivative_count(order)
104 13 : na = nlo_derivative_count(order - 1)
105 52 : CPASSERT(ALL(SHAPE(hamiltonian) == [n, n, nh]))
106 65 : CPASSERT(ALL(SHAPE(connection) == [n, n, 3, na]))
107 52 : CPASSERT(ALL(SHAPE(rho0) == [n, n, nh]))
108 52 : CPASSERT(ALL(SHAPE(curvature) == [n, n, 3]))
109 52 : CPASSERT(ALL(SHAPE(velocity) == [n, n, 3]))
110 39 : CPASSERT(ALL(SHAPE(current) == [3, 3**order]))
111 12649 : selected = .TRUE.
112 13 : IF (PRESENT(elements)) THEN
113 39 : CPASSERT(ALL(SHAPE(elements) == SHAPE(current)))
114 3457 : IF (.NOT. ANY(elements)) THEN
115 0 : CALL timestop(handle)
116 0 : RETURN
117 : END IF
118 3793 : selected(:, 1:3**order) = elements
119 : END IF
120 : ! Restrict every labeled subset to the input directions required by selected elements.
121 13 : required = .FALSE.
122 13 : required(0, 0) = .TRUE.
123 958 : DO i_code = 0, 3**order - 1
124 3531 : IF (.NOT. ANY(selected(:, i_code + 1))) CYCLE
125 573 : DO i_mask = 1, 2**order - 1
126 : i_subset_code = 0
127 : i_radix = 1
128 1931 : DO i_field = 1, order
129 1458 : IF (.NOT. BTEST(i_mask, i_field - 1)) CYCLE
130 836 : i_subset_code = i_subset_code + MOD(i_code/3**(i_field - 1), 3)*i_radix
131 1931 : i_radix = 3*i_radix
132 : END DO
133 1418 : required(i_subset_code, i_mask) = .TRUE.
134 : END DO
135 : END DO
136 : IF (.NOT. ALL(ieee_is_finite(omega)) .OR. .NOT. ALL(ieee_is_finite(energies)) .OR. &
137 : .NOT. ieee_is_finite(rate1) .OR. .NOT. ieee_is_finite(rate2) .OR. &
138 : .NOT. ieee_is_finite(transition_width) .OR. &
139 : .NOT. ALL(ieee_is_finite(REAL(velocity, KIND=dp))) .OR. &
140 : .NOT. ALL(ieee_is_finite(AIMAG(velocity))) .OR. &
141 : .NOT. ALL(ieee_is_finite(REAL(hamiltonian, KIND=dp))) .OR. &
142 : .NOT. ALL(ieee_is_finite(AIMAG(hamiltonian))) .OR. &
143 : .NOT. ALL(ieee_is_finite(REAL(connection, KIND=dp))) .OR. &
144 : .NOT. ALL(ieee_is_finite(AIMAG(connection))) .OR. &
145 : .NOT. ALL(ieee_is_finite(REAL(rho0, KIND=dp))) .OR. &
146 : .NOT. ALL(ieee_is_finite(AIMAG(rho0))) .OR. rate1 <= 0.0_dp .OR. &
147 : rate2 < rate1 .OR. transition_width <= 0.0_dp &
148 1244203 : .OR. ABS(SUM(omega)) <= 1.e-14_dp) THEN
149 0 : CPABORT('NLO current response contains invalid data.')
150 : END IF
151 13 : full_mask = 2**order - 1
152 173 : DO mask = 0, full_mask
153 160 : ncoeff(mask) = nlo_derivative_count(order - POPCNT(mask))
154 173 : ncodes(mask) = 3**POPCNT(mask)
155 : END DO
156 : ALLOCATE (previous_density(n, n, ncoeff(0)), &
157 117 : next_density(n, n, nlo_derivative_count(order - 1)))
158 164330 : previous_density(:, :, :) = rho0(:, :, 1:ncoeff(0))
159 13 : offset_previous = 0
160 13 : offset_previous(0, 0) = 1
161 13 : extra = 0.0_dp
162 53 : DO order_subset = 1, order
163 : ! Omega_S = sum_(j in S) omega_j; retain each required field subset S.
164 : CALL subset_density_offsets(order, order_subset, omega, required, ncoeff, &
165 : offset_current, evaluate, representative_code, &
166 40 : total_coefficients)
167 200 : ALLOCATE (density(n, n, total_coefficients))
168 40 : density = 0.0_dp
169 664 : DO mask = 1, full_mask
170 624 : IF (POPCNT(mask) /= order_subset) CYCLE
171 147 : subset_omega = 0.0_dp
172 771 : DO field = 1, order
173 771 : IF (BTEST(mask, field - 1)) subset_omega = subset_omega + omega(field)
174 : END DO
175 3790 : DO code = 0, ncodes(mask) - 1
176 3603 : IF (.NOT. evaluate(code, mask)) CYCLE
177 106 : current_start = offset_current(code, mask)
178 1028 : DO field = 1, order
179 298 : IF (.NOT. BTEST(mask, field - 1)) CYCLE
180 214 : position = POPCNT(IAND(mask, 2**(field - 1) - 1))
181 214 : radix = 3**position
182 214 : direction = MOD(code/radix, 3) + 1
183 214 : multiplicity = 1
184 614 : DO iother = 1, order
185 469 : IF (iother == field .OR. .NOT. BTEST(mask, iother - 1)) CYCLE
186 231 : IF (ABS(omega(iother) - omega(field)) > 0.0_dp) CYCLE
187 207 : iother_position = POPCNT(IAND(mask, 2**(iother - 1) - 1))
188 207 : iother_direction = MOD(code/3**iother_position, 3) + 1
189 207 : IF (iother_direction /= direction) CYCLE
190 138 : IF (iother < field) EXIT
191 614 : multiplicity = multiplicity + 1
192 : END DO
193 214 : IF (iother <= order) CYCLE
194 145 : previous_code = MOD(code, radix) + (code/(3*radix))*radix
195 145 : previous_mask = IBCLR(mask, field - 1)
196 145 : previous_start = offset_previous(previous_code, previous_mask)
197 145 : IF (order_subset == order) THEN
198 328 : DO ioutput = 1, 3
199 246 : weight = REAL(multiplicity, KIND=dp)/REAL(order, KIND=dp)
200 : curvature_term = curvature_density_trace(curvature, &
201 246 : & previous_density(:, :, previous_start), ioutput, direction)
202 : extra(ioutput, code + 1) = extra(ioutput, code + 1) - &
203 328 : curvature_term*weight
204 : END DO
205 : END IF
206 145 : previous_end = previous_start + ncoeff(previous_mask) - 1
207 : ! (ε_m-ε_n-Ω_S-iγ_mn)ρ_S = source from lower-order responses.
208 : CALL nlo_density_step_derivatives(order - order_subset, direction, &
209 : subset_omega, energies, rate1, rate2, &
210 : transition_width, &
211 : hamiltonian, connection, &
212 : previous_density(:, :, &
213 : previous_start:previous_end), &
214 145 : next_density(:, :, 1:ncoeff(mask)))
215 : density(:, :, current_start:current_start + ncoeff(mask) - 1) = &
216 : density(:, :, current_start:current_start + ncoeff(mask) - 1) + &
217 : next_density(:, :, 1:ncoeff(mask))*REAL(multiplicity, &
218 282548 : KIND=dp)/REAL(order_subset, KIND=dp)
219 : END DO
220 : END DO
221 : END DO
222 40 : CALL MOVE_ALLOC(density, previous_density)
223 53 : offset_previous = offset_current
224 : END DO
225 958 : DO index = 1, 3**order
226 3531 : IF (.NOT. ANY(selected(:, index))) CYCLE
227 87 : current_start = offset_previous(index - 1, full_mask)
228 361 : DO field = 1, 3
229 261 : IF (.NOT. selected(field, index)) CYCLE
230 : ! J_i^(n) = -Tr[v_i rho^(n)] + curvature contribution.
231 : current(field, index) = -SUM(velocity(:, :, field)* &
232 : TRANSPOSE(previous_density(:, :, current_start))) + &
233 158244 : extra(field, representative_code(index - 1, full_mask) + 1)
234 : END DO
235 : END DO
236 7586 : IF (.NOT. ALL(ieee_is_finite(REAL(current, KIND=dp))) .OR. &
237 : .NOT. ALL(ieee_is_finite(AIMAG(current)))) THEN
238 0 : CPABORT('NLO current response contains invalid data.')
239 : END IF
240 :
241 13 : CALL timestop(handle)
242 :
243 26 : END SUBROUTINE nlo_current_tensor_derivatives
244 :
245 : ! **************************************************************************************************
246 : !> \brief Store each required density once per multiset of (frequency, direction) pairs.
247 : !> \param order Number of incident optical fields.
248 : !> \param degree Number of fields in each subset represented at this stage.
249 : !> \param omega Photon energies of the labeled incident fields.
250 : !> \param required Mask of subset densities needed by the requested tensor elements.
251 : !> \param ncoeff Number of retained k-derivative coefficients for each field subset.
252 : !> \param offsets Starting positions of the compact density arrays.
253 : !> \param evaluate Marks one representative density per equal-frequency field multiset.
254 : !> \param representative_code Tensor-direction code of the representative density.
255 : !> \param total Total retained matrix-derivative coefficients.
256 : !> \note The normalized subset recursion is invariant under relabeling equal-frequency
257 : !> fields. Integer multiplicities identify exactly equal multisets, with no frequency
258 : !> tolerance. Offsets address compact matrix-derivative slices of length ncoeff(mask).
259 : ! **************************************************************************************************
260 40 : SUBROUTINE subset_density_offsets(order, degree, omega, required, ncoeff, offsets, evaluate, &
261 : representative_code, total)
262 : INTEGER, INTENT(IN) :: order, degree
263 : REAL(KIND=dp), INTENT(IN) :: omega(:)
264 : LOGICAL, INTENT(IN) :: required(0:242, 0:31)
265 : INTEGER, INTENT(IN) :: ncoeff(0:31)
266 : INTEGER, INTENT(OUT) :: offsets(0:242, 0:31)
267 : LOGICAL, INTENT(OUT) :: evaluate(0:242, 0:31)
268 : INTEGER, INTENT(OUT) :: representative_code(0:242, 0:31), total
269 :
270 : INTEGER :: codes(1024), counts(3, 5), groups(5), icode, idir, ifield, imask, iposition, &
271 : istate, jfield, nstates, signatures(3, 5, 1024), starts(1024)
272 :
273 40 : groups = 0
274 182 : DO ifield = 1, order
275 142 : groups(ifield) = ifield
276 200 : DO jfield = 1, ifield - 1
277 108 : IF (ABS(omega(ifield) - omega(jfield)) > 0.0_dp) CYCLE
278 90 : groups(ifield) = groups(jfield)
279 160 : EXIT
280 : END DO
281 : END DO
282 40 : offsets = 0
283 40 : evaluate = .FALSE.
284 40 : representative_code = 0
285 40 : total = 0
286 40 : nstates = 0
287 664 : DO imask = 1, 2**order - 1
288 624 : IF (POPCNT(imask) /= degree) CYCLE
289 3790 : DO icode = 0, 3**degree - 1
290 3603 : IF (.NOT. required(icode, imask)) CYCLE
291 275 : counts = 0
292 275 : iposition = 0
293 1211 : DO ifield = 1, order
294 936 : IF (.NOT. BTEST(imask, ifield - 1)) CYCLE
295 584 : idir = MOD(icode/3**iposition, 3) + 1
296 584 : counts(idir, groups(ifield)) = counts(idir, groups(ifield)) + 1
297 1211 : iposition = iposition + 1
298 : END DO
299 647 : DO istate = 1, nstates
300 4163 : IF (ALL(signatures(:, :, istate) == counts)) EXIT
301 : END DO
302 275 : IF (istate > nstates) THEN
303 106 : nstates = istate
304 2226 : signatures(:, :, istate) = counts
305 106 : starts(istate) = total + 1
306 106 : codes(istate) = icode
307 106 : total = total + ncoeff(imask)
308 106 : evaluate(icode, imask) = .TRUE.
309 : END IF
310 275 : offsets(icode, imask) = starts(istate)
311 4227 : representative_code(icode, imask) = codes(istate)
312 : END DO
313 : END DO
314 :
315 40 : END SUBROUTINE subset_density_offsets
316 :
317 : ! **************************************************************************************************
318 : !> \brief Trace a precomputed curvature matrix against one lower-order density.
319 : !> \param curvature Three independent Cartesian position-curvature matrices.
320 : !> \param density Lower-order density matrix for the selected field subset.
321 : !> \param ic Cartesian direction of the output current.
322 : !> \param ia Cartesian direction of the incident field.
323 : !> \return Tr[F_(ic,ia) density], with antisymmetry under ic and ia.
324 : ! **************************************************************************************************
325 246 : FUNCTION curvature_density_trace(curvature, density, ic, ia) RESULT(value)
326 : COMPLEX(dp), INTENT(IN) :: curvature(:, :, :), density(:, :)
327 : INTEGER, INTENT(IN) :: ic, ia
328 : COMPLEX(dp) :: value
329 :
330 : INTEGER :: lower, pair, upper
331 :
332 246 : value = 0.0_dp
333 246 : IF (ic == ia) RETURN
334 164 : lower = MIN(ic, ia)
335 164 : upper = MAX(ic, ia)
336 164 : pair = (lower - 1)*(6 - lower)/2 + upper - lower
337 96348 : value = SUM(curvature(:, :, pair)*TRANSPOSE(density))
338 164 : IF (ic > ia) value = -value
339 :
340 : END FUNCTION curvature_density_trace
341 :
342 : ! **************************************************************************************************
343 : !> \brief Advance one optical field in the normalized Cartesian k-derivative recursion.
344 : !> \param max_grade Required highest k-derivative grade, 0:4
345 : !> \param direction Cartesian component of the new input field, 1:3
346 : !> \param omega Cumulative input photon energy in Hartree
347 : !> \param energies Central band energies in Hartree
348 : !> \param rate1 Population damping energy in Hartree
349 : !> \param rate2 Coherence damping energy in Hartree
350 : !> \param transition_width Band-energy width of their smooth transition in Hartree
351 : !> \param hamiltonian Orthogonal-frame Hamiltonian k derivatives through the active response order
352 : !> \param connection Connection k derivatives through one grade below the response order
353 : !> \param previous Previous-order density normalized k derivatives, at least through grade
354 : !> max_grade+1
355 : !> \param response This order's density normalized k derivatives, at least through grade max_grade
356 : !> \note Write X_α=(∂_k^α X)/α! and L_H(X)=[H,X]. At order p, after adding
357 : !> the field component j with cumulative photon energy Ω, solve
358 : !> (L_H0-Ω-iΓ_0)ρ_{p,α} = -i(α_j+1)ρ_{p-1,α+e_j}
359 : !> -Σ_{β+γ=α}[A_{j,β},ρ_{p-1,γ}]
360 : !> -Σ_{0<β≤α}[H_β,ρ_{p,α-β}]
361 : !> +iΣ_{0<β≤α}Γ_β(ρ_{p,α-β}).
362 : !> For equal times Γ_0=1/T1_T2 and all higher Γ derivatives vanish.
363 : !> For separate times, Γ_H(X)=γ₂X+(γ₁−γ₂)exp[−(ad_H/w)²]X.
364 : !> The known positive-grade Γ_H terms enter the right-hand side.
365 : !> The central H_0 is diagonal in this frame.
366 : !> \note If previous has grade <=g+1, this produces the next response through grade g.
367 : !> Denominators use exp(-i ω t): ε_m-ε_n-ω-iγ_mn.
368 : !> The field couples as +E.r for an electron, r=i ∂_k+A.
369 : ! **************************************************************************************************
370 435 : SUBROUTINE nlo_density_step_derivatives(max_grade, direction, omega, energies, rate1, &
371 : rate2, transition_width, &
372 145 : hamiltonian, connection, previous, response)
373 : INTEGER, INTENT(IN) :: max_grade, direction
374 : REAL(KIND=dp), INTENT(IN) :: omega, energies(:), rate1, rate2, &
375 : transition_width
376 : COMPLEX(dp), INTENT(IN) :: hamiltonian(:, :, :), &
377 : connection(:, :, :, :), &
378 : previous(:, :, :)
379 : COMPLEX(dp), INTENT(OUT) :: response(:, :, :)
380 :
381 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_density_step_derivatives'
382 :
383 145 : COMPLEX(dp), ALLOCATABLE :: matrix_term(:, :), rhs(:, :), &
384 145 : transition_term(:, :, :)
385 : INTEGER :: grade, handle, i, ia, ib, j, k, n, &
386 : ngrade, nh, nlower, np, plus(3), &
387 : remain(3)
388 : INTEGER, DIMENSION(3, nlo_n_derivatives) :: alpha
389 : LOGICAL :: separate_times
390 : REAL(KIND=dp) :: pair_rate
391 145 : TYPE(nlo_band_transition_type) :: transition
392 :
393 145 : CALL timeset(routineN, handle)
394 :
395 278792 : response = 0.0_dp
396 145 : n = SIZE(energies)
397 145 : CPASSERT(max_grade >= 0 .AND. max_grade <= 4)
398 145 : CPASSERT(direction >= 1 .AND. direction <= 3)
399 145 : CPASSERT(n >= 2)
400 145 : nh = nlo_derivative_count(max_grade)
401 145 : np = nlo_derivative_count(max_grade + 1)
402 145 : CPASSERT(SIZE(hamiltonian, 1) == n)
403 145 : CPASSERT(SIZE(hamiltonian, 2) == n)
404 145 : CPASSERT(SIZE(hamiltonian, 3) >= nh)
405 145 : CPASSERT(SIZE(connection, 1) == n)
406 145 : CPASSERT(SIZE(connection, 2) == n)
407 145 : CPASSERT(SIZE(connection, 3) == 3)
408 145 : CPASSERT(SIZE(connection, 4) >= nh)
409 145 : CPASSERT(SIZE(previous, 1) == n)
410 145 : CPASSERT(SIZE(previous, 2) == n)
411 145 : CPASSERT(SIZE(previous, 3) >= np)
412 145 : CPASSERT(SIZE(response, 1) == n)
413 145 : CPASSERT(SIZE(response, 2) == n)
414 145 : CPASSERT(SIZE(response, 3) >= nh)
415 : IF (.NOT. ieee_is_finite(omega) .OR. .NOT. ALL(ieee_is_finite(energies)) .OR. &
416 : .NOT. ieee_is_finite(rate1) .OR. .NOT. ieee_is_finite(rate2) .OR. &
417 : .NOT. ieee_is_finite(transition_width) .OR. rate1 <= 0.0_dp .OR. &
418 : rate2 < rate1 .OR. transition_width <= 0.0_dp .OR. &
419 1329185 : .NOT. ALL(ieee_is_finite(REAL(previous, KIND=dp))) .OR. &
420 : .NOT. ALL(ieee_is_finite(AIMAG(previous)))) THEN
421 0 : CPABORT('NLO density response contains invalid data.')
422 : END IF
423 870 : ALLOCATE (matrix_term(n, n), rhs(n, n))
424 145 : separate_times = rate1 < rate2
425 145 : IF (separate_times) THEN
426 25 : ALLOCATE (transition_term(n, n, (max_grade + 1)*(max_grade + 2)/2))
427 5 : CALL nlo_band_transition_prepare(energies, transition_width, rate1, rate2, transition)
428 : END IF
429 145 : CALL nlo_derivative_indices(alpha)
430 145 : nlower = 0
431 714 : DO ia = 1, nlo_derivative_count(max_grade)
432 2276 : grade = SUM(alpha(:, ia))
433 569 : IF (separate_times .AND. grade > 0) THEN
434 65 : IF (ia == nlo_derivative_count(grade - 1) + 1) THEN
435 10 : nlower = ia - 1
436 10 : ngrade = nlo_derivative_count(grade) - nlower
437 : CALL nlo_band_transition_correction(transition, hamiltonian, response, grade, &
438 10 : transition_term(:, :, :ngrade))
439 : END IF
440 : END IF
441 2276 : plus = alpha(:, ia)
442 569 : plus(direction) = plus(direction) + 1
443 569 : ib = derivatives_index(alpha, plus)
444 569 : CPASSERT(ib /= 0)
445 : ! D_a rho = d_a rho - i[A_a,rho]; lower H derivatives enter recursively.
446 278647 : rhs(:, :) = -gaussi*REAL(plus(direction), KIND=dp)*previous(:, :, ib)
447 4526 : DO ib = 1, ia
448 15828 : remain = alpha(:, ia) - alpha(:, ib)
449 3957 : j = derivatives_index(alpha, remain)
450 3957 : IF (j == 0) CYCLE
451 20702103 : matrix_term(:, :) = MATMUL(connection(:, :, direction, ib), previous(:, :, j))
452 838383 : rhs(:, :) = rhs - matrix_term
453 20702103 : matrix_term(:, :) = MATMUL(previous(:, :, j), connection(:, :, direction, ib))
454 841196 : rhs(:, :) = rhs + matrix_term
455 : END DO
456 3957 : DO ib = 2, ia
457 13552 : remain = alpha(:, ia) - alpha(:, ib)
458 3388 : j = derivatives_index(alpha, remain)
459 3388 : IF (j == 0) CYCLE
460 13820792 : matrix_term(:, :) = MATMUL(hamiltonian(:, :, ib), response(:, :, j))
461 559736 : rhs(:, :) = rhs - matrix_term
462 13820792 : matrix_term(:, :) = MATMUL(response(:, :, j), hamiltonian(:, :, ib))
463 562549 : rhs(:, :) = rhs + matrix_term
464 : END DO
465 569 : IF (separate_times .AND. grade > 0) THEN
466 45695 : rhs(:, :) = rhs + gaussi*transition_term(:, :, ia - nlower)
467 : END IF
468 12228 : DO k = 1, n
469 278647 : DO i = 1, n
470 266564 : pair_rate = rate1
471 266564 : IF (separate_times) pair_rate = nlo_band_transition_rate(transition, i, k)
472 : response(i, k, ia) = rhs(i, k)/ &
473 278078 : CMPLX(energies(i) - energies(k) - omega, -pair_rate, dp)
474 : END DO
475 : END DO
476 : END DO
477 557584 : IF (.NOT. ALL(ieee_is_finite(REAL(response, KIND=dp))) .OR. &
478 : .NOT. ALL(ieee_is_finite(AIMAG(response)))) THEN
479 0 : CPABORT('NLO density response contains invalid data.')
480 : END IF
481 :
482 145 : CALL timestop(handle)
483 :
484 290 : END SUBROUTINE nlo_density_step_derivatives
485 :
486 : END MODULE nlo_current
|
