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 : !> \author Jan Wilhelm, with help from ChatGPT
9 :
10 : MODULE nlo_integration
11 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
12 : USE ieee_arithmetic, ONLY: ieee_is_finite
13 : USE kinds, ONLY: dp
14 : USE kubo_projected, ONLY: projected_eigensystem
15 : USE machine, ONLY: m_flush,&
16 : m_walltime
17 : USE mathconstants, ONLY: gaussi,&
18 : twopi
19 : USE nlo_current, ONLY: nlo_current_tensor_derivatives
20 : USE nlo_derivatives, ONLY: nlo_derivative_count,&
21 : nlo_derivative_indices,&
22 : nlo_equilibrium_density_derivatives,&
23 : nlo_n_derivatives,&
24 : nlo_orthonormal_operator_derivatives
25 : USE nlo_read_input, ONLY: nlo_frequency_case
26 : USE nlo_types, ONLY: nlo_type
27 : USE physcon, ONLY: evolt
28 : #include "../base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 :
34 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nlo_integration'
35 :
36 : PUBLIC :: nlo_bz_integration
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief Integrate the density-response current over the Brillouin zone.
42 : !> \param nlo_env NLO settings and data for the calculation
43 : ! **************************************************************************************************
44 22 : SUBROUTINE nlo_bz_integration(nlo_env)
45 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
46 :
47 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_bz_integration'
48 :
49 22 : COMPLEX(dp), ALLOCATABLE :: curvature(:, :, :), matrix_term(:, :)
50 : INTEGER :: handle, i_case, ia, ic, ik, iw, nocc, &
51 : pair, spin
52 : LOGICAL :: zero_output
53 : REAL(KIND=dp) :: energy_ev(5), k(3), local_projection, &
54 : occ_factor, omega(5)
55 :
56 22 : CALL timeset(routineN, handle)
57 :
58 22 : CALL allocate_nlo_quantities(nlo_env)
59 : ALLOCATE (curvature(nlo_env%n_bands, nlo_env%n_bands, 3), &
60 176 : matrix_term(nlo_env%n_bands, nlo_env%n_bands))
61 :
62 22 : nlo_env%projection_error = 0.0_dp
63 22 : nlo_env%min_direct_gap = HUGE(1.0_dp)
64 22 : nlo_env%progress_total = 1 + (nlo_env%n_kpoints - 1)/nlo_env%para_env%num_pe
65 22 : IF (nlo_env%progress_total > HUGE(nlo_env%progress_total)/nlo_env%n_spin) THEN
66 0 : CPABORT('NLO work count exceeds integer indexing.')
67 : END IF
68 22 : nlo_env%progress_total = nlo_env%n_spin*nlo_env%progress_total
69 22 : nlo_env%progress_done = 0
70 22 : nlo_env%progress_next = 1
71 22 : nlo_env%progress_start = m_walltime()
72 :
73 44 : DO spin = 1, nlo_env%n_spin
74 22 : IF (nlo_env%mos(spin)%maxocc <= 0.0_dp) THEN
75 0 : CPABORT('NONLINEAR_OPTICS requires positive band occupation.')
76 : END IF
77 22 : occ_factor = nlo_env%mos(spin)%maxocc
78 22 : nocc = NINT(REAL(nlo_env%mos(spin)%nelectron, KIND=dp)/occ_factor)
79 22 : IF (nocc < 1 .OR. nocc >= nlo_env%n_bands .OR. &
80 : ABS(nocc*occ_factor - nlo_env%mos(spin)%nelectron) > 1.e-8_dp) THEN
81 0 : CPABORT('NLO needs integer occupied-band count and empty bands.')
82 : END IF
83 498 : nlo_env%occupied = 0
84 108 : nlo_env%occupied(1:nocc) = 1
85 68 : DO ik = 1, nlo_env%n_kpoints
86 24 : IF (MOD(ik - 1, nlo_env%para_env%num_pe) /= nlo_env%para_env%mepos) CYCLE
87 : k = REAL([MOD(ik - 1, nlo_env%k_mesh(1)), &
88 : MOD((ik - 1)/nlo_env%k_mesh(1), nlo_env%k_mesh(2)), &
89 : (ik - 1)/(nlo_env%k_mesh(1)*nlo_env%k_mesh(2))], KIND=dp)/ &
90 48 : REAL(nlo_env%k_mesh, KIND=dp)
91 : ! Build H(k), S(k), d(k) and their analytic Cartesian k derivatives.
92 12 : CALL ao_derivatives_at_k(nlo_env, k, spin)
93 12 : CALL hermitize_derivatives(nlo_env%hamiltonian_derivatives, local_projection)
94 :
95 12 : nlo_env%projection_error = MAX(nlo_env%projection_error, local_projection)
96 12 : IF (local_projection > 1.e-5_dp) THEN
97 0 : CPABORT('NLO Hamiltonian images are not Hermitian.')
98 : END IF
99 : ! X = S^(-1/2) orthonormalizes the AO basis; U diagonalizes X H X here.
100 : ! Express all operator derivatives in this same eigenvector basis U.
101 : CALL nlo_ao_to_band_derivatives( &
102 : nlo_env%overlap_derivatives, nlo_env%hamiltonian_derivatives, &
103 : nlo_env%position_derivatives, &
104 : nlo_env%occupied, nlo_env%energies, nlo_env%band_hamiltonian_derivatives, &
105 : nlo_env%band_connection_derivatives, nlo_env%rho0_derivatives, &
106 : nlo_env%velocity, &
107 12 : active_grade=nlo_env%order)
108 :
109 : nlo_env%min_direct_gap = MIN( &
110 : nlo_env%min_direct_gap, &
111 12 : nlo_env%energies(nocc + 1) - nlo_env%energies(nocc))
112 : ! F_ca = d_c A_a - d_a A_c - i[A_c,A_a], reused for every frequency.
113 36 : DO ic = 1, 2
114 72 : DO ia = ic + 1, 3
115 36 : pair = (ic - 1)*(6 - ic)/2 + ia - ic
116 : curvature(:, :, pair) = &
117 : nlo_env%band_connection_derivatives(:, :, ia, ic + 1) - &
118 18204 : nlo_env%band_connection_derivatives(:, :, ic, ia + 1)
119 36 : matrix_term(:, :) = MATMUL(nlo_env%band_connection_derivatives(:, :, ic, 1), &
120 885852 : nlo_env%band_connection_derivatives(:, :, ia, 1))
121 18204 : curvature(:, :, pair) = curvature(:, :, pair) - gaussi*matrix_term
122 36 : matrix_term(:, :) = MATMUL(nlo_env%band_connection_derivatives(:, :, ia, 1), &
123 885852 : nlo_env%band_connection_derivatives(:, :, ic, 1))
124 18228 : curvature(:, :, pair) = curvature(:, :, pair) + gaussi*matrix_term
125 : END DO
126 : END DO
127 25 : DO i_case = 1, nlo_env%n_frequency_combinations
128 : CALL nlo_frequency_case(nlo_env%frequency, nlo_env%order, i_case, &
129 13 : energy_ev, zero_output)
130 13 : IF (zero_output) THEN
131 0 : CPABORT('Invalid NONLINEAR_OPTICS frequency case.')
132 : END IF
133 13 : omega = 0.0_dp
134 53 : omega(1:nlo_env%order) = energy_ev(1:nlo_env%order)/evolt
135 : ! Compute J_i^(n)(k) from rho^(n): -Tr[v_i rho^(n)] plus curvature current.
136 : ! Each rho^(n) follows from lower orders at Omega = sum_j omega_j.
137 : CALL nlo_current_tensor_derivatives( &
138 : nlo_env%order, omega(1:nlo_env%order), nlo_env%energies, &
139 : 1.0_dp/nlo_env%t1_au, 1.0_dp/nlo_env%t2_au, &
140 : nlo_env%transition_width_au, nlo_env%band_hamiltonian_derivatives, &
141 : nlo_env%band_connection_derivatives, curvature, &
142 : nlo_env%rho0_derivatives, &
143 : nlo_env%velocity, &
144 13 : nlo_env%current_flat, elements=nlo_env%elements)
145 :
146 : ! Sum J_i^(n)(k) over spin and k; f_spin is the band occupation factor.
147 : ! I_i^(n) = sum_(spin,k) f_spin J_i^(n)(k) / N_k.
148 : nlo_env%current_spectrum(:, :, i_case) = &
149 : nlo_env%current_spectrum(:, :, i_case) + &
150 3805 : occ_factor*nlo_env%current_flat/REAL(nlo_env%n_kpoints, KIND=dp)
151 : END DO
152 46 : IF (nlo_env%para_env%mepos == 0) CALL report_nlo_progress(nlo_env)
153 :
154 : END DO
155 : END DO
156 22 : CALL nlo_env%para_env%max(nlo_env%projection_error)
157 22 : CALL nlo_env%para_env%min(nlo_env%min_direct_gap)
158 22 : iw = cp_logger_get_default_io_unit()
159 22 : IF (iw > 0) THEN
160 : WRITE (iw, '(T2,A,T69,F9.2,1X,A)') &
161 11 : 'NONLINEAR_OPTICS| Minimum occupied-empty gap:', &
162 22 : nlo_env%min_direct_gap*evolt, 'eV'
163 : END IF
164 22 : CALL release_nlo_quantities(nlo_env)
165 :
166 22 : CALL timestop(handle)
167 :
168 44 : END SUBROUTINE nlo_bz_integration
169 :
170 : ! **************************************************************************************************
171 : !> \brief Project all analytic Hamiltonian derivatives onto the Hermitian image subspace.
172 : !> \param derivatives AO Hamiltonian derivatives over lattice images
173 : !> \param relative_error Largest relative change from enforcing Hermiticity
174 : !> \note This linear projection commutes with all k derivatives.
175 : ! **************************************************************************************************
176 12 : SUBROUTINE hermitize_derivatives(derivatives, relative_error)
177 : COMPLEX(dp), INTENT(INOUT) :: derivatives(:, :, :)
178 : REAL(KIND=dp), INTENT(OUT) :: relative_error
179 :
180 24 : COMPLEX(dp) :: adjoint(SIZE(derivatives, 1), SIZE(derivatives, 2))
181 : INTEGER :: ij
182 : REAL(KIND=dp) :: scale
183 :
184 12 : relative_error = 0.0_dp
185 325 : DO ij = 1, SIZE(derivatives, 3)
186 157287 : adjoint = CONJG(TRANSPOSE(derivatives(:, :, ij)))
187 157287 : scale = MAX(1.0_dp, MAXVAL(ABS(derivatives(:, :, ij))))
188 157287 : relative_error = MAX(relative_error, MAXVAL(ABS(derivatives(:, :, ij) - adjoint))/scale)
189 157299 : derivatives(:, :, ij) = 0.5_dp*(derivatives(:, :, ij) + adjoint)
190 : END DO
191 :
192 12 : END SUBROUTINE hermitize_derivatives
193 :
194 : ! **************************************************************************************************
195 : !> \brief Fourier-transform AO H, S, and position matrices at one optical k point.
196 : !> \param nlo_env NLO settings and data for the calculation
197 : !> \param k Optical k point in fractional coordinates
198 : !> \param spin Spin channel of the converged SCF solution
199 : !> \note For X = H, S, or the AO position matrix d_a, and Cartesian k:
200 : !> dX(k)/dk_x = sum_R i R_x exp(i k·R) X(R).
201 : !> R is a real-space lattice vector. The input k uses fractional coordinates:
202 : !> k·R = 2 pi k_frac·R_int. Higher derivatives add factors i R_a;
203 : !> each stored Taylor coefficient is divided by its factorial.
204 : ! **************************************************************************************************
205 12 : SUBROUTINE ao_derivatives_at_k(nlo_env, k, spin)
206 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
207 : REAL(KIND=dp), INTENT(IN) :: k(3)
208 : INTEGER, INTENT(IN) :: spin
209 :
210 : CHARACTER(len=*), PARAMETER :: routineN = 'ao_derivatives_at_k'
211 :
212 : INTEGER :: alpha(3, nlo_n_derivatives), grade, &
213 : handle, ia, ij, j, n, na, nh, ni
214 : REAL(KIND=dp) :: divisor
215 :
216 12 : CALL timeset(routineN, handle)
217 :
218 : ASSOCIATE (h_image => nlo_env%h_ao_R(spin, :, :, :), &
219 : s_image => nlo_env%s_ao_R(1, :, :, :), &
220 : moment_image => nlo_env%d_ao_R, &
221 : image_cells => nlo_env%images%index_to_cell, &
222 : hmat => nlo_env%cell%hmat, &
223 : h_derivatives => nlo_env%hamiltonian_derivatives, &
224 : s_derivatives => nlo_env%overlap_derivatives, &
225 : moment_derivatives => nlo_env%position_derivatives)
226 12 : h_derivatives = 0.0_dp
227 12 : s_derivatives = 0.0_dp
228 12 : moment_derivatives = 0.0_dp
229 12 : n = SIZE(h_image, 1)
230 12 : ni = SIZE(s_image, 3)
231 0 : CPASSERT(n >= 2)
232 12 : CPASSERT(SIZE(h_image, 2) == n)
233 12 : CPASSERT(SIZE(h_image, 3) >= 1)
234 12 : CPASSERT(SIZE(s_image, 1) == n)
235 12 : CPASSERT(SIZE(s_image, 2) == n)
236 12 : CPASSERT(ni >= 1)
237 60 : CPASSERT(ALL(SHAPE(moment_image) == [3, n, n, ni]))
238 12 : CPASSERT(SIZE(h_image, 3) == ni)
239 36 : CPASSERT(ALL(SHAPE(image_cells) == [3, ni]))
240 12 : nh = nlo_derivative_count(nlo_env%order)
241 12 : na = nlo_derivative_count(nlo_env%order - 1)
242 48 : CPASSERT(ALL(SHAPE(h_derivatives) == [n, n, nh]))
243 48 : CPASSERT(ALL(SHAPE(s_derivatives) == [n, n, nh]))
244 60 : CPASSERT(ALL(SHAPE(moment_derivatives) == [n, n, 3, na]))
245 : IF (.NOT. ALL(ieee_is_finite(h_image)) .OR. .NOT. ALL(ieee_is_finite(s_image)) .OR. &
246 : .NOT. ALL(ieee_is_finite(moment_image)) .OR. &
247 677160 : .NOT. ALL(ieee_is_finite(k)) .OR. .NOT. ALL(ieee_is_finite(hmat))) THEN
248 0 : CPABORT('NONLINEAR_OPTICS numerical derivative construction failed.')
249 : END IF
250 12 : grade = nlo_env%order
251 12 : CPASSERT(grade >= 2 .AND. grade <= 5)
252 12 : CALL nlo_derivative_indices(alpha)
253 325 : DO ij = 1, nlo_derivative_count(grade)
254 : divisor = 1.0_dp
255 1252 : DO ia = 1, 3
256 1672 : DO j = 2, alpha(ia, ij)
257 1359 : divisor = divisor*REAL(j, KIND=dp)
258 : END DO
259 : END DO
260 : CALL rs_to_kp_multi_deriv(h_image, h_derivatives(:, :, ij), image_cells, k, &
261 313 : alpha(:, ij), hmat)
262 : CALL rs_to_kp_multi_deriv(s_image, s_derivatives(:, :, ij), image_cells, k, &
263 313 : alpha(:, ij), hmat)
264 157287 : h_derivatives(:, :, ij) = h_derivatives(:, :, ij)/divisor
265 157287 : s_derivatives(:, :, ij) = s_derivatives(:, :, ij)/divisor
266 1252 : IF (SUM(alpha(:, ij)) >= grade) CYCLE
267 712 : DO ia = 1, 3
268 : CALL rs_to_kp_multi_deriv(moment_image(ia, :, :, :), &
269 : moment_derivatives(:, :, ia, ij), image_cells, k, &
270 525 : alpha(:, ij), hmat)
271 264604 : moment_derivatives(:, :, ia, ij) = moment_derivatives(:, :, ia, ij)/divisor
272 : END DO
273 : END DO
274 : IF (.NOT. ALL(ieee_is_finite(REAL(h_derivatives, KIND=dp))) &
275 : .OR. .NOT. ALL(ieee_is_finite(AIMAG(h_derivatives))) .OR. &
276 : .NOT. ALL(ieee_is_finite(REAL(s_derivatives, KIND=dp))) &
277 : .OR. .NOT. ALL(ieee_is_finite(AIMAG(s_derivatives))) .OR. &
278 1158152 : .NOT. ALL(ieee_is_finite(REAL(moment_derivatives, KIND=dp))) .OR. &
279 12 : .NOT. ALL(ieee_is_finite(AIMAG(moment_derivatives)))) THEN
280 0 : CPABORT('NONLINEAR_OPTICS numerical derivative construction failed.')
281 : END IF
282 :
283 : END ASSOCIATE
284 :
285 12 : CALL timestop(handle)
286 :
287 12 : END SUBROUTINE ao_derivatives_at_k
288 :
289 : ! **************************************************************************************************
290 : !> \brief Convert AO operators to band matrices and build the equilibrium density.
291 : !> \param overlap AO overlap and its k derivatives
292 : !> \param hamiltonian AO Hamiltonian and its k derivatives
293 : !> \param moment AO position-matrix derivatives through one grade below the response order
294 : !> \param occupied One for an occupied band, zero for an empty band
295 : !> \param energies Band energies in Hartree at this k point
296 : !> \param band_h Hamiltonian derivatives in the eigenvector basis
297 : !> \param band_a Position connection derivatives in the eigenvector basis
298 : !> \param rho0 Equilibrium density and its k derivatives for unit occupation
299 : !> \param velocity Band velocity at this k point, in bohr per atomic time
300 : !> \param active_grade Highest derivative grade required for the response
301 : !> \note First X=S^(-1/2) makes the AO basis orthonormal. U diagonalizes
302 : !> X H X at this k point. The same U is applied to every k derivative
303 : !> of H and the position matrix, so no derivative of U is needed.
304 : !> Occupied bands then define rho0; field-driven densities are built later.
305 : ! **************************************************************************************************
306 12 : SUBROUTINE nlo_ao_to_band_derivatives(overlap, hamiltonian, moment, occupied, &
307 12 : energies, band_h, band_a, rho0, velocity, &
308 : active_grade)
309 : COMPLEX(dp), DIMENSION(:, :, :), INTENT(IN) :: overlap, hamiltonian
310 : COMPLEX(dp), INTENT(IN) :: moment(:, :, :, :)
311 : INTEGER, INTENT(IN) :: occupied(:)
312 : REAL(KIND=dp), INTENT(OUT) :: energies(:)
313 : COMPLEX(dp), INTENT(OUT) :: band_h(:, :, :), band_a(:, :, :, :), &
314 : rho0(:, :, :), velocity(:, :, :)
315 : INTEGER, INTENT(IN), OPTIONAL :: active_grade
316 :
317 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_ao_to_band_derivatives'
318 :
319 12 : COMPLEX(dp), ALLOCATABLE :: identity(:, :), matrix_term(:, :), &
320 12 : orth_h(:, :, :), u(:, :)
321 12 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: connection
322 : INTEGER :: alpha(3, nlo_n_derivatives), grade, &
323 : handle, ia, ideriv, ij, info, n, na, nh
324 :
325 12 : CALL timeset(routineN, handle)
326 :
327 260 : energies = 0.0_dp
328 157299 : band_h = 0.0_dp
329 264478 : band_a = 0.0_dp
330 157299 : rho0 = 0.0_dp
331 18216 : velocity = 0.0_dp
332 12 : n = SIZE(energies)
333 12 : CPASSERT(.NOT. (n < 2 .OR. SIZE(occupied) /= n))
334 12 : grade = 5
335 12 : IF (PRESENT(active_grade)) grade = active_grade
336 12 : CPASSERT(grade >= 2 .AND. grade <= 5)
337 12 : nh = nlo_derivative_count(grade)
338 12 : na = nlo_derivative_count(grade - 1)
339 48 : CPASSERT(ALL(SHAPE(overlap) == [n, n, nh]))
340 48 : CPASSERT(ALL(SHAPE(hamiltonian) == [n, n, nh]))
341 60 : CPASSERT(ALL(SHAPE(moment) == [n, n, 3, na]))
342 48 : CPASSERT(ALL(SHAPE(band_h) == [n, n, nh]))
343 60 : CPASSERT(ALL(SHAPE(band_a) == [n, n, 3, na]))
344 48 : CPASSERT(ALL(SHAPE(rho0) == [n, n, nh]))
345 48 : CPASSERT(ALL(SHAPE(velocity) == [n, n, 3]))
346 : ALLOCATE (connection(n, n, 3, na), identity(n, n), matrix_term(n, n), &
347 204 : orth_h(n, n, nh), u(n, n))
348 :
349 : CALL nlo_orthonormal_operator_derivatives(overlap, hamiltonian, moment, orth_h, &
350 12 : connection, active_grade=grade)
351 :
352 12 : identity = 0.0_dp
353 260 : DO ij = 1, n
354 260 : identity(ij, ij) = 1.0_dp
355 : END DO
356 :
357 12 : CALL projected_eigensystem(orth_h(:, :, 1), identity, energies, u, info)
358 12 : IF (info /= 0) CPABORT('NONLINEAR_OPTICS band eigensolver failed.')
359 :
360 325 : DO ij = 1, nh
361 15134215 : matrix_term(:, :) = MATMUL(orth_h(:, :, ij), u)
362 7802725 : band_h(:, :, ij) = MATMUL(CONJG(TRANSPOSE(u)), matrix_term)
363 313 : IF (ij > na) CYCLE
364 712 : DO ia = 1, 3
365 25438083 : matrix_term(:, :) = MATMUL(connection(:, :, ia, ij), u)
366 13115266 : band_a(:, :, ia, ij) = MATMUL(CONJG(TRANSPOSE(u)), matrix_term)
367 : END DO
368 : END DO
369 :
370 : CALL nlo_equilibrium_density_derivatives(energies, occupied, band_h, rho0, &
371 12 : active_grade=grade)
372 :
373 12 : CALL nlo_derivative_indices(alpha)
374 48 : DO ia = 1, 3
375 : ideriv = 0
376 975 : DO ij = 1, nh
377 3792 : IF (alpha(ia, ij) == 1 .AND. SUM(alpha(:, ij)) == 1) ideriv = ij
378 : END DO
379 36 : CPASSERT(ideriv /= 0)
380 452028 : matrix_term(:, :) = MATMUL(band_a(:, :, ia, 1), band_h(:, :, 1))
381 18204 : velocity(:, :, ia) = band_h(:, :, ideriv) - gaussi*matrix_term
382 452028 : matrix_term(:, :) = MATMUL(band_h(:, :, 1), band_a(:, :, ia, 1))
383 18216 : velocity(:, :, ia) = velocity(:, :, ia) + gaussi*matrix_term
384 : END DO
385 36432 : IF (.NOT. ALL(ieee_is_finite(REAL(velocity, KIND=dp))) .OR. &
386 : .NOT. ALL(ieee_is_finite(AIMAG(velocity)))) THEN
387 0 : CPABORT('NONLINEAR_OPTICS numerical derivative construction failed.')
388 : END IF
389 :
390 12 : CALL timestop(handle)
391 :
392 24 : END SUBROUTINE nlo_ao_to_band_derivatives
393 :
394 : ! **************************************************************************************************
395 : !> \brief Analytic Cartesian mixed k derivative of a real-space AO matrix.
396 : !> \param rs_real Real-space matrix for every lattice image
397 : !> \param ks_complex Requested derivative in k space
398 : !> \param index_to_cell Integer lattice vector for each image
399 : !> \param xkp k point in reciprocal fractional coordinates
400 : !> \param deriv_orders Number of derivatives in Cartesian x, y, z
401 : !> \param hmat Direct-cell vectors as columns, in the same length units as k^-1
402 : !> \note The Fourier factor is exp(2πi k_frac dot R_int). Each Cartesian derivative
403 : !> contributes i R_cart. This avoids phase-sensitive finite differences.
404 : ! **************************************************************************************************
405 1151 : SUBROUTINE rs_to_kp_multi_deriv(rs_real, ks_complex, index_to_cell, xkp, deriv_orders, hmat)
406 : REAL(KIND=dp), INTENT(IN) :: rs_real(:, :, :)
407 : COMPLEX(dp), INTENT(OUT) :: ks_complex(:, :)
408 : INTEGER, INTENT(IN) :: index_to_cell(:, :)
409 : REAL(KIND=dp), INTENT(IN) :: xkp(3)
410 : INTEGER, INTENT(IN) :: deriv_orders(3)
411 : REAL(KIND=dp), INTENT(IN) :: hmat(3, 3)
412 :
413 : COMPLEX(dp) :: factor
414 : INTEGER :: i_axis, i_image, n_images
415 : REAL(KIND=dp) :: cell_index(3), cell_vector(3), phase
416 :
417 4604 : IF (ANY(deriv_orders < 0)) CPABORT('Negative k derivative order requested.')
418 1151 : n_images = SIZE(rs_real, 3)
419 1151 : IF (n_images < 1) CPABORT('No real-space lattice images supplied.')
420 1151 : IF (SIZE(index_to_cell, 1) /= 3 .OR. SIZE(index_to_cell, 2) /= n_images) THEN
421 0 : CPABORT('Real-space matrices and lattice-image indices are inconsistent.')
422 : END IF
423 1151 : IF (SIZE(ks_complex, 1) /= SIZE(rs_real, 1) .OR. SIZE(ks_complex, 2) /= SIZE(rs_real, 2)) THEN
424 0 : CPABORT('Real-space and k-space matrix shapes are inconsistent.')
425 : END IF
426 578865 : ks_complex = CMPLX(0.0_dp, 0.0_dp, kind=dp)
427 23020 : DO i_image = 1, n_images
428 87476 : cell_index = REAL(index_to_cell(:, i_image), KIND=dp)
429 284297 : cell_vector = MATMUL(hmat, cell_index)
430 87476 : phase = twopi*SUM(xkp*cell_index)
431 21869 : factor = CMPLX(COS(phase), SIN(phase), kind=dp)
432 87476 : DO i_axis = 1, 3
433 87476 : factor = factor*CMPLX(0.0_dp, cell_vector(i_axis), kind=dp)**deriv_orders(i_axis)
434 : END DO
435 10999586 : ks_complex = ks_complex + factor*rs_real(:, :, i_image)
436 : END DO
437 :
438 1151 : END SUBROUTINE rs_to_kp_multi_deriv
439 :
440 : ! **************************************************************************************************
441 : !> \brief Allocate the per-k matrices and rank-local frequency accumulators once.
442 : !> \param nlo_env NLO settings and data for the calculation
443 : ! **************************************************************************************************
444 22 : SUBROUTINE allocate_nlo_quantities(nlo_env)
445 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
446 :
447 : INTEGER :: na, nh
448 :
449 22 : nh = nlo_derivative_count(nlo_env%order)
450 22 : na = nlo_derivative_count(nlo_env%order - 1)
451 66 : ALLOCATE (nlo_env%energies(nlo_env%n_bands))
452 66 : ALLOCATE (nlo_env%occupied(nlo_env%n_bands))
453 110 : ALLOCATE (nlo_env%hamiltonian_derivatives(nlo_env%n_bands, nlo_env%n_bands, nh))
454 88 : ALLOCATE (nlo_env%overlap_derivatives(nlo_env%n_bands, nlo_env%n_bands, nh))
455 132 : ALLOCATE (nlo_env%position_derivatives(nlo_env%n_bands, nlo_env%n_bands, 3, na))
456 88 : ALLOCATE (nlo_env%band_hamiltonian_derivatives(nlo_env%n_bands, nlo_env%n_bands, nh))
457 88 : ALLOCATE (nlo_env%band_connection_derivatives(nlo_env%n_bands, nlo_env%n_bands, 3, na))
458 88 : ALLOCATE (nlo_env%rho0_derivatives(nlo_env%n_bands, nlo_env%n_bands, nh))
459 88 : ALLOCATE (nlo_env%velocity(nlo_env%n_bands, nlo_env%n_bands, 3))
460 66 : ALLOCATE (nlo_env%current_flat(3, 3**nlo_env%order))
461 0 : ALLOCATE (nlo_env%current_spectrum(3, 3**nlo_env%order, nlo_env%n_frequency_combinations), &
462 7456 : SOURCE=(0.0_dp, 0.0_dp))
463 :
464 22 : END SUBROUTINE allocate_nlo_quantities
465 :
466 : ! **************************************************************************************************
467 : !> \brief Print completed spin/k tasks when the next percentage threshold is crossed.
468 : !> \param nlo_env NLO settings and data for the calculation
469 : ! **************************************************************************************************
470 11 : SUBROUTINE report_nlo_progress(nlo_env)
471 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
472 :
473 : REAL(KIND=dp), PARAMETER :: &
474 : marks(9) = [0.1_dp, 0.3_dp, 1.0_dp, 3.0_dp, 10.0_dp, 30.0_dp, 60.0_dp, 90.0_dp, 100.0_dp]
475 :
476 : INTEGER :: iw
477 : REAL(KIND=dp) :: percent
478 :
479 11 : nlo_env%progress_done = nlo_env%progress_done + 1
480 11 : percent = 100.0_dp*REAL(nlo_env%progress_done, KIND=dp)/REAL(nlo_env%progress_total, KIND=dp)
481 11 : CPASSERT(.NOT. (nlo_env%progress_next > SIZE(marks)))
482 11 : IF (percent < marks(nlo_env%progress_next)) THEN
483 0 : RETURN
484 : END IF
485 11 : iw = cp_logger_get_default_io_unit()
486 11 : IF (iw > 0) THEN
487 11 : IF (nlo_env%progress_total < 10000) THEN
488 : WRITE (iw, '(T2,A,I0,A,I0,7X,F5.1,A,4X,A,T69,F10.2,1X,A)') &
489 11 : 'NONLINEAR_OPTICS| Progress: ', nlo_env%progress_done, ' / ', &
490 11 : nlo_env%progress_total, percent, '%', 'Total time:', &
491 22 : m_walltime() - nlo_env%progress_start, 's'
492 : ELSE
493 : WRITE (iw, '(T2,A,I0,A,I0,7X,F5.1,A)') &
494 0 : 'NONLINEAR_OPTICS| Progress: ', nlo_env%progress_done, ' / ', &
495 0 : nlo_env%progress_total, percent, '%'
496 : WRITE (iw, '(T2,A,T69,F10.2,1X,A)') &
497 0 : 'NONLINEAR_OPTICS| Total time:', m_walltime() - nlo_env%progress_start, 's'
498 : END IF
499 11 : CALL m_flush(iw)
500 : END IF
501 110 : DO WHILE (nlo_env%progress_next <= SIZE(marks))
502 99 : IF (percent < marks(nlo_env%progress_next)) EXIT
503 99 : nlo_env%progress_next = nlo_env%progress_next + 1
504 : END DO
505 :
506 : END SUBROUTINE report_nlo_progress
507 :
508 : ! **************************************************************************************************
509 : !> \brief Free the per-k workspace after the MPI reduction and validation.
510 : !> \param nlo_env NLO settings and data for the calculation
511 : ! **************************************************************************************************
512 22 : SUBROUTINE release_nlo_quantities(nlo_env)
513 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
514 :
515 22 : DEALLOCATE (nlo_env%energies, nlo_env%occupied)
516 22 : DEALLOCATE (nlo_env%hamiltonian_derivatives, nlo_env%overlap_derivatives)
517 22 : DEALLOCATE (nlo_env%position_derivatives)
518 22 : DEALLOCATE (nlo_env%band_hamiltonian_derivatives, nlo_env%band_connection_derivatives)
519 22 : DEALLOCATE (nlo_env%rho0_derivatives, nlo_env%velocity, nlo_env%current_flat)
520 :
521 22 : END SUBROUTINE release_nlo_quantities
522 :
523 : END MODULE nlo_integration
|
