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 Utilities for frontier-orbital basis optimization
9 : !> \par History
10 : !> 28.08.2026 created [Jan Wilhelm]
11 : !> \author Jan Wilhelm
12 : ! **************************************************************************************************
13 : MODULE optbas_frontier_orbitals_utils
14 : USE basis_set_types, ONLY: get_gto_basis_set,&
15 : gto_basis_set_type
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
19 : dbcsr_p_type,&
20 : dbcsr_type
21 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
22 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
23 : USE cp_fm_diag, ONLY: cp_fm_geeig
24 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 : cp_fm_struct_release,&
26 : cp_fm_struct_type
27 : USE cp_fm_types, ONLY: cp_fm_create,&
28 : cp_fm_get_submatrix,&
29 : cp_fm_release,&
30 : cp_fm_set_all,&
31 : cp_fm_set_element,&
32 : cp_fm_type
33 : USE ieee_arithmetic, ONLY: ieee_is_finite
34 : USE kinds, ONLY: default_string_length,&
35 : dp
36 : USE parallel_gemm_api, ONLY: parallel_gemm
37 : USE qs_condnum, ONLY: overlap_condnum
38 : USE qs_environment_types, ONLY: get_qs_env,&
39 : qs_environment_type
40 : USE qs_kind_types, ONLY: get_qs_kind,&
41 : qs_kind_type
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 : PRIVATE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optbas_frontier_orbitals_utils'
48 :
49 : TYPE, PUBLIC :: frontier_orbitals_result_type
50 : REAL(KIND=dp) :: loss_occupied_subspace = 0.0_dp
51 : REAL(KIND=dp) :: loss_virtual_subspace = 0.0_dp
52 : REAL(KIND=dp) :: loss_empty_subspace = 0.0_dp
53 : REAL(KIND=dp) :: loss_gap = 0.0_dp
54 : REAL(KIND=dp) :: occupied_subspace_overlap = 0.0_dp
55 : REAL(KIND=dp) :: virtual_subspace_similarity = 0.0_dp
56 : REAL(KIND=dp) :: empty_subspace_overlap = 0.0_dp
57 : REAL(KIND=dp) :: gap_reference = 0.0_dp
58 : REAL(KIND=dp) :: gap_candidate = 0.0_dp
59 : REAL(KIND=dp) :: condition_number = 0.0_dp
60 : REAL(KIND=dp) :: reference_condition_number = 0.0_dp
61 : REAL(KIND=dp) :: calculation_time = 0.0_dp
62 : INTEGER :: number_candidate_orbitals = 0
63 : INTEGER :: number_reference_orbitals = 0
64 : END TYPE frontier_orbitals_result_type
65 :
66 : TYPE, PUBLIC :: frontier_orbitals_reference_type
67 : INTEGER :: number_occupied_orbitals = 0
68 : INTEGER :: number_reference_orbitals = 0
69 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues
70 : TYPE(cp_fm_type) :: eigenvectors
71 : REAL(KIND=dp) :: condition_number = 0.0_dp
72 : INTEGER :: number_kinds = 0
73 : CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: element_symbol
74 : CHARACTER(LEN=default_string_length), DIMENSION(:), ALLOCATABLE :: reference_basis_name
75 : CHARACTER(LEN=default_string_length), DIMENSION(:), ALLOCATABLE :: candidate_basis_name
76 : END TYPE frontier_orbitals_reference_type
77 :
78 : PUBLIC :: compute_frontier_orbital_loss_function, &
79 : evaluate_frontier_orbitals_objective, &
80 : frontier_orbitals_reference_init, frontier_orbitals_reference_release, &
81 : print_frontier_orbital_basis_sets, &
82 : smooth_frontier_orbital_weight
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief Print reference and initial basis-set names for every fitted atom kind.
88 : !> \param unit_nr output unit
89 : !> \param element_symbols fitted atom-kind symbols
90 : !> \param reference_basis_names reference basis-set names
91 : !> \param initial_basis_names initial small-basis names
92 : ! **************************************************************************************************
93 2 : SUBROUTINE print_frontier_orbital_basis_sets(unit_nr, element_symbols, reference_basis_names, &
94 2 : initial_basis_names)
95 : INTEGER, INTENT(IN) :: unit_nr
96 : CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: element_symbols, reference_basis_names, &
97 : initial_basis_names
98 :
99 : INTEGER :: ikind
100 :
101 2 : CPASSERT(SIZE(element_symbols) == SIZE(reference_basis_names))
102 2 : CPASSERT(SIZE(element_symbols) == SIZE(initial_basis_names))
103 2 : WRITE (unit_nr, '(1X,A,T14,A,T51,A)') "Atom kind", "Ref. basis set", "Initial basis set"
104 6 : DO ikind = 1, SIZE(element_symbols)
105 4 : WRITE (unit_nr, '(1X,A,T14,A,T51,A)') TRIM(element_symbols(ikind)), &
106 10 : TRIM(reference_basis_names(ikind)), TRIM(initial_basis_names(ikind))
107 : END DO
108 :
109 2 : END SUBROUTINE print_frontier_orbital_basis_sets
110 :
111 : ! **************************************************************************************************
112 : !> \brief Diagonalize and store the frozen reference Hamiltonian.
113 : !> \param qs_env initialized reference Quickstep environment
114 : !> \param reference reference eigenvalues and eigenvectors
115 : ! **************************************************************************************************
116 4 : SUBROUTINE frontier_orbitals_reference_init(qs_env, reference)
117 : TYPE(qs_environment_type), POINTER :: qs_env
118 : TYPE(frontier_orbitals_reference_type), &
119 : INTENT(OUT) :: reference
120 :
121 : INTEGER :: nspins
122 : INTEGER, DIMENSION(2) :: nelectron_spin
123 : LOGICAL :: do_kpoints
124 : REAL(KIND=dp), DIMENSION(2) :: condnum
125 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
126 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
127 : TYPE(cp_fm_type) :: reference_hamiltonian, &
128 : reference_overlap, work
129 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
130 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: smat
131 : TYPE(dft_control_type), POINTER :: dft_control
132 :
133 4 : NULLIFY (blacs_env, dft_control, fm_struct, matrix_ks, matrix_s, smat)
134 : CALL get_qs_env(qs_env, blacs_env=blacs_env, dft_control=dft_control, do_kpoints=do_kpoints, &
135 4 : matrix_ks=matrix_ks, matrix_s=matrix_s, nelectron_spin=nelectron_spin)
136 :
137 4 : nspins = SIZE(matrix_ks)
138 4 : IF (nspins /= 1 .OR. dft_control%multiplicity /= 1) THEN
139 0 : CPABORT("FRONTIER_ORBITALS currently supports only closed-shell calculations")
140 : END IF
141 4 : IF (do_kpoints) THEN
142 0 : CPABORT("FRONTIER_ORBITALS currently supports only Γ-point calculations")
143 : END IF
144 4 : IF (ASSOCIATED(qs_env%x_data)) THEN
145 0 : CPABORT("FRONTIER_ORBITALS does not yet support hybrid functionals")
146 : END IF
147 :
148 : CALL dbcsr_get_info(matrix_s(1)%matrix, &
149 4 : nfullrows_total=reference%number_reference_orbitals)
150 12 : ALLOCATE (smat(1, 1))
151 4 : smat(1, 1)%matrix => matrix_s(1)%matrix
152 4 : CALL overlap_condnum(smat, condnum, -1, .FALSE., .TRUE., .FALSE., blacs_env)
153 4 : reference%condition_number = condnum(2)
154 4 : DEALLOCATE (smat)
155 :
156 4 : CALL store_basis_metadata(qs_env, reference)
157 :
158 : ! A restricted spatial orbital contains two electrons.
159 4 : reference%number_occupied_orbitals = nelectron_spin(1)/2
160 4 : IF (2*reference%number_occupied_orbitals /= nelectron_spin(1)) THEN
161 0 : CPABORT("FRONTIER_ORBITALS requires an even closed-shell electron count")
162 : END IF
163 4 : IF (reference%number_occupied_orbitals < 1 .OR. &
164 : reference%number_occupied_orbitals >= reference%number_reference_orbitals) THEN
165 0 : CPABORT("FRONTIER_ORBITALS requires occupied and virtual reference orbitals")
166 : END IF
167 :
168 12 : ALLOCATE (reference%eigenvalues(reference%number_reference_orbitals))
169 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, &
170 : nrow_global=reference%number_reference_orbitals, &
171 4 : ncol_global=reference%number_reference_orbitals)
172 : CALL cp_fm_create(reference_hamiltonian, fm_struct, &
173 4 : name="frontier_orbitals_reference_hamiltonian")
174 : CALL cp_fm_create(reference_overlap, fm_struct, &
175 4 : name="frontier_orbitals_reference_overlap")
176 : CALL cp_fm_create(reference%eigenvectors, fm_struct, &
177 4 : name="frontier_orbitals_reference_eigenvectors")
178 4 : CALL cp_fm_create(work, fm_struct, name="frontier_orbitals_reference_work")
179 4 : CALL cp_fm_struct_release(fm_struct)
180 :
181 4 : CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, reference_hamiltonian)
182 4 : CALL cp_fm_uplo_to_full(reference_hamiltonian, work)
183 4 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, reference_overlap)
184 4 : CALL cp_fm_uplo_to_full(reference_overlap, work)
185 :
186 : ! H_R C_R = S_R C_R diag(ε_R), with C_R^T S_R C_R = I.
187 : CALL cp_fm_geeig(reference_hamiltonian, reference_overlap, reference%eigenvectors, &
188 4 : reference%eigenvalues, work)
189 :
190 4 : CALL cp_fm_release(reference_hamiltonian)
191 4 : CALL cp_fm_release(reference_overlap)
192 4 : CALL cp_fm_release(work)
193 :
194 12 : END SUBROUTINE frontier_orbitals_reference_init
195 :
196 : ! **************************************************************************************************
197 : !> \brief Release the stored reference eigenvalues and eigenvectors.
198 : !> \param reference reference eigenvalues and eigenvectors
199 : ! **************************************************************************************************
200 4 : SUBROUTINE frontier_orbitals_reference_release(reference)
201 : TYPE(frontier_orbitals_reference_type), &
202 : INTENT(INOUT) :: reference
203 :
204 4 : IF (ALLOCATED(reference%eigenvalues)) DEALLOCATE (reference%eigenvalues)
205 4 : IF (ALLOCATED(reference%element_symbol)) DEALLOCATE (reference%element_symbol)
206 4 : IF (ALLOCATED(reference%reference_basis_name)) DEALLOCATE (reference%reference_basis_name)
207 4 : IF (ALLOCATED(reference%candidate_basis_name)) DEALLOCATE (reference%candidate_basis_name)
208 4 : IF (ASSOCIATED(reference%eigenvectors%matrix_struct)) THEN
209 4 : CALL cp_fm_release(reference%eigenvectors)
210 : END IF
211 4 : reference%number_occupied_orbitals = 0
212 4 : reference%number_reference_orbitals = 0
213 4 : reference%number_kinds = 0
214 4 : reference%condition_number = 0.0_dp
215 :
216 4 : END SUBROUTINE frontier_orbitals_reference_release
217 :
218 : ! **************************************************************************************************
219 : !> \brief Store the reference and candidate basis assignment for every atomic kind.
220 : !> \param qs_env initialized Quickstep environment
221 : !> \param reference reference data receiving the metadata
222 : ! **************************************************************************************************
223 4 : SUBROUTINE store_basis_metadata(qs_env, reference)
224 : TYPE(qs_environment_type), POINTER :: qs_env
225 : TYPE(frontier_orbitals_reference_type), &
226 : INTENT(INOUT) :: reference
227 :
228 : INTEGER :: ikind
229 : TYPE(gto_basis_set_type), POINTER :: basis_set
230 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
231 :
232 4 : NULLIFY (basis_set, qs_kind_set)
233 4 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
234 4 : reference%number_kinds = SIZE(qs_kind_set)
235 8 : ALLOCATE (reference%element_symbol(reference%number_kinds))
236 12 : ALLOCATE (reference%reference_basis_name(reference%number_kinds))
237 8 : ALLOCATE (reference%candidate_basis_name(reference%number_kinds))
238 12 : reference%reference_basis_name = ""
239 12 : reference%candidate_basis_name = ""
240 :
241 12 : DO ikind = 1, reference%number_kinds
242 8 : CALL get_qs_kind(qs_kind_set(ikind), element_symbol=reference%element_symbol(ikind))
243 8 : NULLIFY (basis_set)
244 8 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type="ORB")
245 8 : IF (ASSOCIATED(basis_set)) THEN
246 8 : CALL get_gto_basis_set(basis_set, name=reference%reference_basis_name(ikind))
247 : END IF
248 8 : NULLIFY (basis_set)
249 8 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type="AUX_OPT")
250 12 : IF (ASSOCIATED(basis_set)) THEN
251 0 : CALL get_gto_basis_set(basis_set, name=reference%candidate_basis_name(ikind))
252 : END IF
253 : END DO
254 :
255 4 : END SUBROUTINE store_basis_metadata
256 :
257 : ! **************************************************************************************************
258 : !> \brief Construct and diagonalize the candidate Hamiltonian and evaluate the loss function.
259 : !> \param reference reference eigenvalues and eigenvectors C_R
260 : !> \param matrix_s_candidate candidate AO overlap matrix S_C
261 : !> \param matrix_s_candidate_reference candidate/reference AO overlap matrix S_CR
262 : !> \param virtual_cutoff virtual-orbital energy cutoff
263 : !> \param virtual_smoothing virtual-orbital energy smoothing width
264 : !> \param gap_scale normalization energy for the HOMO-LUMO gap loss
265 : !> \param objective_result loss-function terms and result information
266 : !> \note The reference Hamiltonian is transformed to the candidate AO basis using
267 : !> B = C_R^T S_RC
268 : !> H_C = B^T diag(ε_R) B
269 : !> The candidate orbitals are obtained from
270 : !> H_C C_C = S_C C_C diag(ε_C)
271 : !> C_C^T S_C C_C = I. The reference/candidate orbital overlap matrix is
272 : !> M = C_R^T S_RC C_C = B C_C.
273 : ! **************************************************************************************************
274 408 : SUBROUTINE evaluate_frontier_orbitals_objective(reference, matrix_s_candidate, &
275 : matrix_s_candidate_reference, virtual_cutoff, &
276 : virtual_smoothing, gap_scale, objective_result)
277 : TYPE(frontier_orbitals_reference_type), INTENT(IN) :: reference
278 : TYPE(dbcsr_type), POINTER :: matrix_s_candidate, &
279 : matrix_s_candidate_reference
280 : REAL(KIND=dp), INTENT(IN) :: virtual_cutoff, virtual_smoothing, &
281 : gap_scale
282 : TYPE(frontier_orbitals_result_type), INTENT(OUT) :: objective_result
283 :
284 : INTEGER :: istate, naux, nref
285 408 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_candidate
286 408 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orbital_overlap_dense
287 : REAL(KIND=dp), DIMENSION(2) :: condnum
288 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
289 : TYPE(cp_fm_struct_type), POINTER :: ao_cross_struct, aux_struct, &
290 : cross_struct, ref_struct
291 : TYPE(cp_fm_type) :: candidate_eigenvectors, candidate_hamiltonian, candidate_overlap, &
292 : candidate_work, reference_candidate_ao_overlap, reference_candidate_orbital_overlap, &
293 : reference_eigenvalue_matrix, reference_orbital_candidate_ao_overlap
294 408 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: smat
295 :
296 408 : NULLIFY (ao_cross_struct, aux_struct, blacs_env, cross_struct, ref_struct, smat)
297 408 : nref = reference%number_reference_orbitals
298 408 : CALL dbcsr_get_info(matrix_s_candidate, nfullrows_total=naux)
299 408 : IF (naux <= reference%number_occupied_orbitals) THEN
300 0 : CPABORT("FRONTIER_ORBITALS candidate basis has no virtual orbitals")
301 : END IF
302 :
303 408 : blacs_env => reference%eigenvectors%matrix_struct%context
304 1224 : ALLOCATE (smat(1, 1))
305 408 : smat(1, 1)%matrix => matrix_s_candidate
306 408 : CALL overlap_condnum(smat, condnum, -1, .FALSE., .TRUE., .FALSE., blacs_env)
307 408 : DEALLOCATE (smat)
308 408 : objective_result%condition_number = condnum(2)
309 408 : objective_result%number_candidate_orbitals = naux
310 : IF (.NOT. ieee_is_finite(objective_result%condition_number) .OR. &
311 408 : objective_result%condition_number > 1.0E14_dp .OR. &
312 : objective_result%condition_number <= 0.0_dp) THEN
313 0 : objective_result%loss_occupied_subspace = 1.0E6_dp
314 0 : objective_result%loss_virtual_subspace = 1.0E6_dp
315 0 : objective_result%loss_empty_subspace = 1.0E6_dp
316 0 : objective_result%loss_gap = 1.0E6_dp
317 : RETURN
318 : END IF
319 :
320 408 : CALL cp_fm_struct_create(aux_struct, context=blacs_env, nrow_global=naux, ncol_global=naux)
321 : CALL cp_fm_struct_create(ao_cross_struct, context=blacs_env, &
322 408 : nrow_global=naux, ncol_global=nref)
323 408 : CALL cp_fm_struct_create(cross_struct, context=blacs_env, nrow_global=nref, ncol_global=naux)
324 408 : CALL cp_fm_struct_create(ref_struct, context=blacs_env, nrow_global=nref, ncol_global=nref)
325 :
326 : CALL cp_fm_create(candidate_hamiltonian, aux_struct, &
327 408 : name="frontier_orbitals_candidate_hamiltonian")
328 : CALL cp_fm_create(candidate_overlap, aux_struct, &
329 408 : name="frontier_orbitals_candidate_overlap")
330 : CALL cp_fm_create(candidate_eigenvectors, aux_struct, &
331 408 : name="frontier_orbitals_candidate_eigenvectors")
332 408 : CALL cp_fm_create(candidate_work, aux_struct, name="frontier_orbitals_candidate_work")
333 : CALL cp_fm_create(reference_candidate_ao_overlap, matrix_struct=ao_cross_struct, &
334 408 : name="frontier_orbitals_reference_candidate_ao_overlap")
335 : CALL cp_fm_create(reference_orbital_candidate_ao_overlap, matrix_struct=cross_struct, &
336 408 : name="frontier_orbitals_reference_orbital_candidate_ao_overlap")
337 : CALL cp_fm_create(reference_candidate_orbital_overlap, matrix_struct=cross_struct, &
338 408 : name="frontier_orbitals_reference_candidate_orbital_overlap")
339 : CALL cp_fm_create(reference_eigenvalue_matrix, matrix_struct=ref_struct, &
340 408 : name="frontier_orbitals_reference_eigenvalues")
341 408 : CALL cp_fm_struct_release(ao_cross_struct)
342 408 : CALL cp_fm_struct_release(aux_struct)
343 408 : CALL cp_fm_struct_release(cross_struct)
344 408 : CALL cp_fm_struct_release(ref_struct)
345 :
346 : ! B = C_R^T S_RC. The supplied cross-overlap matrix is stored as S_CR.
347 408 : CALL copy_dbcsr_to_fm(matrix_s_candidate_reference, reference_candidate_ao_overlap)
348 : CALL parallel_gemm('T', 'T', nref, naux, nref, 1.0_dp, reference%eigenvectors, &
349 : reference_candidate_ao_overlap, 0.0_dp, &
350 408 : reference_orbital_candidate_ao_overlap)
351 :
352 : ! H_C = B^T diag(ε_R) B.
353 408 : CALL cp_fm_set_all(reference_eigenvalue_matrix, 0.0_dp)
354 19584 : DO istate = 1, nref
355 : CALL cp_fm_set_element(reference_eigenvalue_matrix, istate, istate, &
356 19584 : reference%eigenvalues(istate))
357 : END DO
358 : CALL parallel_gemm('N', 'N', nref, naux, nref, 1.0_dp, reference_eigenvalue_matrix, &
359 : reference_orbital_candidate_ao_overlap, 0.0_dp, &
360 408 : reference_candidate_orbital_overlap)
361 : CALL parallel_gemm('T', 'N', naux, naux, nref, 1.0_dp, &
362 : reference_orbital_candidate_ao_overlap, &
363 408 : reference_candidate_orbital_overlap, 0.0_dp, candidate_hamiltonian)
364 :
365 408 : CALL copy_dbcsr_to_fm(matrix_s_candidate, candidate_overlap)
366 408 : CALL cp_fm_uplo_to_full(candidate_overlap, candidate_work)
367 1224 : ALLOCATE (eigenvalues_candidate(naux))
368 :
369 : ! H_C C_C = S_C C_C diag(ε_C).
370 : CALL cp_fm_geeig(candidate_hamiltonian, candidate_overlap, candidate_eigenvectors, &
371 408 : eigenvalues_candidate, candidate_work)
372 :
373 : ! M = B C_C, where M_ij = <ψ_i^R | ψ_j^C>.
374 : CALL parallel_gemm('N', 'N', nref, naux, naux, 1.0_dp, &
375 : reference_orbital_candidate_ao_overlap, candidate_eigenvectors, &
376 408 : 0.0_dp, reference_candidate_orbital_overlap)
377 1632 : ALLOCATE (orbital_overlap_dense(nref, naux))
378 408 : CALL cp_fm_get_submatrix(reference_candidate_orbital_overlap, orbital_overlap_dense)
379 : CALL compute_frontier_orbital_loss_function( &
380 : orbital_overlap_dense, reference%eigenvalues, eigenvalues_candidate, &
381 408 : reference%number_occupied_orbitals, virtual_cutoff, virtual_smoothing, gap_scale, objective_result)
382 408 : objective_result%condition_number = condnum(2)
383 408 : objective_result%number_candidate_orbitals = naux
384 :
385 408 : DEALLOCATE (eigenvalues_candidate, orbital_overlap_dense)
386 408 : CALL cp_fm_release(candidate_eigenvectors)
387 408 : CALL cp_fm_release(candidate_hamiltonian)
388 408 : CALL cp_fm_release(candidate_overlap)
389 408 : CALL cp_fm_release(candidate_work)
390 408 : CALL cp_fm_release(reference_candidate_ao_overlap)
391 408 : CALL cp_fm_release(reference_candidate_orbital_overlap)
392 408 : CALL cp_fm_release(reference_eigenvalue_matrix)
393 408 : CALL cp_fm_release(reference_orbital_candidate_ao_overlap)
394 :
395 1224 : END SUBROUTINE evaluate_frontier_orbitals_objective
396 :
397 : ! **************************************************************************************************
398 : !> \brief Compute the frontier-orbital loss function.
399 : !> \param orbital_overlap reference/candidate orbital overlaps
400 : !> \param eigenvalues_reference reference eigenvalues
401 : !> \param eigenvalues_candidate candidate eigenvalues
402 : !> \param number_occupied_orbitals number of occupied spatial orbitals
403 : !> \param virtual_cutoff virtual-orbital energy cutoff E_cut
404 : !> \param virtual_smoothing virtual-orbital energy smoothing width ΔE
405 : !> \param gap_scale normalization energy E_scale for the gap loss
406 : !> \param objective_result loss-function terms and result information
407 : !> \note Let i,j denote occupied orbitals, a,b virtual orbitals, and N_occ the number of occupied
408 : !> orbitals. R and C label reference and candidate quantities, with
409 : !> M_ij = <ψ_i^R | ψ_j^C> and M_ab = <ψ_a^R | ψ_b^C>. The occupied-subspace loss is
410 : !> L_occ = 1 - sum_(i,j in occupied) |M_ij|^2 / N_occ.
411 : !>
412 : !> For virtual orbitals, define the smooth energy-window weights and their norms as
413 : !> w_a^X = 1 / {1 + exp[(ε_a^X - ε_LUMO^X - E_cut) / ΔE]}, for X = R or C,
414 : !> N_X = sum_(a in virtual) (w_a^X)^2.
415 : !> The virtual-subspace loss is the normalized squared Frobenius distance between the
416 : !> energy-weighted virtual-space operators. Its expanded form is
417 : !> L_vir = [N_R + N_C - 2 sum_(a,b in virtual) w_a^R w_b^C |M_ab|^2] / (2 N_R).
418 : !>
419 : !> The empty-subspace loss is
420 : !> L_empty = 1 - sum_(a,b in virtual) (w_a^R)^2 |M_ab|^2 / N_R.
421 : !>
422 : !> With E_gap^X = ε_LUMO^X - ε_HOMO^X, the gap loss is
423 : !> L_gap = [(E_gap^C - E_gap^R) / E_scale]^2.
424 : ! **************************************************************************************************
425 816 : PURE SUBROUTINE compute_frontier_orbital_loss_function( &
426 408 : orbital_overlap, eigenvalues_reference, eigenvalues_candidate, number_occupied_orbitals, &
427 : virtual_cutoff, virtual_smoothing, gap_scale, objective_result)
428 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: orbital_overlap
429 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues_reference, &
430 : eigenvalues_candidate
431 : INTEGER, INTENT(IN) :: number_occupied_orbitals
432 : REAL(KIND=dp), INTENT(IN) :: virtual_cutoff, virtual_smoothing, &
433 : gap_scale
434 : TYPE(frontier_orbitals_result_type), INTENT(OUT) :: objective_result
435 :
436 : INTEGER :: candidate_orbital, reference_orbital
437 : REAL(KIND=dp) :: candidate_virtual_weight_norm, occupied_subspace_overlap, &
438 : reference_virtual_subspace_overlap, reference_virtual_weight_norm, &
439 : virtual_subspace_overlap, weight_candidate, weight_reference
440 :
441 : ! O_occ = sum_{i,j in occupied} |M_ij|^2 / N_occ.
442 : occupied_subspace_overlap = &
443 : SUM(orbital_overlap(1:number_occupied_orbitals, &
444 : 1:number_occupied_orbitals)**2)/ &
445 8568 : REAL(number_occupied_orbitals, KIND=dp)
446 408 : objective_result%loss_occupied_subspace = MAX(0.0_dp, 1.0_dp - occupied_subspace_overlap)
447 408 : objective_result%occupied_subspace_overlap = occupied_subspace_overlap
448 :
449 : ! N_R = sum_{a in virtual} (w_a^R)^2.
450 408 : reference_virtual_weight_norm = 0.0_dp
451 17952 : DO reference_orbital = number_occupied_orbitals + 1, SIZE(eigenvalues_reference)
452 : weight_reference = smooth_frontier_orbital_weight( &
453 : eigenvalues_reference(reference_orbital) - &
454 : eigenvalues_reference(number_occupied_orbitals + 1), &
455 17544 : virtual_cutoff, virtual_smoothing)
456 17952 : reference_virtual_weight_norm = reference_virtual_weight_norm + weight_reference**2
457 : END DO
458 :
459 : ! N_C = sum_{b in virtual} (w_b^C)^2.
460 408 : candidate_virtual_weight_norm = 0.0_dp
461 8976 : DO candidate_orbital = number_occupied_orbitals + 1, SIZE(eigenvalues_candidate)
462 : weight_candidate = smooth_frontier_orbital_weight( &
463 : eigenvalues_candidate(candidate_orbital) - &
464 : eigenvalues_candidate(number_occupied_orbitals + 1), &
465 8568 : virtual_cutoff, virtual_smoothing)
466 8976 : candidate_virtual_weight_norm = candidate_virtual_weight_norm + weight_candidate**2
467 : END DO
468 :
469 : ! O_vir = sum_{a,b in virtual} w_a^R w_b^C |M_ab|^2.
470 8976 : virtual_subspace_overlap = 0.0_dp
471 8976 : reference_virtual_subspace_overlap = 0.0_dp
472 8976 : DO candidate_orbital = number_occupied_orbitals + 1, SIZE(eigenvalues_candidate)
473 : weight_candidate = smooth_frontier_orbital_weight( &
474 : eigenvalues_candidate(candidate_orbital) - &
475 : eigenvalues_candidate(number_occupied_orbitals + 1), &
476 8568 : virtual_cutoff, virtual_smoothing)
477 377400 : DO reference_orbital = number_occupied_orbitals + 1, SIZE(eigenvalues_reference)
478 : weight_reference = smooth_frontier_orbital_weight( &
479 : eigenvalues_reference(reference_orbital) - &
480 : eigenvalues_reference(number_occupied_orbitals + 1), &
481 368424 : virtual_cutoff, virtual_smoothing)
482 : virtual_subspace_overlap = virtual_subspace_overlap + &
483 : weight_reference*weight_candidate* &
484 368424 : orbital_overlap(reference_orbital, candidate_orbital)**2
485 : reference_virtual_subspace_overlap = reference_virtual_subspace_overlap + &
486 : weight_reference**2* &
487 376992 : orbital_overlap(reference_orbital, candidate_orbital)**2
488 : END DO
489 : END DO
490 :
491 408 : IF (reference_virtual_weight_norm > TINY(1.0_dp)) THEN
492 : ! L_vir = (N_R + N_C - 2 O_vir)/(2 N_R).
493 : objective_result%loss_virtual_subspace = &
494 : MAX(0.0_dp, (reference_virtual_weight_norm + candidate_virtual_weight_norm - &
495 408 : 2.0_dp*virtual_subspace_overlap)/(2.0_dp*reference_virtual_weight_norm))
496 : ! Bound the computed overlap against small roundoff excursions from its exact range [0, 1].
497 : objective_result%empty_subspace_overlap = &
498 408 : MIN(1.0_dp, MAX(0.0_dp, reference_virtual_subspace_overlap/reference_virtual_weight_norm))
499 408 : objective_result%loss_empty_subspace = 1.0_dp - objective_result%empty_subspace_overlap
500 408 : IF (candidate_virtual_weight_norm > TINY(1.0_dp)) THEN
501 : objective_result%virtual_subspace_similarity = &
502 : MIN(1.0_dp, MAX(0.0_dp, &
503 : virtual_subspace_overlap/ &
504 408 : SQRT(reference_virtual_weight_norm*candidate_virtual_weight_norm)))
505 : END IF
506 : END IF
507 :
508 : ! E_gap^X = ε_LUMO^X - ε_HOMO^X.
509 : objective_result%gap_reference = eigenvalues_reference(number_occupied_orbitals + 1) - &
510 408 : eigenvalues_reference(number_occupied_orbitals)
511 : objective_result%gap_candidate = eigenvalues_candidate(number_occupied_orbitals + 1) - &
512 408 : eigenvalues_candidate(number_occupied_orbitals)
513 :
514 : ! L_gap = [(E_gap^C - E_gap^R)/E_scale]^2.
515 : objective_result%loss_gap = &
516 408 : ((objective_result%gap_candidate - objective_result%gap_reference)/gap_scale)**2
517 :
518 408 : END SUBROUTINE compute_frontier_orbital_loss_function
519 :
520 : ! **************************************************************************************************
521 : !> \brief Return the smooth energy-window weight for a virtual orbital.
522 : !> \param energy_from_lumo orbital energy relative to the LUMO
523 : !> \param cutoff outer energy of the virtual-orbital window
524 : !> \param smoothing smoothing width at the outer boundary
525 : !> \return energy-window weight in the interval [0, 1]
526 : ! **************************************************************************************************
527 403104 : PURE FUNCTION smooth_frontier_orbital_weight(energy_from_lumo, cutoff, smoothing) RESULT(weight)
528 : REAL(KIND=dp), INTENT(IN) :: energy_from_lumo, cutoff, smoothing
529 : REAL(KIND=dp) :: weight
530 :
531 : REAL(KIND=dp) :: argument
532 :
533 403104 : argument = (energy_from_lumo - cutoff)/smoothing
534 403104 : IF (argument >= 40.0_dp) THEN
535 : weight = 0.0_dp
536 88640 : ELSE IF (argument <= -40.0_dp) THEN
537 : weight = 1.0_dp
538 : ELSE
539 88640 : weight = 1.0_dp/(1.0_dp + EXP(argument))
540 : END IF
541 :
542 403104 : END FUNCTION smooth_frontier_orbital_weight
543 :
544 0 : END MODULE optbas_frontier_orbitals_utils
|
