(git:5e7fe52)
Loading...
Searching...
No Matches
optbas_frontier_orbitals_utils.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE 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
84CONTAINS
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 SUBROUTINE print_frontier_orbital_basis_sets(unit_nr, element_symbols, reference_basis_names, &
94 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 cpassert(SIZE(element_symbols) == SIZE(reference_basis_names))
102 cpassert(SIZE(element_symbols) == SIZE(initial_basis_names))
103 WRITE (unit_nr, '(1X,A,T14,A,T51,A)') "Atom kind", "Ref. basis set", "Initial basis set"
104 DO ikind = 1, SIZE(element_symbols)
105 WRITE (unit_nr, '(1X,A,T14,A,T51,A)') trim(element_symbols(ikind)), &
106 trim(reference_basis_names(ikind)), trim(initial_basis_names(ikind))
107 END DO
108
109 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
130 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: smat
131 TYPE(dft_control_type), POINTER :: dft_control
132
133 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 matrix_ks=matrix_ks, matrix_s=matrix_s, nelectron_spin=nelectron_spin)
136
137 nspins = SIZE(matrix_ks)
138 IF (nspins /= 1 .OR. dft_control%multiplicity /= 1) THEN
139 cpabort("FRONTIER_ORBITALS currently supports only closed-shell calculations")
140 END IF
141 IF (do_kpoints) THEN
142 cpabort("FRONTIER_ORBITALS currently supports only Γ-point calculations")
143 END IF
144 IF (ASSOCIATED(qs_env%x_data)) THEN
145 cpabort("FRONTIER_ORBITALS does not yet support hybrid functionals")
146 END IF
147
148 CALL dbcsr_get_info(matrix_s(1)%matrix, &
149 nfullrows_total=reference%number_reference_orbitals)
150 ALLOCATE (smat(1, 1))
151 smat(1, 1)%matrix => matrix_s(1)%matrix
152 CALL overlap_condnum(smat, condnum, -1, .false., .true., .false., blacs_env)
153 reference%condition_number = condnum(2)
154 DEALLOCATE (smat)
155
156 CALL store_basis_metadata(qs_env, reference)
157
158 ! A restricted spatial orbital contains two electrons.
159 reference%number_occupied_orbitals = nelectron_spin(1)/2
160 IF (2*reference%number_occupied_orbitals /= nelectron_spin(1)) THEN
161 cpabort("FRONTIER_ORBITALS requires an even closed-shell electron count")
162 END IF
163 IF (reference%number_occupied_orbitals < 1 .OR. &
164 reference%number_occupied_orbitals >= reference%number_reference_orbitals) THEN
165 cpabort("FRONTIER_ORBITALS requires occupied and virtual reference orbitals")
166 END IF
167
168 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 ncol_global=reference%number_reference_orbitals)
172 CALL cp_fm_create(reference_hamiltonian, fm_struct, &
173 name="frontier_orbitals_reference_hamiltonian")
174 CALL cp_fm_create(reference_overlap, fm_struct, &
175 name="frontier_orbitals_reference_overlap")
176 CALL cp_fm_create(reference%eigenvectors, fm_struct, &
177 name="frontier_orbitals_reference_eigenvectors")
178 CALL cp_fm_create(work, fm_struct, name="frontier_orbitals_reference_work")
179 CALL cp_fm_struct_release(fm_struct)
180
181 CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, reference_hamiltonian)
182 CALL cp_fm_uplo_to_full(reference_hamiltonian, work)
183 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, reference_overlap)
184 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 reference%eigenvalues, work)
189
190 CALL cp_fm_release(reference_hamiltonian)
191 CALL cp_fm_release(reference_overlap)
192 CALL cp_fm_release(work)
193
194 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 SUBROUTINE frontier_orbitals_reference_release(reference)
201 TYPE(frontier_orbitals_reference_type), &
202 INTENT(INOUT) :: reference
203
204 IF (ALLOCATED(reference%eigenvalues)) DEALLOCATE (reference%eigenvalues)
205 IF (ALLOCATED(reference%element_symbol)) DEALLOCATE (reference%element_symbol)
206 IF (ALLOCATED(reference%reference_basis_name)) DEALLOCATE (reference%reference_basis_name)
207 IF (ALLOCATED(reference%candidate_basis_name)) DEALLOCATE (reference%candidate_basis_name)
208 IF (ASSOCIATED(reference%eigenvectors%matrix_struct)) THEN
209 CALL cp_fm_release(reference%eigenvectors)
210 END IF
211 reference%number_occupied_orbitals = 0
212 reference%number_reference_orbitals = 0
213 reference%number_kinds = 0
214 reference%condition_number = 0.0_dp
215
216 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 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
231
232 NULLIFY (basis_set, qs_kind_set)
233 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
234 reference%number_kinds = SIZE(qs_kind_set)
235 ALLOCATE (reference%element_symbol(reference%number_kinds))
236 ALLOCATE (reference%reference_basis_name(reference%number_kinds))
237 ALLOCATE (reference%candidate_basis_name(reference%number_kinds))
238 reference%reference_basis_name = ""
239 reference%candidate_basis_name = ""
240
241 DO ikind = 1, reference%number_kinds
242 CALL get_qs_kind(qs_kind_set(ikind), element_symbol=reference%element_symbol(ikind))
243 NULLIFY (basis_set)
244 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type="ORB")
245 IF (ASSOCIATED(basis_set)) THEN
246 CALL get_gto_basis_set(basis_set, name=reference%reference_basis_name(ikind))
247 END IF
248 NULLIFY (basis_set)
249 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, basis_type="AUX_OPT")
250 IF (ASSOCIATED(basis_set)) THEN
251 CALL get_gto_basis_set(basis_set, name=reference%candidate_basis_name(ikind))
252 END IF
253 END DO
254
255 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_candidate
286 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 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: smat
295
296 NULLIFY (ao_cross_struct, aux_struct, blacs_env, cross_struct, ref_struct, smat)
297 nref = reference%number_reference_orbitals
298 CALL dbcsr_get_info(matrix_s_candidate, nfullrows_total=naux)
299 IF (naux <= reference%number_occupied_orbitals) THEN
300 cpabort("FRONTIER_ORBITALS candidate basis has no virtual orbitals")
301 END IF
302
303 blacs_env => reference%eigenvectors%matrix_struct%context
304 ALLOCATE (smat(1, 1))
305 smat(1, 1)%matrix => matrix_s_candidate
306 CALL overlap_condnum(smat, condnum, -1, .false., .true., .false., blacs_env)
307 DEALLOCATE (smat)
308 objective_result%condition_number = condnum(2)
309 objective_result%number_candidate_orbitals = naux
310 IF (.NOT. ieee_is_finite(objective_result%condition_number) .OR. &
311 objective_result%condition_number > 1.0e14_dp .OR. &
312 objective_result%condition_number <= 0.0_dp) THEN
313 objective_result%loss_occupied_subspace = 1.0e6_dp
314 objective_result%loss_virtual_subspace = 1.0e6_dp
315 objective_result%loss_empty_subspace = 1.0e6_dp
316 objective_result%loss_gap = 1.0e6_dp
317 RETURN
318 END IF
319
320 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 nrow_global=naux, ncol_global=nref)
323 CALL cp_fm_struct_create(cross_struct, context=blacs_env, nrow_global=nref, ncol_global=naux)
324 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 name="frontier_orbitals_candidate_hamiltonian")
328 CALL cp_fm_create(candidate_overlap, aux_struct, &
329 name="frontier_orbitals_candidate_overlap")
330 CALL cp_fm_create(candidate_eigenvectors, aux_struct, &
331 name="frontier_orbitals_candidate_eigenvectors")
332 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 name="frontier_orbitals_reference_candidate_ao_overlap")
335 CALL cp_fm_create(reference_orbital_candidate_ao_overlap, matrix_struct=cross_struct, &
336 name="frontier_orbitals_reference_orbital_candidate_ao_overlap")
337 CALL cp_fm_create(reference_candidate_orbital_overlap, matrix_struct=cross_struct, &
338 name="frontier_orbitals_reference_candidate_orbital_overlap")
339 CALL cp_fm_create(reference_eigenvalue_matrix, matrix_struct=ref_struct, &
340 name="frontier_orbitals_reference_eigenvalues")
341 CALL cp_fm_struct_release(ao_cross_struct)
342 CALL cp_fm_struct_release(aux_struct)
343 CALL cp_fm_struct_release(cross_struct)
344 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 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 reference_orbital_candidate_ao_overlap)
351
352 ! H_C = B^T diag(ε_R) B.
353 CALL cp_fm_set_all(reference_eigenvalue_matrix, 0.0_dp)
354 DO istate = 1, nref
355 CALL cp_fm_set_element(reference_eigenvalue_matrix, istate, istate, &
356 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 reference_candidate_orbital_overlap)
361 CALL parallel_gemm('T', 'N', naux, naux, nref, 1.0_dp, &
362 reference_orbital_candidate_ao_overlap, &
363 reference_candidate_orbital_overlap, 0.0_dp, candidate_hamiltonian)
364
365 CALL copy_dbcsr_to_fm(matrix_s_candidate, candidate_overlap)
366 CALL cp_fm_uplo_to_full(candidate_overlap, candidate_work)
367 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 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 0.0_dp, reference_candidate_orbital_overlap)
377 ALLOCATE (orbital_overlap_dense(nref, naux))
378 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 reference%number_occupied_orbitals, virtual_cutoff, virtual_smoothing, gap_scale, objective_result)
382 objective_result%condition_number = condnum(2)
383 objective_result%number_candidate_orbitals = naux
384
385 DEALLOCATE (eigenvalues_candidate, orbital_overlap_dense)
386 CALL cp_fm_release(candidate_eigenvectors)
387 CALL cp_fm_release(candidate_hamiltonian)
388 CALL cp_fm_release(candidate_overlap)
389 CALL cp_fm_release(candidate_work)
390 CALL cp_fm_release(reference_candidate_ao_overlap)
391 CALL cp_fm_release(reference_candidate_orbital_overlap)
392 CALL cp_fm_release(reference_eigenvalue_matrix)
393 CALL cp_fm_release(reference_orbital_candidate_ao_overlap)
394
395 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 PURE SUBROUTINE compute_frontier_orbital_loss_function( &
426 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 REAL(number_occupied_orbitals, kind=dp)
446 objective_result%loss_occupied_subspace = max(0.0_dp, 1.0_dp - occupied_subspace_overlap)
447 objective_result%occupied_subspace_overlap = occupied_subspace_overlap
448
449 ! N_R = sum_{a in virtual} (w_a^R)^2.
450 reference_virtual_weight_norm = 0.0_dp
451 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 virtual_cutoff, virtual_smoothing)
456 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 candidate_virtual_weight_norm = 0.0_dp
461 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 virtual_cutoff, virtual_smoothing)
466 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 virtual_subspace_overlap = 0.0_dp
471 reference_virtual_subspace_overlap = 0.0_dp
472 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 virtual_cutoff, virtual_smoothing)
477 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 virtual_cutoff, virtual_smoothing)
482 virtual_subspace_overlap = virtual_subspace_overlap + &
483 weight_reference*weight_candidate* &
484 orbital_overlap(reference_orbital, candidate_orbital)**2
485 reference_virtual_subspace_overlap = reference_virtual_subspace_overlap + &
486 weight_reference**2* &
487 orbital_overlap(reference_orbital, candidate_orbital)**2
488 END DO
489 END DO
490
491 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 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 min(1.0_dp, max(0.0_dp, reference_virtual_subspace_overlap/reference_virtual_weight_norm))
499 objective_result%loss_empty_subspace = 1.0_dp - objective_result%empty_subspace_overlap
500 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 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 eigenvalues_reference(number_occupied_orbitals)
511 objective_result%gap_candidate = eigenvalues_candidate(number_occupied_orbitals + 1) - &
512 eigenvalues_candidate(number_occupied_orbitals)
513
514 ! L_gap = [(E_gap^C - E_gap^R)/E_scale]^2.
515 objective_result%loss_gap = &
516 ((objective_result%gap_candidate - objective_result%gap_reference)/gap_scale)**2
517
518 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 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 argument = (energy_from_lumo - cutoff)/smoothing
534 IF (argument >= 40.0_dp) THEN
535 weight = 0.0_dp
536 ELSE IF (argument <= -40.0_dp) THEN
537 weight = 1.0_dp
538 ELSE
539 weight = 1.0_dp/(1.0_dp + exp(argument))
540 END IF
541
542 END FUNCTION smooth_frontier_orbital_weight
543
544END MODULE optbas_frontier_orbitals_utils
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE. Use cuSOLVERMp directly when requested and large enough; otherwi...
Definition cp_fm_diag.F:1503
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition cp_fm_types.F:684
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
optbas_frontier_orbitals_utils
Utilities for frontier-orbital basis optimization.
Definition optbas_frontier_orbitals_utils.F:13
optbas_frontier_orbitals_utils::frontier_orbitals_reference_init
subroutine, public frontier_orbitals_reference_init(qs_env, reference)
Diagonalize and store the frozen reference Hamiltonian.
Definition optbas_frontier_orbitals_utils.F:117
optbas_frontier_orbitals_utils::frontier_orbitals_reference_release
subroutine, public frontier_orbitals_reference_release(reference)
Release the stored reference eigenvalues and eigenvectors.
Definition optbas_frontier_orbitals_utils.F:201
optbas_frontier_orbitals_utils::evaluate_frontier_orbitals_objective
subroutine, public evaluate_frontier_orbitals_objective(reference, matrix_s_candidate, matrix_s_candidate_reference, virtual_cutoff, virtual_smoothing, gap_scale, objective_result)
Construct and diagonalize the candidate Hamiltonian and evaluate the loss function.
Definition optbas_frontier_orbitals_utils.F:277
optbas_frontier_orbitals_utils::smooth_frontier_orbital_weight
pure real(kind=dp) function, public smooth_frontier_orbital_weight(energy_from_lumo, cutoff, smoothing)
Return the smooth energy-window weight for a virtual orbital.
Definition optbas_frontier_orbitals_utils.F:528
optbas_frontier_orbitals_utils::compute_frontier_orbital_loss_function
pure subroutine, public compute_frontier_orbital_loss_function(orbital_overlap, eigenvalues_reference, eigenvalues_candidate, number_occupied_orbitals, virtual_cutoff, virtual_smoothing, gap_scale, objective_result)
Compute the frontier-orbital loss function.
Definition optbas_frontier_orbitals_utils.F:428
optbas_frontier_orbitals_utils::print_frontier_orbital_basis_sets
subroutine, public print_frontier_orbital_basis_sets(unit_nr, element_symbols, reference_basis_names, initial_basis_names)
Print reference and initial basis-set names for every fitted atom kind.
Definition optbas_frontier_orbitals_utils.F:95
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
qs_condnum
Calculation of overlap matrix condition numbers.
Definition qs_condnum.F:13
qs_condnum::overlap_condnum
subroutine, public overlap_condnum(matrixkp_s, condnum, iunit, norml1, norml2, use_arnoldi, blacs_env)
Calculation of the overlap matrix Condition Number.
Definition qs_condnum.F:66
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:758
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
optbas_frontier_orbitals_utils::frontier_orbitals_reference_type
Definition optbas_frontier_orbitals_utils.F:66
optbas_frontier_orbitals_utils::frontier_orbitals_result_type
Definition optbas_frontier_orbitals_utils.F:49
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190