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
11 : USE cp_control_types, ONLY: dft_control_type
12 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
13 : dbcsr_p_type
14 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
15 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
16 : USE input_constants, ONLY: do_method_gapw,&
17 : do_method_gpw
18 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
19 : section_vals_type,&
20 : section_vals_val_get
21 : USE kinds, ONLY: dp
22 : USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices
23 : USE kpoint_methods, ONLY: kpoint_init_cell_index
24 : USE kpoint_types, ONLY: kpoint_create,&
25 : kpoint_type
26 : USE mathlib, ONLY: vector_product
27 : USE nlo_integration, ONLY: nlo_bz_integration
28 : USE nlo_read_input, ONLY: create_nonlinear_optics_section,&
29 : nlo_frequency_case,&
30 : read_nlo_input
31 : USE nlo_types, ONLY: nlo_type,&
32 : release_nlo_environment
33 : USE physcon, ONLY: a_bohr,&
34 : c_light,&
35 : e_charge,&
36 : evolt,&
37 : femtoseconds,&
38 : joule,&
39 : mu_perm
40 : USE qs_environment_types, ONLY: get_qs_env,&
41 : qs_environment_type
42 : USE qs_moments, ONLY: build_local_moment_matrix_rs_img
43 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
44 : USE scf_control_types, ONLY: scf_control_type
45 : #include "../base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nlo'
52 :
53 : PUBLIC :: create_nonlinear_optics_section, run_nlo
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Calculate frequency-dependent χ^(2) to χ^(5) after the ground-state SCF.
59 : !> \param qs_env Converged k-point SCF environment
60 : !> \note The independent-particle approximation uses the converged SCF Hamiltonian.
61 : !> Excitonic effects and time-dependent screening (TDDFT, GW-BSE) are absent.
62 : !>
63 : !> Method: length-gauge density-matrix response in the spirit of
64 : !> Aversa and Sipe, Phys. Rev. B 52, 14636 (1995),
65 : !> DOI: 10.1103/PhysRevB.52.14636; analytic covariant k derivatives follow
66 : !> Ventura et al., Phys. Rev. B 96, 035431 (2017),
67 : !> DOI: 10.1103/PhysRevB.96.035431. The finite AO projection, analytic
68 : !> position curvature, and numerical implementation are defined below
69 : !> and in nlo_current_tensor_derivatives.
70 : !>
71 : !> 1. Read the optical k mesh, incident frequencies, tensor elements,
72 : !> and either a common time T₁ = T₂ or separate relaxation times.
73 : !>
74 : !> 2. Build lattice-image matrices for the AO Hamiltonian H, overlap S,
75 : !> dipole d.
76 : !> d_μν(R) = ⟨ϕ_(μ0)|r|ϕ_(νR)⟩ is a position matrix; the electron
77 : !> charge enters the later SI conversion. At each optical k point:
78 : !>
79 : !> X(k) = Σ_R exp(2π i k_frac · R) X(R), with X = H, S, d.
80 : !>
81 : !> Differentiating the phase gives every required k derivative.
82 : !> The AO overlap is removed, then H and position matrices are
83 : !> expressed in the eigenvector basis of H at the current k point.
84 : !>
85 : !> 3. At each optical k point, diagonalize the Hamiltonian once.
86 : !> The occupied bands define the equilibrium density matrix.
87 : !> A light field changes that density matrix: first order uses the
88 : !> equilibrium matrix, second order uses first-order responses,
89 : !> and so on through fifth order. Each new response is obtained
90 : !> from the previous one by an analytic k derivative and a
91 : !> band-energy denominator ε_m - ε_n - Ω - iΓ_mn. Here ε_m and ε_n
92 : !> are band energies, Ω is the sum of the incident photon energies,
93 : !> and Γ is the band-pair relaxation rate. For T₁ = T₂ it is the
94 : !> constant 1/T₁. For separate times, with γᵢ = 1/Tᵢ and
95 : !> Δε = ε_m - ε_n, it is
96 : !>
97 : !> Γ_mn = γ₁ + (γ₂ - γ₁)[1 - exp(-(Δε/w)²)].
98 : !>
99 : !> The positive width w sets the smooth transition at nearly
100 : !> degenerate bands. Its basis-covariant operator form is
101 : !> Γ_H(X) = γ₂X + (γ₁ - γ₂)exp[-(ad_H/w)²]X,
102 : !> where ad_H(X) = [H,X]. Its analytic k derivatives enter the
103 : !> response. All possible orders of
104 : !> the incident fields are averaged. No numerical derivatives of
105 : !> separately diagonalized bands or their phases are used.
106 : !>
107 : !> The position operator in the band representation is
108 : !>
109 : !> r_a = i ∂/∂k_a + A_a .
110 : !>
111 : !> The subscript a denotes x, y or z. A_a is the matrix obtained
112 : !> from the AO dipole matrix after orthonormalizing and rotating
113 : !> into the central band basis. It is the finite-basis part of
114 : !> the position operator. When a field changes the current,
115 : !> an additional term from this matrix is required:
116 : !>
117 : !> F_ia = ∂A_a/∂k_i - ∂A_i/∂k_a - i [A_i, A_a] .
118 : !>
119 : !> Here i is the output-current direction. F_ia is the position
120 : !> curvature; its contribution to the current is always included.
121 : !> The other current contribution is the band-velocity matrix
122 : !> multiplied by the density response.
123 : !> Summing over optical k points and spin gives the
124 : !> cell-integrated current response I^(n) in atomic units,
125 : !> stored locally on each MPI rank.
126 : !>
127 : !> 4. Sum I^(n) across MPI ranks, convert it to SI units and
128 : !> divide by the simulation-cell volume V. The result is the
129 : !> nonlinear conductivity σ^(n) for each field order. For an
130 : !> output frequency Ω, the corresponding polarization obeys
131 : !> J(Ω) = -i Ω P(Ω). Hence:
132 : !>
133 : !> χ²_iab = σ²_iab / [-i (ω₁+ω₂) ε₀]
134 : !> χ³_iabc = σ³_iabc / [-i (ω₁+ω₂+ω₃) ε₀]
135 : !> χ⁴_iabcd = σ⁴_iabcd / [-i (ω₁+ω₂+ω₃+ω₄) ε₀]
136 : !> χ⁵_iabcde = σ⁵_iabcde / [-i (ω₁+ω₂+ω₃+ω₄+ω₅) ε₀] .
137 : !>
138 : !> The SI conversion is in nlo_write_susceptibility. For a
139 : !> two-dimensional slab of in-plane area A, the sheet response
140 : !> is χ_sheet = (V/A) χ_volume. Ω must be nonzero.
141 : !>
142 : !> 5. Release the NLO matrices.
143 : ! **************************************************************************************************
144 14227 : SUBROUTINE run_nlo(qs_env)
145 : TYPE(qs_environment_type), POINTER :: qs_env
146 :
147 : CHARACTER(len=*), PARAMETER :: routineN = 'run_nlo'
148 :
149 : INTEGER :: handle
150 : LOGICAL :: nlo_present
151 128043 : TYPE(nlo_type) :: nlo_env
152 : TYPE(section_vals_type), POINTER :: input, nlo_section
153 :
154 14227 : CALL get_qs_env(qs_env, input=input)
155 14227 : nlo_section => section_vals_get_subs_vals(input, 'PROPERTIES%NONLINEAR_OPTICS')
156 14227 : CALL section_vals_val_get(nlo_section, '_SECTION_PARAMETERS_', l_val=nlo_present)
157 14227 : IF (.NOT. nlo_present) RETURN
158 :
159 22 : CALL timeset(routineN, handle)
160 :
161 : ! 1. Read and validate the nonlinear optical input.
162 22 : CALL read_nlo_input(qs_env, nlo_env)
163 :
164 : ! 2. Prepare the AO matrices from the converged SCF calculation.
165 22 : CALL nlo_prepare_environment(qs_env, nlo_env)
166 :
167 : ! 3. Calculate the density response and integrate the current over k.
168 : ! I^(n) = -sum_(spin,k) w_k Tr[v(E) rho^(n)] in atomic units.
169 22 : CALL nlo_bz_integration(nlo_env)
170 :
171 : ! 4. Convert the current to susceptibility and print the requested tensors.
172 : ! chi^(n) = I^(n) / [V epsilon_0 (-i Omega)] after SI conversion.
173 22 : CALL nlo_write_susceptibility(nlo_env)
174 :
175 : ! 5. Release the NLO matrices.
176 22 : CALL release_nlo_environment(nlo_env)
177 :
178 22 : CALL timestop(handle)
179 :
180 14227 : END SUBROUTINE run_nlo
181 :
182 : ! **************************************************************************************************
183 : !> \brief Build the real-space Hamiltonian, overlap and position matrices from the SCF result.
184 : !> \param qs_env Converged k-point SCF environment
185 : !> \param nlo_env NLO settings and data for the calculation
186 : ! **************************************************************************************************
187 22 : SUBROUTINE nlo_prepare_environment(qs_env, nlo_env)
188 : TYPE(qs_environment_type), POINTER :: qs_env
189 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
190 :
191 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_prepare_environment'
192 :
193 : INTEGER :: handle, ik, iw
194 : LOGICAL :: scf_has_kpoints
195 22 : TYPE(dbcsr_p_type), POINTER :: dr(:, :), hr(:, :), sr(:, :)
196 : TYPE(dft_control_type), POINTER :: dft_control
197 : TYPE(kpoint_type), POINTER :: scf
198 22 : TYPE(neighbor_list_set_p_type), POINTER :: all_pairs(:)
199 : TYPE(scf_control_type), POINTER :: scf_control
200 :
201 22 : CALL timeset(routineN, handle)
202 :
203 : CALL get_qs_env(qs_env, do_kpoints=scf_has_kpoints, dft_control=dft_control, &
204 22 : scf_control=scf_control)
205 22 : IF (scf_control%smear%do_smear) CPABORT('NONLINEAR_OPTICS requires SCF without smearing.')
206 22 : IF (.NOT. scf_has_kpoints) THEN
207 0 : CPABORT('NONLINEAR_OPTICS requires k-point SCF; Gamma-only SCF is unsupported.')
208 : END IF
209 22 : SELECT CASE (dft_control%qs_control%method_id)
210 : CASE (do_method_gpw, do_method_gapw)
211 : CASE DEFAULT
212 22 : CPABORT('NONLINEAR_OPTICS currently supports GPW/GAPW only.')
213 : END SELECT
214 22 : nlo_env%n_kpoints = 1
215 88 : DO ik = 1, 3
216 66 : IF (nlo_env%n_kpoints > HUGE(nlo_env%n_kpoints)/nlo_env%k_mesh(ik)) THEN
217 0 : CPABORT('NONLINEAR_OPTICS K_MESH exceeds integer indexing.')
218 : END IF
219 88 : nlo_env%n_kpoints = nlo_env%n_kpoints*nlo_env%k_mesh(ik)
220 : END DO
221 22 : NULLIFY (nlo_env%images, dr)
222 : CALL get_qs_env(qs_env, cell=nlo_env%cell, para_env=nlo_env%para_env, matrix_ks_kp=hr, &
223 22 : matrix_s_kp=sr, sab_all=all_pairs, kpoints=scf, mos=nlo_env%mos)
224 88 : IF (ANY(nlo_env%cell%perd == 0 .AND. nlo_env%k_mesh /= 1)) THEN
225 0 : CPABORT('NONLINEAR_OPTICS K_MESH must be one in each nonperiodic direction.')
226 : END IF
227 22 : IF (nlo_env%cell%deth <= 0.0_dp) CPABORT('NONLINEAR_OPTICS requires positive cell volume.')
228 22 : CALL dbcsr_get_info(sr(1, 1)%matrix, nfullrows_total=nlo_env%n_ao)
229 22 : nlo_env%n_spin = SIZE(hr, 1)
230 22 : nlo_env%n_bands = nlo_env%n_ao
231 22 : CALL kpoint_create(nlo_env%images)
232 22 : CALL kpoint_init_cell_index(nlo_env%images, all_pairs, nlo_env%para_env, nlo_env%n_images)
233 1738 : ALLOCATE (dr(3, nlo_env%n_images))
234 22 : CALL build_local_moment_matrix_rs_img(qs_env, dr)
235 0 : ALLOCATE (nlo_env%h_ao_R(nlo_env%n_spin, nlo_env%n_ao, nlo_env%n_ao, nlo_env%n_images), &
236 443402 : SOURCE=0.0_dp)
237 443380 : ALLOCATE (nlo_env%s_ao_R(1, nlo_env%n_ao, nlo_env%n_ao, nlo_env%n_images), SOURCE=0.0_dp)
238 877188 : ALLOCATE (nlo_env%d_ao_R(3, nlo_env%n_ao, nlo_env%n_ao, nlo_env%n_images), SOURCE=0.0_dp)
239 : ! H(R), S(R), d(R) supply X(k) = sum_R exp(2 pi i k_frac.R) X(R).
240 22 : CALL replicate_rs_matrices(hr, scf, nlo_env%h_ao_R, nlo_env%images%cell_to_index)
241 22 : CALL replicate_rs_matrices(sr, scf, nlo_env%s_ao_R, nlo_env%images%cell_to_index)
242 22 : CALL replicate_rs_matrices(dr, nlo_env%images, nlo_env%d_ao_R, nlo_env%images%cell_to_index)
243 22 : CALL dbcsr_deallocate_matrix_set(dr)
244 22 : iw = cp_logger_get_default_io_unit()
245 22 : IF (iw > 0) THEN
246 11 : WRITE (iw, '()')
247 11 : WRITE (iw, '()')
248 11 : WRITE (iw, '(T2,A,T69,I12)') 'NONLINEAR_OPTICS| Selected tensor elements:', &
249 3658 : COUNT(nlo_env%elements)
250 11 : WRITE (iw, '(T2,A,T69,I12)') 'NONLINEAR_OPTICS| Bands in response:', nlo_env%n_bands
251 11 : IF (nlo_env%t1_au == nlo_env%t2_au) THEN
252 10 : WRITE (iw, '(T2,A,T70,F8.1,1X,A)') 'NONLINEAR_OPTICS| T1 = T2:', &
253 20 : nlo_env%t1_au*femtoseconds, 'fs'
254 10 : WRITE (iw, '(T2,A,T73,F8.1,1X,A)') 'NONLINEAR_OPTICS| ℏ/T1 = ℏ/T2:', &
255 20 : 1000.0_dp*evolt/nlo_env%t1_au, 'meV'
256 : ELSE
257 1 : WRITE (iw, '(T2,A,T70,F8.1,1X,A)') 'NONLINEAR_OPTICS| T1:', &
258 2 : nlo_env%t1_au*femtoseconds, 'fs'
259 1 : WRITE (iw, '(T2,A,T70,F8.1,1X,A)') 'NONLINEAR_OPTICS| T2:', &
260 2 : nlo_env%t2_au*femtoseconds, 'fs'
261 1 : WRITE (iw, '(T2,A,T71,F8.1,1X,A)') 'NONLINEAR_OPTICS| ℏ/T1:', &
262 2 : 1000.0_dp*evolt/nlo_env%t1_au, 'meV'
263 1 : WRITE (iw, '(T2,A,T71,F8.1,1X,A)') 'NONLINEAR_OPTICS| ℏ/T2:', &
264 2 : 1000.0_dp*evolt/nlo_env%t2_au, 'meV'
265 1 : WRITE (iw, '(T2,A,T69,F8.1,1X,A)') 'NONLINEAR_OPTICS| Band transition width:', &
266 2 : 1000.0_dp*nlo_env%transition_width_au*evolt, 'meV'
267 : END IF
268 : END IF
269 :
270 22 : CALL timestop(handle)
271 :
272 22 : END SUBROUTINE nlo_prepare_environment
273 :
274 : ! **************************************************************************************************
275 : !> \brief Sum the rank-local currents and convert them to volume and sheet susceptibilities.
276 : !> \param nlo_env NLO settings and data for the calculation
277 : ! **************************************************************************************************
278 22 : SUBROUTINE nlo_write_susceptibility(nlo_env)
279 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
280 :
281 : CHARACTER(len=*), PARAMETER :: routineN = 'nlo_write_susceptibility'
282 : CHARACTER(len=3), PARAMETER :: &
283 : omega_subscript(5) = [CHARACTER(len=3) :: '₁', '₂', '₃', '₄', '₅']
284 :
285 : CHARACTER(len=80) :: heading_2d, heading_3d
286 22 : COMPLEX(dp), ALLOCATABLE :: chi_flat(:, :), current_flat_sum(:, :)
287 : INTEGER :: a, c, code, handle, i_case, iw
288 : LOGICAL :: zero_output
289 : REAL(KIND=dp) :: area, energy_ev(5), factor, omega(5), &
290 : sheet_height
291 :
292 22 : CALL timeset(routineN, handle)
293 :
294 88 : ALLOCATE (current_flat_sum(3, 3**nlo_env%order), chi_flat(3, 3**nlo_env%order))
295 10 : SELECT CASE (nlo_env%order)
296 : CASE (2)
297 10 : heading_3d = 'NONLINEAR_OPTICS| χ² 3D (m/V)'
298 10 : heading_2d = 'NONLINEAR_OPTICS| χ² 2D (m²/V)'
299 : CASE (3)
300 4 : heading_3d = 'NONLINEAR_OPTICS| χ³ 3D (m²/V²)'
301 4 : heading_2d = 'NONLINEAR_OPTICS| χ³ 2D (m³/V²)'
302 : CASE (4)
303 2 : heading_3d = 'NONLINEAR_OPTICS| χ⁴ 3D (m³/V³)'
304 2 : heading_2d = 'NONLINEAR_OPTICS| χ⁴ 2D (m⁴/V³)'
305 : CASE (5)
306 6 : heading_3d = 'NONLINEAR_OPTICS| χ⁵ 3D (m⁴/V⁴)'
307 22 : heading_2d = 'NONLINEAR_OPTICS| χ⁵ 2D (m⁵/V⁴)'
308 : END SELECT
309 46 : DO i_case = 1, nlo_env%n_frequency_combinations
310 : CALL nlo_frequency_case(nlo_env%frequency, nlo_env%order, i_case, &
311 24 : energy_ev, zero_output)
312 24 : IF (zero_output) CPABORT('Invalid NONLINEAR_OPTICS frequency case.')
313 24 : omega = 0.0_dp
314 98 : omega(1:nlo_env%order) = energy_ev(1:nlo_env%order)/evolt
315 7368 : current_flat_sum(:, :) = nlo_env%current_spectrum(:, :, i_case)
316 24 : CALL nlo_env%para_env%sum(current_flat_sum)
317 : factor = e_charge**(nlo_env%order + 1)*a_bohr**(nlo_env%order - 2)* &
318 24 : mu_perm*c_light**2/(joule**nlo_env%order*nlo_env%cell%deth)
319 7488 : chi_flat(:, :) = current_flat_sum*factor/CMPLX(0.0_dp, -SUM(omega), dp)
320 : area = SQRT(SUM(vector_product(nlo_env%cell%hmat(:, 1), &
321 96 : nlo_env%cell%hmat(:, 2))**2))
322 24 : IF (area <= 0.0_dp) CPABORT('NONLINEAR_OPTICS requires nonzero in-plane area.')
323 24 : sheet_height = nlo_env%cell%deth/area*a_bohr
324 24 : iw = cp_logger_get_default_io_unit()
325 24 : IF (iw <= 0) CYCLE
326 12 : WRITE (iw, '()')
327 12 : WRITE (iw, '(T2,A,I0,A,I0)') 'NONLINEAR_OPTICS| Frequency case ', &
328 24 : i_case, ' / ', nlo_env%n_frequency_combinations
329 : ! Adjust tab columns for the multibyte ω and subscript characters.
330 49 : DO a = 1, nlo_env%order
331 37 : WRITE (iw, '(T2,A,A,A,T72,F9.3,1X,A)') 'NONLINEAR_OPTICS| ω', &
332 86 : omega_subscript(a), ' =', energy_ev(a), 'eV'
333 : END DO
334 12 : WRITE (iw, '(T2,A,T70,F9.3,1X,A)') 'NONLINEAR_OPTICS| ω_out =', &
335 61 : SUM(energy_ev(1:nlo_env%order)), 'eV'
336 12 : WRITE (iw, '(T2,A,T50,A,T72,A)') TRIM(heading_3d), 'Real', 'Imag'
337 930 : DO code = 1, 3**nlo_env%order
338 3684 : DO c = 1, 3
339 2754 : IF (.NOT. nlo_env%elements(c, code)) CYCLE
340 248 : WRITE (iw, '(T2,A,A6,T37,2ES22.12)') 'NONLINEAR_OPTICS| 3D ', &
341 248 : nlo_tensor_indices(nlo_env%order, c, code), &
342 3920 : REAL(chi_flat(c, code), KIND=dp), AIMAG(chi_flat(c, code))
343 : END DO
344 : END DO
345 58 : IF (ALL(nlo_env%cell%perd(1:2) /= 0) .AND. nlo_env%cell%perd(3) == 0) THEN
346 : WRITE (iw, '(T2,A,T62,ES15.6,1X,A)') &
347 0 : 'NONLINEAR_OPTICS| Effective sheet height:', sheet_height, 'm'
348 0 : WRITE (iw, '(T2,A,T50,A,T72,A)') TRIM(heading_2d), 'Real', 'Imag'
349 0 : DO code = 1, 3**nlo_env%order
350 0 : DO c = 1, 3
351 0 : IF (.NOT. nlo_env%elements(c, code)) CYCLE
352 0 : WRITE (iw, '(T2,A,A6,T37,2ES22.12)') 'NONLINEAR_OPTICS| 2D ', &
353 0 : nlo_tensor_indices(nlo_env%order, c, code), &
354 0 : REAL(chi_flat(c, code)*sheet_height, KIND=dp), &
355 0 : AIMAG(chi_flat(c, code)*sheet_height)
356 : END DO
357 : END DO
358 : END IF
359 : END DO
360 :
361 22 : CALL timestop(handle)
362 :
363 44 : END SUBROUTINE nlo_write_susceptibility
364 :
365 : ! **************************************************************************************************
366 : !> \brief Decode the output and input Cartesian directions of one tensor element.
367 : !> \param order Rank of the susceptibility tensor
368 : !> \param output_direction Cartesian direction of the induced current
369 : !> \param code Base-three code for the field directions
370 : !> \return Cartesian indices for one tensor component
371 : ! **************************************************************************************************
372 248 : PURE FUNCTION nlo_tensor_indices(order, output_direction, code) RESULT(indices)
373 : INTEGER, INTENT(IN) :: order, output_direction, code
374 : CHARACTER(len=6) :: indices
375 :
376 : CHARACTER(len=*), PARAMETER :: axes = 'XYZ'
377 :
378 : INTEGER :: field, input_direction, remaining_code
379 :
380 248 : indices = ' '
381 248 : indices(1:1) = axes(output_direction:output_direction)
382 248 : remaining_code = code - 1
383 837 : DO field = 1, order
384 589 : input_direction = MOD(remaining_code, 3) + 1
385 589 : indices(field + 1:field + 1) = axes(input_direction:input_direction)
386 837 : remaining_code = remaining_code/3
387 : END DO
388 :
389 248 : END FUNCTION nlo_tensor_indices
390 :
391 : END MODULE nlo
|
