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_read_input
11 : !&<
12 : USE bibliography, ONLY: Cheng2015
13 : USE cp_units, ONLY: cp_unit_to_cp2k
14 : USE ieee_arithmetic, ONLY: ieee_is_finite
15 : USE input_keyword_types, ONLY: keyword_create, &
16 : keyword_release, &
17 : keyword_type
18 : USE input_section_types, ONLY: section_add_keyword, &
19 : section_create, &
20 : section_type, &
21 : section_vals_get_subs_vals, &
22 : section_vals_type, &
23 : section_vals_val_get
24 : USE input_val_types, ONLY: char_t, &
25 : integer_t, &
26 : real_t, &
27 : val_type
28 : USE kinds, ONLY: default_string_length, &
29 : dp
30 : USE nlo_types, ONLY: nlo_frequency_type, &
31 : nlo_type
32 : USE qs_environment_types, ONLY: get_qs_env, &
33 : qs_environment_type
34 : USE string_utilities, ONLY: uppercase
35 : !&>
36 : #include "../base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nlo_read_input'
43 :
44 : PUBLIC :: create_nonlinear_optics_section, read_nlo_input
45 : PUBLIC :: nlo_frequency_case_count, nlo_frequency_case
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Define the nonlinear susceptibility model, conventions, and input parameters.
51 : !> \param section Created NONLINEAR_OPTICS section of the input tree.
52 : ! **************************************************************************************************
53 1462 : SUBROUTINE create_nonlinear_optics_section(section)
54 : TYPE(section_type), POINTER :: section
55 :
56 : CHARACTER(LEN=6) :: frequency_name
57 : INTEGER :: i_frequency
58 : TYPE(keyword_type), POINTER :: keyword
59 :
60 1462 : NULLIFY (keyword)
61 1462 : CPASSERT(.NOT. ASSOCIATED(section))
62 : CALL section_create(section, __LOCATION__, name="NONLINEAR_OPTICS", &
63 : description= &
64 : "Frequency-dependent nonlinear optical response of an insulating "// &
65 : "periodic Kohn-Sham system in the independent-particle "// &
66 : "approximation (no explicit excited-state treatment, in particular"// &
67 : " no TDDFT, GW-BSE or excitonic effects). The induced polarization"// &
68 : " is expanded in powers of the electric field. With $E(t)=\sum_ω "// &
69 : "E(ω)e^{-iωt}"//"$, the second- through fifth-order "// &
70 : "contributions are "// &
71 : ACHAR(10)//ACHAR(10)// &
72 : "$$P_i^{(2)}(Ω)=ε_0\sum_{ab}\sum_{ω_1+ω_2=Ω} "// &
73 : "χ_{i;ab}^{(2)}(ω_1,ω_2)E_a(ω_1)E_b(ω_2),$$"// &
74 : ACHAR(10)//ACHAR(10)// &
75 : ACHAR(10)//ACHAR(10)// &
76 : "$$P_i^{(3)}(Ω)=ε_0\sum_{abc} \sum_{ω_1+ω_2+ω_3=Ω} "// &
77 : "χ_{i;abc}^{(3)}(ω_1,ω_2,ω_3) E_a(ω_1)E_b(ω_2)E_c(ω_3),$$"// &
78 : ACHAR(10)// &
79 : ACHAR(10)//ACHAR(10)//ACHAR(10)// &
80 : "$$P_i^{(4)}(Ω)=ε_0\sum_{abcd} \sum_{ω_1+ω_2+ω_3+ω_4=Ω} "// &
81 : "χ_{i;abcd}^{(4)}(ω_1,ω_2,ω_3,ω_4) "// &
82 : "E_a(ω_1)E_b(ω_2)E_c(ω_3)E_d(ω_4),$$"// &
83 : ACHAR(10)//ACHAR(10)//ACHAR(10)// &
84 : ACHAR(10)// &
85 : "$$P_i^{(5)}(Ω)=ε_0\sum_{abcde} \sum_{ω_1+ω_2+ω_3+ω_4+ω_5=Ω} "// &
86 : "χ_{i;abcde}^{(5)}(ω_1,ω_2,ω_3,ω_4,ω_5) "// &
87 : "E_a(ω_1)E_b(ω_2)E_c(ω_3)E_d(ω_4)E_e(ω_5).$$"// &
88 : ACHAR(10)//ACHAR(10)// &
89 : "Here i,a,... denote Cartesian directions and ε_0 is the vacuum "// &
90 : "permittivity. Each tensor is averaged over all simultaneous "// &
91 : "permutations of its incident frequency-direction pairs. The "// &
92 : "frequency sums include ordered tuples, so no extra factorial "// &
93 : "occurs in these formulas. Negative frequencies describe "// &
94 : "complex-conjugate fields; harmonic generation uses equal positive"// &
95 : " frequencies. The input specifies photon energies ℏω_j in eV, "// &
96 : "with Ω≠0. The density matrix is expanded recursively in the "// &
97 : "fields, using all AO bands and analytic derivatives with respect "// &
98 : "to k. In the energy eigenbasis the usual "// &
99 : "population and coherence relaxation terms are "//ACHAR(10)//ACHAR(10)// &
100 : "$$\left.\dot{ρ}_{mn}\right|_{\mathrm{rel}}= "// &
101 : "-δ_{mn}(ρ_{nn}-f_n)/T_1 -(1-δ_{mn})ρ_{mn}/T_2.$$"//ACHAR(10)// &
102 : ACHAR(10)// &
103 : "Here $f_n$ is the equilibrium occupation. With T1_T2, "// &
104 : "$T_1=T_2=T$ and the relaxation is "// &
105 : "$-(ρ-ρ_{\mathrm{eq}})/T$. With separate T1 and T2, "// &
106 : "the rate for bands m and n changes smoothly with their energy "// &
107 : "difference Δε: "//ACHAR(10)//ACHAR(10)// &
108 : "$$γ_{mn}=1/T_1+(1/T_2-1/T_1) "// &
109 : "[1-\exp(-(Δε/w)^2)],\qquad Δε=ε_m-ε_n.$$"//ACHAR(10)//ACHAR(10)// &
110 : "The width w is T1_T2_SMOOTHING_WIDTH. At exact degeneracy the "// &
111 : "rate is 1/T1 for all matrices within the degenerate band "// &
112 : "subspace. The Gaussian function of the Hamiltonian commutator "// &
113 : "gives a basis-covariant transition without band grouping. "// &
114 : "At each response step, denominators are "// &
115 : "$ε_m-ε_n-ℏΩ-iγ_{mn}$. The "// &
116 : "current includes the field-dependent contribution "// &
117 : "required by the projected position operator. In "// &
118 : "atomic units, with r_i=i∂_{k_i}+A_i, "//ACHAR(10)//ACHAR(10)// &
119 : "$$F_{ia}=∂_{k_i}A_a-∂_{k_a}A_i-i[A_i,A_a], \qquad "// &
120 : "v_i(E)=i[H+E_ar_a, r_i]=v_i(0)+E_aF_{ia}.$$"//ACHAR(10)//ACHAR(10)// &
121 : "Let I^{(n)} denote the cell-integrated current response, "// &
122 : "including the Brillouin-zone weights and spin "// &
123 : "occupations, in SI units. The volume normalization and "// &
124 : "conversion to susceptibility are "//ACHAR(10)//ACHAR(10)// &
125 : "$$V=|\det(\mathbf a_1,\mathbf a_2,\mathbf a_3)|, \qquad "// &
126 : "j^{(n)}=I^{(n)}/V,\qquad χ_{\mathrm{3D}}^{(n)}= "// &
127 : "σ_{\mathrm{3D}}^{(n)}/(-iΩε_0),$$"//ACHAR(10)//ACHAR(10)// &
128 : "where $σ^{(n)}"//"$ is the coefficient of the n field factors in "// &
129 : "$j^{(n)}"//"$. For a two-dimensional slab with vacuum along "// &
130 : "nonperiodic z, the following sheet quantities remove the "// &
131 : "dependence on the chosen vacuum height: "//ACHAR(10)//ACHAR(10)// &
132 : "$$A=|\mathbf a_1\times\mathbf a_2|,\qquad L_z=V/A, \qquad "// &
133 : "P_{\mathrm{2D}}^{(n)}=L_zP_{\mathrm{3D}}^{(n)}, \qquad "// &
134 : "χ_{\mathrm{2D}}^{(n)}=L_z χ_{\mathrm{3D}}^{(n)}.$$"//ACHAR(10)// &
135 : ACHAR(10)// &
136 : "The units are $(m/V)^{n-1}"//"$ for $χ_{3D}"//"$ and "// &
137 : "$m^n/V^{n-1}"//"$ for "// &
138 : "$χ_{2D}"//"$. Converge the basis set, the k mesh in the SCF "// &
139 : "calculation, and the k mesh used for the NLO susceptibility "// &
140 : "separately. Higher-order susceptibilities "// &
141 : "usually require large optical k meshes.", &
142 1462 : n_keywords=10, n_subsections=0, repeats=.FALSE.)
143 :
144 : CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
145 : description= &
146 : "Activates the nonlinear optical susceptibility calculation.", &
147 1462 : default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
148 1462 : CALL section_add_keyword(section, keyword)
149 1462 : CALL keyword_release(keyword)
150 :
151 : CALL keyword_create(keyword, __LOCATION__, name="ORDER", &
152 : description= &
153 : "Order n of the electric-field response $χ^{(n)}"//"$; "// &
154 : "supported values are 2, 3, 4, and 5. Exactly "// &
155 : "FREQ_1 through FREQ_n must be set.", &
156 1462 : usage="ORDER 2", type_of_var=integer_t, n_var=1, default_i_val=2)
157 1462 : CALL section_add_keyword(section, keyword)
158 1462 : CALL keyword_release(keyword)
159 :
160 : CALL keyword_create(keyword, __LOCATION__, name="ELEMENTS", &
161 : description= &
162 : "Susceptibility elements to compute and print. Use -1 alone for "// &
163 : "all elements (default), or a list of Cartesian index strings. "// &
164 : "Each string has ORDER+1 letters X, Y, or Z: the output direction "// &
165 : "followed by the input directions in FREQ_1 through FREQ_ORDER "// &
166 : "order. For example, YXX selects $χ^{(2)}_{y;xx}"//"$, and XYYYX "// &
167 : "selects $χ^{(4)}_{x;yyyx}"//"$. Uppercase and "// &
168 : "lowercase are accepted."// &
169 : " Repeat the keyword to continue a long list. Duplicate entries "// &
170 : "are computed once. Only selected elements and their required "// &
171 : "intermediate responses are evaluated; unselected "// &
172 : "elements are omitted from the output. Do not "// &
173 : "combine -1 with explicit elements.", &
174 : usage="ELEMENTS YXX YYY", type_of_var=char_t, n_var=-1, &
175 2924 : repeats=.TRUE., default_c_vals=["-1"])
176 1462 : CALL section_add_keyword(section, keyword)
177 1462 : CALL keyword_release(keyword)
178 :
179 8772 : DO i_frequency = 1, 5
180 7310 : WRITE (frequency_name, '(A,I1)') "FREQ_", i_frequency
181 : CALL keyword_create(keyword, __LOCATION__, name=TRIM(frequency_name), &
182 : description= &
183 : "Input photon energy in eV. VALUE specifies one energy; RANGE "// &
184 : "gives inclusive start, end, and positive step; SAME_AS and "// &
185 : "MINUS_OF reference an earlier input frequency. "// &
186 : "Independent ranges produce every combination. The sum "// &
187 : "of input energies must be nonzero.", &
188 : usage=TRIM(frequency_name)//" RANGE 0.5 1.0 0.01", &
189 7310 : type_of_var=char_t, n_var=-1)
190 7310 : CALL section_add_keyword(section, keyword)
191 8772 : CALL keyword_release(keyword)
192 : END DO
193 :
194 : CALL keyword_create(keyword, __LOCATION__, name="T1_T2", &
195 : description= &
196 : "Common relaxation time $T_1=T_2=T$. Usually T₁ relaxes band "// &
197 : "occupations (the diagonal density-matrix elements $ρ_{nn}"//"$) "// &
198 : "towards equilibrium; T₂ describes decoherence, damping "// &
199 : "off-diagonal elements $ρ_{mn}"//"$ to zero in the energy "// &
200 : "eigenbasis. "// &
201 : "Using one time gives $-(ρ-ρ_{eq})/T$ and removes the "// &
202 : "basis-dependent distinction between populations and coherences "// &
203 : "inside degenerate subspaces. This provides a simple, numerically "// &
204 : "robust first approximation near band degeneracies and close "// &
205 : "bands. It does not establish that the physical lifetimes are "// &
206 : "equal. For relaxation-dependent singularities near the graphene "// &
207 : "Dirac point, see Cheng, Vermeulen and Sipe (2015). The "// &
208 : "broadening is $Γ=ℏ/T_2$. Must be positive; the default is 30 fs.", &
209 : usage="T1_T2 [fs] 30", type_of_var=real_t, n_var=1, unit_str="fs", &
210 : default_r_val=cp_unit_to_cp2k(value=30.0_dp, unit_str="fs"), &
211 2924 : citations=[Cheng2015])
212 1462 : CALL section_add_keyword(section, keyword)
213 1462 : CALL keyword_release(keyword)
214 :
215 : CALL keyword_create(keyword, __LOCATION__, name="T1", &
216 : description= &
217 : "Relaxation time for diagonal band occupations. Set together with "// &
218 : "T2 to use separate population and coherence times. Must satisfy "// &
219 : "T1 >= T2 > 0. Do not combine with an explicitly set T1_T2.", &
220 1462 : usage="T1 [fs] 300", type_of_var=real_t, n_var=1, unit_str="fs")
221 1462 : CALL section_add_keyword(section, keyword)
222 1462 : CALL keyword_release(keyword)
223 :
224 : CALL keyword_create(keyword, __LOCATION__, name="T2", &
225 : description= &
226 : "Decoherence time for well-separated bands. Set together with T1. "// &
227 : "Near degenerate bands the relaxation rate changes smoothly to 1/T1.", &
228 1462 : usage="T2 [fs] 30", type_of_var=real_t, n_var=1, unit_str="fs")
229 1462 : CALL section_add_keyword(section, keyword)
230 1462 : CALL keyword_release(keyword)
231 :
232 : CALL keyword_create(keyword, __LOCATION__, name="T1_T2_SMOOTHING_WIDTH", &
233 : description= &
234 : "Positive energy width w of the smooth transition between population "// &
235 : "and coherence relaxation. The band-pair rate is "// &
236 : "1/T1 + (1/T2 - 1/T1)[1 - exp(-((epsilon_m-epsilon_n)/w)^2)]. "// &
237 : "The default 0.025 eV is a model choice to converge for each system.", &
238 : usage="T1_T2_SMOOTHING_WIDTH [eV] 0.025", type_of_var=real_t, &
239 : n_var=1, unit_str="eV", &
240 1462 : default_r_val=cp_unit_to_cp2k(value=0.025_dp, unit_str="eV"))
241 1462 : CALL section_add_keyword(section, keyword)
242 1462 : CALL keyword_release(keyword)
243 :
244 : CALL keyword_create(keyword, __LOCATION__, name="K_MESH", &
245 : description= &
246 : "Monkhorst-Pack-type uniform k mesh, shifted to include Γ, for the"// &
247 : " Brillouin-zone integration required for $χ^{(n)}"//"$. Fractional "// &
248 : "coordinates are k_i=j_i/N_i with j_i=0,...,N_i-1; no symmetry "// &
249 : "reduction is used. This mesh can differ from the required k-point"// &
250 : " SCF mesh. Higher-order nonlinear susceptibilities usually "// &
251 : "require denser k meshes. Every entry must be "// &
252 : "positive; use one in each nonperiodic direction.", &
253 1462 : usage="K_MESH 24 24 1", type_of_var=integer_t, n_var=3)
254 1462 : CALL section_add_keyword(section, keyword)
255 1462 : CALL keyword_release(keyword)
256 :
257 1462 : END SUBROUTINE create_nonlinear_optics_section
258 :
259 : ! **************************************************************************************************
260 : !> \brief Read the activated input into nlo_env and reject invalid frequency combinations.
261 : !> \param qs_env Converged QS environment containing the parsed input.
262 : !> \param nlo_env Validated settings for the nonlinear response.
263 : ! **************************************************************************************************
264 242 : SUBROUTINE read_nlo_input(qs_env, nlo_env)
265 : TYPE(qs_environment_type), POINTER :: qs_env
266 : TYPE(nlo_type), INTENT(OUT) :: nlo_env
267 :
268 : CHARACTER(len=*), PARAMETER :: routineN = 'read_nlo_input'
269 :
270 : INTEGER :: handle, i_case
271 : LOGICAL :: zero_output
272 : REAL(KIND=dp) :: energy_ev(5)
273 : TYPE(section_vals_type), POINTER :: input, section
274 :
275 22 : CALL timeset(routineN, handle)
276 :
277 22 : CALL get_qs_env(qs_env, input=input)
278 22 : section => section_vals_get_subs_vals(input, 'PROPERTIES%NONLINEAR_OPTICS')
279 22 : CALL read_nlo_values(section, nlo_env)
280 46 : DO i_case = 1, nlo_env%n_frequency_combinations
281 24 : CALL nlo_frequency_case(nlo_env%frequency, nlo_env%order, i_case, energy_ev, zero_output)
282 46 : IF (zero_output) CPABORT('NLO output frequency must be nonzero.')
283 : END DO
284 :
285 22 : CALL timestop(handle)
286 :
287 22 : END SUBROUTINE read_nlo_input
288 :
289 : ! **************************************************************************************************
290 : !> \brief Read and validate the NONLINEAR_OPTICS input before the first response calculation.
291 : !> \param section Activated NONLINEAR_OPTICS input section.
292 : !> \param nlo_env Validated response settings and selected elements.
293 : ! **************************************************************************************************
294 242 : SUBROUTINE read_nlo_values(section, nlo_env)
295 : TYPE(section_vals_type), POINTER :: section
296 : TYPE(nlo_type), INTENT(OUT) :: nlo_env
297 :
298 : CHARACTER(len=6) :: frequency_name
299 22 : CHARACTER(len=default_string_length), POINTER :: words(:)
300 : INTEGER :: i_frequency
301 : LOGICAL :: common_explicit, t1_explicit, &
302 : t2_explicit
303 22 : INTEGER, POINTER :: mesh(:)
304 : TYPE(val_type), POINTER :: input_value
305 :
306 220 : nlo_env = nlo_type()
307 22 : IF (.NOT. ASSOCIATED(section)) CPABORT('NLO input section is missing.')
308 22 : NULLIFY (input_value, mesh, words)
309 22 : CALL section_vals_val_get(section, 'ORDER', i_val=nlo_env%order)
310 22 : IF (nlo_env%order < 2 .OR. nlo_env%order > 5) THEN
311 0 : CPABORT('NLO ORDER must be between 2 and 5.')
312 : END IF
313 22 : CALL read_nlo_elements(section, nlo_env)
314 22 : CALL section_vals_val_get(section, 'T1_T2', r_val=nlo_env%t1_au)
315 22 : IF (.NOT. ieee_is_finite(nlo_env%t1_au) .OR. nlo_env%t1_au <= 0.0_dp) THEN
316 0 : CPABORT('NONLINEAR_OPTICS T1_T2 must be a positive, finite time.')
317 : END IF
318 22 : CALL section_vals_val_get(section, 'T1_T2', explicit=common_explicit)
319 22 : CALL section_vals_val_get(section, 'T1', explicit=t1_explicit)
320 22 : CALL section_vals_val_get(section, 'T2', explicit=t2_explicit)
321 22 : IF (t1_explicit .NEQV. t2_explicit) THEN
322 0 : CPABORT('NONLINEAR_OPTICS requires both T1 and T2 when either is set.')
323 : END IF
324 22 : IF (t1_explicit .AND. common_explicit) THEN
325 0 : CPABORT('NONLINEAR_OPTICS cannot combine T1_T2 with T1 and T2.')
326 : END IF
327 22 : nlo_env%t2_au = nlo_env%t1_au
328 22 : IF (t1_explicit) THEN
329 4 : CALL section_vals_val_get(section, 'T1', r_val=nlo_env%t1_au)
330 4 : CALL section_vals_val_get(section, 'T2', r_val=nlo_env%t2_au)
331 : END IF
332 22 : IF (.NOT. ieee_is_finite(nlo_env%t1_au) .OR. &
333 : .NOT. ieee_is_finite(nlo_env%t2_au)) THEN
334 0 : CPABORT('NONLINEAR_OPTICS T1 and T2 must be finite times.')
335 : END IF
336 22 : IF (nlo_env%t2_au <= 0.0_dp .OR. nlo_env%t1_au < nlo_env%t2_au) THEN
337 0 : CPABORT('NONLINEAR_OPTICS requires T1 >= T2 > 0.')
338 : END IF
339 : CALL section_vals_val_get(section, 'T1_T2_SMOOTHING_WIDTH', &
340 22 : r_val=nlo_env%transition_width_au)
341 22 : IF (.NOT. ieee_is_finite(nlo_env%transition_width_au) .OR. &
342 : nlo_env%transition_width_au <= 0.0_dp) THEN
343 0 : CPABORT('NONLINEAR_OPTICS requires a positive, finite transition width.')
344 : END IF
345 22 : CALL section_vals_val_get(section, 'K_MESH', val=input_value)
346 22 : IF (.NOT. ASSOCIATED(input_value)) CPABORT('NONLINEAR_OPTICS requires K_MESH.')
347 22 : CALL section_vals_val_get(section, 'K_MESH', i_vals=mesh)
348 22 : IF (.NOT. ASSOCIATED(mesh)) CPABORT('NLO requires three K_MESH entries.')
349 22 : IF (SIZE(mesh) /= 3) CPABORT('NLO requires three K_MESH entries.')
350 88 : IF (ANY(mesh <= 0)) CPABORT('NONLINEAR_OPTICS K_MESH entries must be positive.')
351 88 : nlo_env%k_mesh = mesh
352 132 : DO i_frequency = 1, 5
353 110 : WRITE (frequency_name, '(A,I1)') 'FREQ_', i_frequency
354 110 : NULLIFY (input_value)
355 110 : CALL section_vals_val_get(section, frequency_name, val=input_value)
356 110 : IF (i_frequency > nlo_env%order) THEN
357 40 : IF (ASSOCIATED(input_value)) THEN
358 0 : CPABORT('Do not set input frequencies beyond NONLINEAR_OPTICS ORDER.')
359 : END IF
360 : CYCLE
361 : END IF
362 70 : IF (.NOT. ASSOCIATED(input_value)) CPABORT('Set all FREQ_1 through FREQ_ORDER.')
363 70 : NULLIFY (words)
364 70 : CALL section_vals_val_get(section, frequency_name, c_vals=words)
365 70 : IF (.NOT. ASSOCIATED(words)) CPABORT('NLO frequency specification is empty.')
366 92 : CALL parse_nlo_frequency(words, i_frequency, nlo_env%frequency(i_frequency))
367 : END DO
368 : ! Number of independent frequency tuples: product_j N(FREQ_j).
369 : CALL nlo_frequency_case_count(nlo_env%frequency, nlo_env%order, &
370 22 : nlo_env%n_frequency_combinations)
371 :
372 22 : END SUBROUTINE read_nlo_values
373 :
374 : ! **************************************************************************************************
375 : !> \brief Combine repeated ELEMENTS lists without imposing a limit on their length.
376 : !> \param section Input section containing the ELEMENTS keywords.
377 : !> \param nlo_env Response settings receiving the selected elements.
378 : ! **************************************************************************************************
379 22 : SUBROUTINE read_nlo_elements(section, nlo_env)
380 : TYPE(section_vals_type), POINTER :: section
381 : TYPE(nlo_type), INTENT(INOUT) :: nlo_env
382 :
383 22 : CHARACTER(len=default_string_length), POINTER :: words(:)
384 : INTEGER :: i_rep, n_rep
385 : LOGICAL, ALLOCATABLE :: elements(:, :)
386 :
387 88 : ALLOCATE (nlo_env%elements(3, 3**nlo_env%order), elements(3, 3**nlo_env%order))
388 7294 : nlo_env%elements = .FALSE.
389 22 : CALL section_vals_val_get(section, 'ELEMENTS', n_rep_val=n_rep)
390 44 : DO i_rep = 1, MAX(1, n_rep)
391 22 : NULLIFY (words)
392 22 : CALL section_vals_val_get(section, 'ELEMENTS', i_rep_val=i_rep, c_vals=words)
393 22 : IF (.NOT. ASSOCIATED(words)) CPABORT('NLO ELEMENTS list is empty.')
394 32 : IF (ANY(words == '-1') .AND. n_rep > 1) THEN
395 0 : CPABORT('NLO ELEMENTS: -1 cannot be repeated.')
396 : END IF
397 22 : CALL parse_nlo_elements(words, nlo_env%order, elements)
398 7316 : nlo_env%elements(:, :) = nlo_env%elements .OR. elements
399 : END DO
400 :
401 66 : END SUBROUTINE read_nlo_elements
402 :
403 : ! **************************************************************************************************
404 : !> \brief Decode Cartesian output/input strings into the selected tensor elements.
405 : !> \param tokens Cartesian tensor-index strings from an ELEMENTS keyword.
406 : !> \param order Number of incident electric fields.
407 : !> \param elements Unique, validated tensor-index selections.
408 : !> \note elements(c,1+sum_j(a_j-1)*3**(j-1)) selects χ_(c;a_1,...,a_order).
409 : ! **************************************************************************************************
410 22 : SUBROUTINE parse_nlo_elements(tokens, order, elements)
411 : CHARACTER(len=*), INTENT(IN) :: tokens(:)
412 : INTEGER, INTENT(IN) :: order
413 : LOGICAL, INTENT(OUT) :: elements(:, :)
414 :
415 : CHARACTER(len=6) :: element
416 : INTEGER :: i_code, i_direction, i_letter, i_out, &
417 : i_token
418 :
419 7294 : elements = .FALSE.
420 22 : IF (order < 2 .OR. order > 5) THEN
421 0 : CPABORT('NONLINEAR_OPTICS ELEMENTS requires ORDER 2 through 5.')
422 : END IF
423 66 : IF (ANY(SHAPE(elements) /= [3, 3**order]) .OR. SIZE(tokens) == 0) THEN
424 0 : CPABORT('NLO ELEMENTS list is empty.')
425 : END IF
426 32 : IF (ANY(tokens == '-1')) THEN
427 12 : IF (SIZE(tokens) /= 1) CPABORT('NLO ELEMENTS: -1 must appear alone.')
428 588 : elements = .TRUE.
429 12 : RETURN
430 : END IF
431 20 : DO i_token = 1, SIZE(tokens)
432 10 : IF (LEN_TRIM(tokens(i_token)) /= order + 1) THEN
433 0 : CPABORT('NLO ELEMENTS: use ORDER+1 letters, output first.')
434 : END IF
435 10 : element = tokens(i_token)
436 10 : CALL uppercase(element)
437 10 : i_out = 0
438 10 : i_code = 1
439 64 : DO i_letter = 1, order + 1
440 54 : i_direction = INDEX('XYZ', element(i_letter:i_letter))
441 54 : IF (i_direction == 0) THEN
442 0 : CPABORT('NONLINEAR_OPTICS ELEMENTS: use only X, Y, Z in each element.')
443 : END IF
444 64 : IF (i_letter == 1) THEN
445 : i_out = i_direction
446 : ELSE
447 44 : i_code = i_code + (i_direction - 1)*3**(i_letter - 2)
448 : END IF
449 : END DO
450 20 : elements(i_out, i_code) = .TRUE.
451 : END DO
452 :
453 : END SUBROUTINE parse_nlo_elements
454 :
455 : ! **************************************************************************************************
456 : !> \brief Parse VALUE, RANGE, SAME_AS or MINUS_OF for input frequency i_frequency.
457 : !> \param tokens Values of the incident-frequency keyword.
458 : !> \param i_frequency One-based index of the incident frequency.
459 : !> \param frequency Parsed single, ranged, or linked incident frequency.
460 : ! **************************************************************************************************
461 70 : SUBROUTINE parse_nlo_frequency(tokens, i_frequency, frequency)
462 : CHARACTER(len=*), INTENT(IN) :: tokens(:)
463 : INTEGER, INTENT(IN) :: i_frequency
464 : TYPE(nlo_frequency_type), INTENT(OUT) :: frequency
465 :
466 : CHARACTER(len=32) :: mode, reference_name
467 : INTEGER :: i_reference, io_error
468 : REAL(KIND=dp) :: end_ev, interval_count
469 :
470 : frequency = nlo_frequency_type()
471 70 : IF (i_frequency < 1 .OR. i_frequency > 5 .OR. SIZE(tokens) < 1) THEN
472 0 : CPABORT('Invalid FREQ index or empty frequency specification.')
473 : END IF
474 70 : mode = tokens(1)
475 70 : CALL uppercase(mode)
476 92 : SELECT CASE (TRIM(mode))
477 : CASE ('VALUE')
478 22 : IF (SIZE(tokens) /= 2) CPABORT('FREQ VALUE needs one energy in eV.')
479 22 : READ (tokens(2), *, IOSTAT=io_error) frequency%start_ev
480 22 : IF (io_error /= 0) CPABORT('FREQ VALUE requires a finite photon energy in eV.')
481 22 : IF (.NOT. ieee_is_finite(frequency%start_ev)) THEN
482 0 : CPABORT('FREQ VALUE requires a finite photon energy in eV.')
483 : END IF
484 22 : frequency%n_values = 1
485 : CASE ('RANGE')
486 2 : IF (SIZE(tokens) /= 4) CPABORT('FREQ RANGE requires start, end and step in eV.')
487 2 : READ (tokens(2), *, IOSTAT=io_error) frequency%start_ev
488 2 : IF (io_error /= 0) CPABORT('Invalid FREQ RANGE start energy.')
489 2 : READ (tokens(3), *, IOSTAT=io_error) end_ev
490 2 : IF (io_error /= 0) CPABORT('Invalid FREQ RANGE end energy.')
491 2 : READ (tokens(4), *, IOSTAT=io_error) frequency%step_ev
492 2 : IF (io_error /= 0) CPABORT('Invalid FREQ RANGE step energy.')
493 2 : IF (.NOT. ieee_is_finite(frequency%start_ev) .OR. .NOT. ieee_is_finite(end_ev) .OR. &
494 : .NOT. ieee_is_finite(frequency%step_ev)) THEN
495 0 : CPABORT('FREQ RANGE energies and step must be finite.')
496 : END IF
497 2 : IF (frequency%step_ev <= 0.0_dp .OR. end_ev < frequency%start_ev) THEN
498 0 : CPABORT('FREQ RANGE requires positive step and end >= start.')
499 : END IF
500 2 : interval_count = (end_ev - frequency%start_ev)/frequency%step_ev
501 : IF (.NOT. ieee_is_finite(interval_count) &
502 2 : .OR. interval_count >= REAL(HUGE(frequency%n_values) - 1, KIND=dp)) THEN
503 0 : CPABORT('FREQ RANGE contains too many values.')
504 : END IF
505 : frequency%n_values = 1 + FLOOR(interval_count + &
506 2 : 32.0_dp*EPSILON(interval_count)*MAX(1.0_dp, interval_count))
507 : CASE ('SAME_AS', 'MINUS_OF')
508 46 : IF (SIZE(tokens) /= 2) CPABORT('FREQ relation needs one earlier index.')
509 46 : reference_name = tokens(2)
510 46 : CALL uppercase(reference_name)
511 46 : IF (LEN_TRIM(reference_name) /= 6 .OR. reference_name(1:5) /= 'FREQ_') THEN
512 0 : CPABORT('Frequency reference must be FREQ_1, FREQ_2, etc.')
513 : END IF
514 46 : READ (reference_name(6:6), *, IOSTAT=io_error) i_reference
515 46 : IF (io_error /= 0) CPABORT('Invalid frequency reference index.')
516 46 : IF (i_reference < 1 .OR. i_reference >= i_frequency) THEN
517 0 : CPABORT('Frequency references must point to an earlier FREQ.')
518 : END IF
519 46 : frequency%reference = i_reference
520 46 : IF (TRIM(mode) == 'MINUS_OF') frequency%sign = -1
521 : CASE DEFAULT
522 70 : CPABORT('FREQ must use VALUE, RANGE, SAME_AS or MINUS_OF.')
523 : END SELECT
524 :
525 70 : END SUBROUTINE parse_nlo_frequency
526 :
527 : ! **************************************************************************************************
528 : !> \brief Count independent frequency combinations without materializing the Cartesian product.
529 : !> \param frequencies Parsed incident-frequency definitions.
530 : !> \param order Number of active incident frequencies.
531 : !> \param n_frequency_combinations Product of independent range lengths.
532 : ! **************************************************************************************************
533 83 : SUBROUTINE nlo_frequency_case_count(frequencies, order, n_frequency_combinations)
534 : TYPE(nlo_frequency_type), INTENT(IN) :: frequencies(5)
535 : INTEGER, INTENT(IN) :: order
536 : INTEGER, INTENT(OUT) :: n_frequency_combinations
537 :
538 : INTEGER :: i_frequency
539 :
540 83 : n_frequency_combinations = 1
541 83 : IF (order < 2 .OR. order > 5) CPABORT('NLO ORDER must be between 2 and 5.')
542 341 : DO i_frequency = 1, order
543 341 : IF (frequencies(i_frequency)%reference > 0) THEN
544 167 : IF (frequencies(i_frequency)%reference >= i_frequency) THEN
545 0 : CPABORT('Frequency reference must point to an earlier FREQ.')
546 : END IF
547 : ELSE
548 91 : IF (frequencies(i_frequency)%n_values < 1) THEN
549 0 : CPABORT('Every input frequency through ORDER must be specified.')
550 : END IF
551 91 : IF (n_frequency_combinations > &
552 : HUGE(n_frequency_combinations)/frequencies(i_frequency)%n_values) THEN
553 0 : CPABORT('Too many independent frequency combinations.')
554 : END IF
555 91 : n_frequency_combinations = n_frequency_combinations*frequencies(i_frequency)%n_values
556 : END IF
557 : END DO
558 :
559 83 : END SUBROUTINE nlo_frequency_case_count
560 :
561 : ! **************************************************************************************************
562 : !> \brief Resolve one Cartesian-product case, preserving SAME_AS and MINUS_OF links.
563 : !> \param frequencies Parsed incident-frequency definitions.
564 : !> \param order Number of active incident frequencies.
565 : !> \param i_case One-based index of the requested range combination.
566 : !> \param energy_ev Resolved incident photon energies in eV.
567 : !> \param zero_output Whether the sum of incident photon energies is zero.
568 : ! **************************************************************************************************
569 61 : SUBROUTINE nlo_frequency_case(frequencies, order, i_case, energy_ev, zero_output)
570 : TYPE(nlo_frequency_type), INTENT(IN) :: frequencies(5)
571 : INTEGER, INTENT(IN) :: order, i_case
572 : REAL(KIND=dp), INTENT(OUT) :: energy_ev(5)
573 : LOGICAL, INTENT(OUT) :: zero_output
574 :
575 : INTEGER :: case_offset, i_frequency, i_value, &
576 : n_frequency_combinations
577 : REAL(KIND=dp) :: scale
578 :
579 61 : energy_ev = 0.0_dp
580 61 : zero_output = .FALSE.
581 61 : CALL nlo_frequency_case_count(frequencies, order, n_frequency_combinations)
582 61 : IF (i_case < 1 .OR. i_case > n_frequency_combinations) THEN
583 0 : CPABORT('Frequency case index is out of range.')
584 : END IF
585 61 : case_offset = i_case - 1
586 249 : DO i_frequency = 1, order
587 249 : IF (frequencies(i_frequency)%reference > 0) THEN
588 : energy_ev(i_frequency) &
589 121 : = frequencies(i_frequency)%sign*energy_ev(frequencies(i_frequency)%reference)
590 : ELSE
591 67 : i_value = MOD(case_offset, frequencies(i_frequency)%n_values)
592 67 : case_offset = case_offset/frequencies(i_frequency)%n_values
593 : energy_ev(i_frequency) = frequencies(i_frequency)%start_ev &
594 67 : + i_value*frequencies(i_frequency)%step_ev
595 : END IF
596 : END DO
597 249 : scale = MAX(1.0_dp, SUM(ABS(energy_ev(1:order))))
598 249 : zero_output = ABS(SUM(energy_ev(1:order))) <= 64.0_dp*EPSILON(scale)*scale
599 :
600 61 : END SUBROUTINE nlo_frequency_case
601 :
602 : END MODULE nlo_read_input
|
