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 Automatic screening of frontier-orbital basis optimizations
9 : !> \par History
10 : !> 28.08.2026 created [Jan Wilhelm]
11 : !> \author Jan Wilhelm
12 : ! **************************************************************************************************
13 : MODULE optbas_frontier_orbital_screening
14 : USE cp_blacs_env, ONLY: cp_blacs_env_type
15 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
16 : USE cp_files, ONLY: close_file,&
17 : discover_file,&
18 : open_file
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_get_default_unit_nr,&
21 : cp_logger_type,&
22 : cp_to_string
23 : USE cp_output_handling, ONLY: debug_print_level,&
24 : silent_print_level
25 : USE global_types, ONLY: global_environment_type
26 : USE input_constants, ONLY: atomic_guess,&
27 : do_qs
28 : USE input_section_types, ONLY: section_type,&
29 : section_vals_duplicate,&
30 : section_vals_get,&
31 : section_vals_get_subs_vals,&
32 : section_vals_release,&
33 : section_vals_type,&
34 : section_vals_val_get,&
35 : section_vals_val_set
36 : USE kinds, ONLY: default_path_length,&
37 : default_string_length,&
38 : dp
39 : USE machine, ONLY: m_getcwd
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE optbas_fenv_manipulation, ONLY: update_basis_set
42 : USE optbas_frontier_orbitals_utils, ONLY: frontier_orbitals_reference_init,&
43 : frontier_orbitals_reference_release,&
44 : frontier_orbitals_reference_type,&
45 : print_frontier_orbital_basis_sets
46 : USE optimize_basis_types, ONLY: basis_optimization_type,&
47 : deallocate_basis_optimization_type,&
48 : method_mo_fit_occ_virtual
49 : USE optimize_basis_utils, ONLY: optimize_basis_init_read_input,&
50 : write_basis
51 : USE physcon, ONLY: evolt
52 : USE qs_condnum, ONLY: overlap_condnum
53 : USE qs_energy, ONLY: qs_energies
54 : USE qs_energy_types, ONLY: qs_energy_type
55 : USE qs_environment, ONLY: qs_init
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_env_create,&
58 : qs_env_release,&
59 : qs_environment_type
60 : USE qs_mo_types, ONLY: get_mo_set,&
61 : mo_set_type
62 : USE string_utilities, ONLY: uppercase
63 : #include "./base/base_uses.f90"
64 :
65 : IMPLICIT NONE
66 : PRIVATE
67 :
68 : PUBLIC :: run_frontier_orbital_screening
69 :
70 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optbas_frontier_orbital_screening'
71 :
72 : TYPE frontier_orbital_screening_candidate_type
73 : REAL(KIND=dp) :: condition_number = 0.0_dp
74 : REAL(KIND=dp) :: gap = 0.0_dp
75 : REAL(KIND=dp) :: energy = 0.0_dp
76 : TYPE(basis_optimization_type) :: opt_bas
77 : END TYPE frontier_orbital_screening_candidate_type
78 :
79 : ABSTRACT INTERFACE
80 : ! **************************************************************************************************
81 : !> \brief ...
82 : !> \param opt_bas ...
83 : !> \param input_declaration ...
84 : !> \param para_env ...
85 : !> \param globenv ...
86 : !> \param training_input ...
87 : !> \param reference ...
88 : ! **************************************************************************************************
89 : SUBROUTINE basis_optimization_driver(opt_bas, input_declaration, para_env, globenv, training_input, reference)
90 : IMPORT :: basis_optimization_type, frontier_orbitals_reference_type, global_environment_type, &
91 : mp_para_env_type, &
92 : section_type, section_vals_type
93 : TYPE(basis_optimization_type) :: opt_bas
94 : TYPE(section_type), POINTER :: input_declaration
95 : TYPE(mp_para_env_type), POINTER :: para_env
96 : TYPE(global_environment_type), POINTER :: globenv
97 : TYPE(section_vals_type), OPTIONAL, POINTER :: training_input
98 : TYPE(frontier_orbitals_reference_type), &
99 : DIMENSION(:), INTENT(IN), OPTIONAL, TARGET :: reference
100 : END SUBROUTINE basis_optimization_driver
101 : END INTERFACE
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief Run serial frontier-orbital optimizations and evaluate their HOMO-LUMO gaps with new SCFs.
107 : !> \param input_declaration ...
108 : !> \param section ...
109 : !> \param para_env ...
110 : !> \param globenv ...
111 : !> \param frontier_screening_section ...
112 : !> \param optimize_basis_driver ...
113 : ! **************************************************************************************************
114 2 : SUBROUTINE run_frontier_orbital_screening(input_declaration, section, para_env, globenv, frontier_screening_section, &
115 : optimize_basis_driver)
116 : TYPE(section_type), POINTER :: input_declaration
117 : TYPE(section_vals_type), POINTER :: section
118 : TYPE(mp_para_env_type), POINTER :: para_env
119 : TYPE(global_environment_type), POINTER :: globenv
120 : TYPE(section_vals_type), POINTER :: frontier_screening_section
121 :
122 : PROCEDURE(basis_optimization_driver) :: optimize_basis_driver
123 :
124 : INTEGER, PARAMETER :: number_parameter_sets = 5
125 : REAL(KIND=dp), DIMENSION(6, number_parameter_sets), PARAMETER :: &
126 : pre_tabulated_parameters = RESHAPE( &
127 : [10.0_dp, 10.0_dp, 10.0_dp, 10.0_dp, 0.001_dp, 0.001_dp, &
128 : 30.0_dp, 30.0_dp, 30.0_dp, 10.0_dp, 0.001_dp, 0.001_dp, &
129 : 10.0_dp, 30.0_dp, 30.0_dp, 10.0_dp, 0.001_dp, 0.001_dp, &
130 : 10.0_dp, 30.0_dp, 30.0_dp, 3.0_dp, 0.001_dp, 0.001_dp, &
131 : 10.0_dp, 30.0_dp, 30.0_dp, 10.0_dp, 0.001_dp, 0.01_dp], &
132 : [6, number_parameter_sets])
133 :
134 : CHARACTER(LEN=default_path_length) :: final_basis_file, main_dir
135 : INTEGER :: i, n_candidates, unit_nr
136 : REAL(KIND=dp) :: initial_condition_number, initial_energy, initial_gap, &
137 : reference_condition_number, reference_energy, reference_gap
138 : TYPE(cp_logger_type), POINTER :: logger
139 : TYPE(frontier_orbitals_reference_type), &
140 2 : ALLOCATABLE, DIMENSION(:), TARGET :: reference
141 : TYPE(frontier_orbital_screening_candidate_type), ALLOCATABLE, &
142 2 : DIMENSION(:) :: candidates
143 :
144 : ! Rows are occupied, virtual, empty, gap, condition, and coefficient weights.
145 :
146 2 : CALL check_frontier_orbital_screening_input(section)
147 2 : CALL section_vals_val_get(frontier_screening_section, "NUMBER_OF_OPTIMIZATIONS", i_val=n_candidates)
148 2 : IF (n_candidates < 1 .OR. n_candidates > number_parameter_sets) THEN
149 : CALL cp_abort(__LOCATION__, &
150 : "NUMBER_OF_OPTIMIZATIONS must be between 1 and "// &
151 0 : TRIM(cp_to_string(number_parameter_sets)))
152 : END IF
153 :
154 2 : CALL m_getcwd(main_dir)
155 10 : ALLOCATE (candidates(n_candidates), reference(1))
156 2 : logger => cp_get_default_logger()
157 2 : unit_nr = -1
158 2 : IF (para_env%is_source()) unit_nr = cp_logger_get_default_unit_nr(logger)
159 :
160 : CALL run_frontier_orbital_reference_scf(section, para_env, globenv, reference(1), &
161 : reference_gap, reference_energy, &
162 2 : reference_condition_number)
163 :
164 4 : DO i = 1, n_candidates
165 : CALL optimize_basis_init_read_input(candidates(i)%opt_bas, section, para_env, &
166 2 : quiet_output=.TRUE., embedded_training=.TRUE.)
167 2 : IF (candidates(i)%opt_bas%method /= method_mo_fit_occ_virtual) THEN
168 0 : CPABORT("Internal error while initializing FRONTIER_ORBITAL_SCREENING")
169 : END IF
170 2 : IF (candidates(i)%opt_bas%ntraining_sets /= 1 .OR. candidates(i)%opt_bas%ncombinations /= 1) THEN
171 : CALL cp_abort(__LOCATION__, &
172 : "FRONTIER_ORBITAL_SCREENING currently supports exactly one training system "// &
173 0 : "and one candidate basis assignment")
174 : END IF
175 :
176 2 : IF (i == 1) THEN
177 2 : final_basis_file = work_path(main_dir, candidates(i)%opt_bas%output_basis_file)
178 2 : IF (unit_nr > 0) THEN
179 1 : WRITE (unit_nr, '(A)') ""
180 1 : WRITE (unit_nr, '(1X,A)') REPEAT("=", 79)
181 1 : WRITE (unit_nr, '(1X,A)') "Optimize basis set to improve frontier orbitals:"
182 1 : WRITE (unit_nr, '(1X,A)') "Using pre-tabulated parameters in the optimization"
183 1 : WRITE (unit_nr, '(1X,A)') REPEAT("=", 79)
184 1 : WRITE (unit_nr, '(A)') " "
185 : END IF
186 : END IF
187 :
188 2 : candidates(i)%opt_bas%method = method_mo_fit_occ_virtual
189 2 : candidates(i)%opt_bas%use_condition_number = .TRUE.
190 2 : candidates(i)%opt_bas%write_frequency = 5000
191 6 : candidates(i)%opt_bas%fval_weight = 1.0_dp
192 2 : candidates(i)%opt_bas%occupied_weight = pre_tabulated_parameters(1, i)
193 2 : candidates(i)%opt_bas%virtual_weight = pre_tabulated_parameters(2, i)
194 2 : candidates(i)%opt_bas%empty_overlap_weight = pre_tabulated_parameters(3, i)
195 2 : candidates(i)%opt_bas%gap_weight = pre_tabulated_parameters(4, i)
196 6 : candidates(i)%opt_bas%condition_weight = pre_tabulated_parameters(5, i)
197 2 : candidates(i)%opt_bas%coefficient_weight = pre_tabulated_parameters(6, i)
198 2 : candidates(i)%opt_bas%gap_energy_scale = 1.0_dp/evolt
199 2 : candidates(i)%opt_bas%virtual_energy_cutoff = 3.0_dp/evolt
200 2 : candidates(i)%opt_bas%virtual_energy_smoothing = 0.2_dp/evolt
201 2 : IF (unit_nr > 0) THEN
202 1 : WRITE (unit_nr, '(1X,A,I0,A,I0,A)') "Candidate ", i, " of ", n_candidates, &
203 2 : ": optimizing basis set ..."
204 : END IF
205 : CALL optimize_basis_driver(candidates(i)%opt_bas, input_declaration, para_env, globenv, &
206 4 : training_input=section, reference=reference)
207 : END DO
208 :
209 : CALL run_frontier_orbital_scf(candidates(1)%opt_bas, section, para_env, globenv, &
210 : optimized_basis=.FALSE., &
211 : gap=initial_gap, energy=initial_energy, &
212 2 : condition_number=initial_condition_number)
213 :
214 4 : DO i = 1, n_candidates
215 : CALL run_frontier_orbital_scf(candidates(i)%opt_bas, section, para_env, globenv, &
216 : optimized_basis=.TRUE., &
217 : gap=candidates(i)%gap, energy=candidates(i)%energy, &
218 4 : condition_number=candidates(i)%condition_number)
219 : END DO
220 :
221 : CALL write_frontier_orbital_basis_file(final_basis_file, candidates, &
222 2 : reference_gap, reference_energy, reference_condition_number, para_env)
223 : CALL print_frontier_orbital_summary(unit_nr, final_basis_file, section, candidates, &
224 : reference_gap, reference_energy, reference_condition_number, &
225 2 : initial_gap, initial_energy, initial_condition_number)
226 :
227 4 : DO i = 1, n_candidates
228 4 : CALL deallocate_basis_optimization_type(candidates(i)%opt_bas)
229 : END DO
230 2 : CALL frontier_orbitals_reference_release(reference(1))
231 6 : DEALLOCATE (candidates, reference)
232 :
233 4 : END SUBROUTINE run_frontier_orbital_screening
234 :
235 : ! **************************************************************************************************
236 : !> \brief Check the single-file input contract of frontier-orbital screening.
237 : !> \param section ...
238 : ! **************************************************************************************************
239 12 : SUBROUTINE check_frontier_orbital_screening_input(section)
240 : TYPE(section_vals_type), POINTER :: section
241 :
242 : INTEGER :: method, nforce_eval, ntraining
243 : LOGICAL :: dft_explicit, keyword_explicit
244 : TYPE(section_vals_type), POINTER :: dft_section, force_eval_section, &
245 : optbas_section, training_section
246 :
247 2 : optbas_section => section_vals_get_subs_vals(section, "OPTIMIZE_BASIS")
248 :
249 2 : CALL reject_controlled_keyword(optbas_section, "BASIS_WORK_FILE")
250 2 : CALL reject_controlled_keyword(optbas_section, "WRITE_FREQUENCY")
251 2 : CALL reject_controlled_keyword(optbas_section, "USE_CONDITION_NUMBER")
252 2 : CALL reject_controlled_keyword(optbas_section, "RESIDUUM_WEIGHT")
253 2 : CALL reject_controlled_keyword(optbas_section, "CONDITION_WEIGHT")
254 2 : CALL reject_controlled_keyword(optbas_section, "BASIS_COMBINATIONS")
255 2 : CALL reject_controlled_keyword(optbas_section, "GROUP_PARTITION")
256 2 : CALL section_vals_val_get(optbas_section, "BASIS_TEMPLATE_FILE", explicit=keyword_explicit)
257 2 : IF (.NOT. keyword_explicit) THEN
258 0 : CPABORT("FRONTIER_ORBITAL_SCREENING requires BASIS_TEMPLATE_FILE")
259 : END IF
260 :
261 2 : training_section => section_vals_get_subs_vals(optbas_section, "TRAINING_FILES")
262 2 : CALL section_vals_get(training_section, n_repetition=ntraining)
263 2 : IF (ntraining /= 0) CALL controlled_parameter_error("TRAINING_FILES")
264 :
265 2 : force_eval_section => section_vals_get_subs_vals(section, "FORCE_EVAL")
266 2 : CALL section_vals_get(force_eval_section, n_repetition=nforce_eval)
267 2 : IF (nforce_eval /= 1) THEN
268 : CALL cp_abort(__LOCATION__, &
269 0 : "FRONTIER_ORBITAL_SCREENING requires exactly one FORCE_EVAL in the same input")
270 : END IF
271 2 : CALL section_vals_val_get(force_eval_section, "METHOD", i_val=method)
272 2 : IF (method /= do_qs) THEN
273 0 : CPABORT("FRONTIER_ORBITAL_SCREENING requires a Quickstep FORCE_EVAL")
274 : END IF
275 2 : dft_section => section_vals_get_subs_vals(force_eval_section, "DFT")
276 2 : CALL section_vals_get(dft_section, explicit=dft_explicit)
277 2 : IF (.NOT. dft_explicit) THEN
278 : CALL cp_abort(__LOCATION__, &
279 0 : "FRONTIER_ORBITAL_SCREENING requires an explicit FORCE_EVAL%DFT section")
280 : END IF
281 2 : END SUBROUTINE check_frontier_orbital_screening_input
282 :
283 : ! **************************************************************************************************
284 : !> \brief Reject an OPTIMIZE_BASIS keyword controlled by frontier-orbital screening.
285 : !> \param section ...
286 : !> \param keyword_name ...
287 : ! **************************************************************************************************
288 14 : SUBROUTINE reject_controlled_keyword(section, keyword_name)
289 : TYPE(section_vals_type), POINTER :: section
290 : CHARACTER(LEN=*), INTENT(IN) :: keyword_name
291 :
292 : LOGICAL :: explicit
293 :
294 14 : CALL section_vals_val_get(section, keyword_name, explicit=explicit)
295 14 : IF (explicit) CALL controlled_parameter_error("OPTIMIZE_BASIS%"//TRIM(keyword_name))
296 :
297 14 : END SUBROUTINE reject_controlled_keyword
298 :
299 : ! **************************************************************************************************
300 : !> \brief Abort for an input parameter replaced by frontier-orbital screening.
301 : !> \param parameter_name ...
302 : ! **************************************************************************************************
303 0 : SUBROUTINE controlled_parameter_error(parameter_name)
304 : CHARACTER(LEN=*), INTENT(IN) :: parameter_name
305 :
306 : CALL cp_abort(__LOCATION__, &
307 : TRIM(parameter_name)//" is not allowed together with FRONTIER_ORBITAL_SCREENING "// &
308 0 : "because it is overwritten automatically")
309 :
310 0 : END SUBROUTINE controlled_parameter_error
311 :
312 : ! **************************************************************************************************
313 : !> \brief Run and retain the large-basis reference from the embedded FORCE_EVAL.
314 : !> \param section ...
315 : !> \param para_env ...
316 : !> \param globenv ...
317 : !> \param reference ...
318 : !> \param gap ...
319 : !> \param energy ...
320 : !> \param condition_number ...
321 : ! **************************************************************************************************
322 2 : SUBROUTINE run_frontier_orbital_reference_scf(section, para_env, globenv, reference, &
323 : gap, energy, condition_number)
324 : TYPE(section_vals_type), POINTER :: section
325 : TYPE(mp_para_env_type), POINTER :: para_env
326 : TYPE(global_environment_type), POINTER :: globenv
327 : TYPE(frontier_orbitals_reference_type), &
328 : INTENT(OUT) :: reference
329 : REAL(KIND=dp), INTENT(OUT) :: gap, energy, condition_number
330 :
331 : INTEGER :: outer_print_level
332 : TYPE(cp_logger_type), POINTER :: logger
333 : TYPE(qs_environment_type), POINTER :: qs_env
334 : TYPE(section_vals_type), POINTER :: input_file
335 :
336 2 : NULLIFY (input_file, logger, qs_env)
337 2 : CALL section_vals_duplicate(section, input_file)
338 2 : CALL section_vals_val_set(input_file, "GLOBAL%PRINT_LEVEL", i_val=silent_print_level)
339 2 : CALL section_vals_val_set(input_file, "FORCE_EVAL%DFT%SCF%SCF_GUESS", i_val=atomic_guess)
340 : CALL section_vals_val_set(input_file, "GLOBAL%PROGRAM_RUN_INFO%_SECTION_PARAMETERS_", &
341 2 : i_val=debug_print_level)
342 : CALL section_vals_val_set(input_file, "GLOBAL%REFERENCES%_SECTION_PARAMETERS_", &
343 2 : i_val=debug_print_level)
344 : CALL section_vals_val_set(input_file, "GLOBAL%TIMINGS%_SECTION_PARAMETERS_", &
345 2 : i_val=debug_print_level)
346 : CALL section_vals_val_set(input_file, &
347 : "FORCE_EVAL%DFT%SCF%PRINT%RESTART%_SECTION_PARAMETERS_", &
348 2 : i_val=debug_print_level + 1)
349 : CALL section_vals_val_set(input_file, &
350 : "FORCE_EVAL%DFT%SCF%PRINT%RESTART_HISTORY%_SECTION_PARAMETERS_", &
351 2 : i_val=debug_print_level + 1)
352 2 : logger => cp_get_default_logger()
353 2 : outer_print_level = logger%iter_info%print_level
354 2 : logger%iter_info%print_level = silent_print_level
355 2 : CALL create_silent_qs_env(input_file, para_env, globenv, qs_env)
356 2 : CALL qs_energies(qs_env)
357 2 : CALL get_frontier_orbital_scf_result(qs_env, gap, energy, condition_number)
358 2 : CALL frontier_orbitals_reference_init(qs_env, reference)
359 2 : CALL qs_env_release(qs_env)
360 2 : DEALLOCATE (qs_env)
361 2 : logger%iter_info%print_level = outer_print_level
362 2 : CALL section_vals_release(input_file)
363 :
364 2 : END SUBROUTINE run_frontier_orbital_reference_scf
365 :
366 : ! **************************************************************************************************
367 : !> \brief Run one silent SCF using a compact basis held in memory.
368 : !> \param opt_bas ...
369 : !> \param section ...
370 : !> \param para_env ...
371 : !> \param globenv ...
372 : !> \param optimized_basis ...
373 : !> \param gap ...
374 : !> \param energy ...
375 : !> \param condition_number ...
376 : ! **************************************************************************************************
377 8 : SUBROUTINE run_frontier_orbital_scf(opt_bas, section, para_env, globenv, &
378 : optimized_basis, gap, energy, condition_number)
379 : TYPE(basis_optimization_type), INTENT(IN) :: opt_bas
380 : TYPE(section_vals_type), POINTER :: section
381 : TYPE(mp_para_env_type), POINTER :: para_env
382 : TYPE(global_environment_type), POINTER :: globenv
383 : LOGICAL, INTENT(IN) :: optimized_basis
384 : REAL(KIND=dp), INTENT(OUT) :: gap, energy, condition_number
385 :
386 : INTEGER :: outer_print_level
387 : TYPE(cp_logger_type), POINTER :: logger
388 : TYPE(qs_environment_type), POINTER :: qs_env
389 : TYPE(section_vals_type), POINTER :: input_file
390 :
391 4 : NULLIFY (input_file, logger, qs_env)
392 4 : gap = 0.0_dp
393 4 : energy = 0.0_dp
394 4 : condition_number = 0.0_dp
395 4 : CALL section_vals_duplicate(section, input_file)
396 4 : CALL prepare_frontier_orbital_scf_input(input_file, opt_bas)
397 4 : logger => cp_get_default_logger()
398 4 : outer_print_level = logger%iter_info%print_level
399 4 : logger%iter_info%print_level = silent_print_level
400 4 : CALL create_silent_qs_env(input_file, para_env, globenv, qs_env)
401 4 : IF (optimized_basis) CALL update_basis_set(opt_bas, 1, "ORB", qs_env)
402 4 : CALL qs_energies(qs_env)
403 4 : CALL get_frontier_orbital_scf_result(qs_env, gap, energy, condition_number)
404 4 : CALL qs_env_release(qs_env)
405 4 : DEALLOCATE (qs_env)
406 4 : logger%iter_info%print_level = outer_print_level
407 4 : CALL section_vals_release(input_file)
408 :
409 4 : END SUBROUTINE run_frontier_orbital_scf
410 :
411 : ! **************************************************************************************************
412 : !> \brief Initialize Quickstep directly under the active CP2K environment without nested startup.
413 : !> \param input_file ...
414 : !> \param para_env ...
415 : !> \param globenv ...
416 : !> \param qs_env ...
417 : ! **************************************************************************************************
418 12 : SUBROUTINE create_silent_qs_env(input_file, para_env, globenv, qs_env)
419 : TYPE(section_vals_type), POINTER :: input_file
420 : TYPE(mp_para_env_type), POINTER :: para_env
421 : TYPE(global_environment_type), POINTER :: globenv
422 : TYPE(qs_environment_type), POINTER :: qs_env
423 :
424 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
425 :
426 6 : force_env_section => section_vals_get_subs_vals(input_file, "FORCE_EVAL")
427 6 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
428 6 : ALLOCATE (qs_env)
429 6 : CALL qs_env_create(qs_env, globenv)
430 : CALL qs_init(qs_env, para_env, input_file, globenv=globenv, &
431 : force_env_section=force_env_section, subsys_section=subsys_section, &
432 6 : use_motion_section=.TRUE., silent=.TRUE.)
433 :
434 6 : END SUBROUTINE create_silent_qs_env
435 :
436 : ! **************************************************************************************************
437 : !> \brief Extract the gap, energy, and overlap condition number from a completed SCF.
438 : !> \param qs_env ...
439 : !> \param gap ...
440 : !> \param energy ...
441 : !> \param condition_number ...
442 : ! **************************************************************************************************
443 6 : SUBROUTINE get_frontier_orbital_scf_result(qs_env, gap, energy, condition_number)
444 : TYPE(qs_environment_type), POINTER :: qs_env
445 : REAL(KIND=dp), INTENT(OUT) :: gap, energy, condition_number
446 :
447 : INTEGER :: homo, nmo
448 : REAL(KIND=dp), DIMENSION(2) :: condnum
449 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
450 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
451 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
452 6 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: smat
453 6 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
454 : TYPE(qs_energy_type), POINTER :: qs_energy
455 :
456 6 : NULLIFY (blacs_env, eigenvalues, matrix_s, mos, qs_energy, smat)
457 6 : gap = 0.0_dp
458 6 : CALL get_qs_env(qs_env, blacs_env=blacs_env, energy=qs_energy, matrix_s=matrix_s, mos=mos)
459 6 : energy = qs_energy%total
460 18 : ALLOCATE (smat(1, 1))
461 6 : smat(1, 1)%matrix => matrix_s(1)%matrix
462 6 : CALL overlap_condnum(smat, condnum, -1, .FALSE., .TRUE., .FALSE., blacs_env)
463 6 : condition_number = condnum(2)
464 6 : DEALLOCATE (smat)
465 6 : IF (SIZE(mos) /= 1) THEN
466 0 : CPABORT("FRONTIER_ORBITAL_SCREENING requires a restricted calculation")
467 : END IF
468 6 : CALL get_mo_set(mos(1), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
469 6 : IF (homo < 1 .OR. homo >= nmo) THEN
470 0 : CPABORT("FRONTIER_ORBITAL_SCREENING could not determine a HOMO-LUMO gap")
471 : END IF
472 6 : gap = eigenvalues(homo + 1) - eigenvalues(homo)
473 :
474 12 : END SUBROUTINE get_frontier_orbital_scf_result
475 :
476 : ! **************************************************************************************************
477 : !> \brief Make a parsed reference input suitable for a quiet, independent screening SCF.
478 : !> \param input_file ...
479 : !> \param opt_bas ...
480 : ! **************************************************************************************************
481 12 : SUBROUTINE prepare_frontier_orbital_scf_input(input_file, opt_bas)
482 : TYPE(section_vals_type), POINTER :: input_file
483 : TYPE(basis_optimization_type), INTENT(IN) :: opt_bas
484 :
485 : CHARACTER(LEN=default_string_length) :: element
486 : CHARACTER(LEN=default_string_length), &
487 4 : DIMENSION(:), POINTER :: basis_info
488 : INTEGER :: ifit, ikind, nbasis_files, nkinds
489 : LOGICAL :: element_explicit, found
490 : TYPE(section_vals_type), POINTER :: dft_section, kind_section, subsys_section
491 :
492 4 : CALL section_vals_val_set(input_file, "GLOBAL%PRINT_LEVEL", i_val=silent_print_level)
493 4 : CALL section_vals_val_set(input_file, "FORCE_EVAL%DFT%SCF%SCF_GUESS", i_val=atomic_guess)
494 : CALL section_vals_val_set(input_file, &
495 : "FORCE_EVAL%DFT%SCF%PRINT%RESTART%_SECTION_PARAMETERS_", &
496 4 : i_val=debug_print_level + 1)
497 : CALL section_vals_val_set(input_file, &
498 : "FORCE_EVAL%DFT%SCF%PRINT%RESTART_HISTORY%_SECTION_PARAMETERS_", &
499 4 : i_val=debug_print_level + 1)
500 :
501 4 : dft_section => section_vals_get_subs_vals(input_file, "FORCE_EVAL%DFT")
502 4 : CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", n_rep_val=nbasis_files)
503 4 : IF (nbasis_files < 1) THEN
504 0 : CPABORT("The screening reference input has no basis-set file")
505 : END IF
506 : CALL section_vals_val_set( &
507 : dft_section, "BASIS_SET_FILE_NAME", i_rep_val=1, &
508 4 : c_val=discover_file(opt_bas%template_basis_file))
509 :
510 4 : subsys_section => section_vals_get_subs_vals(input_file, "FORCE_EVAL%SUBSYS")
511 4 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
512 4 : CALL section_vals_get(kind_section, n_repetition=nkinds)
513 12 : DO ikind = 1, nkinds
514 : CALL section_vals_val_get(kind_section, "ELEMENT", explicit=element_explicit, &
515 8 : i_rep_section=ikind)
516 8 : IF (element_explicit) THEN
517 0 : CALL section_vals_val_get(kind_section, "ELEMENT", c_val=element, i_rep_section=ikind)
518 : ELSE
519 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
520 8 : c_val=element, i_rep_section=ikind)
521 : END IF
522 8 : CALL uppercase(element)
523 8 : found = .FALSE.
524 12 : DO ifit = 1, opt_bas%nkind
525 12 : IF (TRIM(element) == TRIM(opt_bas%kind_basis(ifit)%element)) THEN
526 8 : NULLIFY (basis_info)
527 : CALL section_vals_val_get(kind_section, "BASIS_SET", c_vals=basis_info, &
528 8 : i_rep_val=1, i_rep_section=ikind)
529 8 : basis_info(SIZE(basis_info)) = opt_bas%kind_basis(ifit)%flex_basis(0)%basis_name
530 : found = .TRUE.
531 : EXIT
532 : END IF
533 : END DO
534 12 : IF (.NOT. found) THEN
535 : CALL cp_abort(__LOCATION__, &
536 : "Every KIND in the screening reference input must have a corresponding "// &
537 0 : "FIT_KIND: "//TRIM(element))
538 : END IF
539 : END DO
540 :
541 4 : END SUBROUTINE prepare_frontier_orbital_scf_input
542 :
543 : ! **************************************************************************************************
544 : !> \brief Write all optimized candidates to one CP2K basis file.
545 : !> \param filename ...
546 : !> \param candidates ...
547 : !> \param reference_gap ...
548 : !> \param reference_energy ...
549 : !> \param reference_condition_number ...
550 : !> \param para_env ...
551 : ! **************************************************************************************************
552 2 : SUBROUTINE write_frontier_orbital_basis_file(filename, candidates, &
553 : reference_gap, reference_energy, reference_condition_number, para_env)
554 : CHARACTER(LEN=*), INTENT(IN) :: filename
555 : TYPE(frontier_orbital_screening_candidate_type), &
556 : DIMENSION(:) :: candidates
557 : REAL(KIND=dp), INTENT(IN) :: reference_gap, reference_energy, &
558 : reference_condition_number
559 : TYPE(mp_para_env_type), POINTER :: para_env
560 :
561 : INTEGER :: candidate, ibasis, ikind, unit_nr
562 :
563 2 : IF (para_env%is_source()) THEN
564 1 : CALL open_file(file_name=filename, file_status="UNKNOWN", file_action="WRITE", unit_number=unit_nr)
565 : ELSE
566 1 : unit_nr = -999
567 : END IF
568 4 : DO candidate = 1, SIZE(candidates)
569 2 : IF (unit_nr > 0) THEN
570 1 : WRITE (unit_nr, '(A)') "#"
571 1 : WRITE (unit_nr, '(A,I0)') "# Optimized basis set ", candidate
572 1 : WRITE (unit_nr, '(A,F12.6,A)') "# Reference HOMO-LUMO gap: ", reference_gap*evolt, " eV"
573 1 : WRITE (unit_nr, '(A,F12.6,A)') "# Basis HOMO-LUMO gap: ", candidates(candidate)%gap*evolt, " eV"
574 1 : WRITE (unit_nr, '(A,F12.6,A)') "# Absolute gap error: ", &
575 2 : ABS(candidates(candidate)%gap - reference_gap)*evolt, " eV"
576 1 : WRITE (unit_nr, '(A,ES15.6)') "# Reference condition number: ", reference_condition_number
577 : WRITE (unit_nr, '(A,ES15.6)') &
578 1 : "# Basis condition number: ", candidates(candidate)%condition_number
579 1 : WRITE (unit_nr, '(A,F22.12,A)') "# Reference total energy: ", reference_energy, " hartree"
580 1 : WRITE (unit_nr, '(A,F22.12,A)') "# Basis total energy: ", candidates(candidate)%energy, " hartree"
581 1 : WRITE (unit_nr, '(A)') "#"
582 : END IF
583 8 : DO ikind = 1, candidates(candidate)%opt_bas%nkind
584 10 : DO ibasis = 0, candidates(candidate)%opt_bas%kind_basis(ikind)%nbasis_deriv
585 : candidates(candidate)%opt_bas%kind_basis(ikind)%flex_basis(ibasis)%basis_name = &
586 : TRIM(candidates(candidate)%opt_bas%kind_basis(ikind)%flex_basis(ibasis)%basis_name)// &
587 4 : "_"//TRIM(ADJUSTL(cp_to_string(candidate)))
588 : CALL write_basis(candidates(candidate)%opt_bas%kind_basis(ikind)%flex_basis(ibasis), &
589 8 : candidates(candidate)%opt_bas%kind_basis(ikind)%element, unit_nr)
590 : END DO
591 : END DO
592 : END DO
593 2 : IF (para_env%is_source()) CALL close_file(unit_number=unit_nr)
594 :
595 2 : END SUBROUTINE write_frontier_orbital_basis_file
596 :
597 : ! **************************************************************************************************
598 : !> \brief Print the two-line column headings shared by screening result tables.
599 : !> \param unit_nr ...
600 : ! **************************************************************************************************
601 1 : SUBROUTINE print_frontier_orbital_table_header(unit_nr)
602 : INTEGER, INTENT(IN) :: unit_nr
603 :
604 : WRITE (unit_nr, '(1X,A,T14,A,T24,A,T46,A,T59,A)') &
605 1 : "Basis set", "Gap / eV", "Abs. error to", "Condition", "Total energy / Ha"
606 1 : WRITE (unit_nr, '(T24,A,T46,A)') "ref. basis / meV", "number"
607 :
608 1 : END SUBROUTINE print_frontier_orbital_table_header
609 :
610 : ! **************************************************************************************************
611 : !> \brief Collect reference and initial basis-set names for every fitted atom kind.
612 : !> \param section ...
613 : !> \param opt_bas ...
614 : !> \param element_symbols ...
615 : !> \param reference_basis_names ...
616 : !> \param initial_basis_names ...
617 : ! **************************************************************************************************
618 1 : SUBROUTINE get_frontier_orbital_basis_sets(section, opt_bas, element_symbols, &
619 : reference_basis_names, initial_basis_names)
620 : TYPE(section_vals_type), POINTER :: section
621 : TYPE(basis_optimization_type), INTENT(IN) :: opt_bas
622 : CHARACTER(LEN=default_string_length), &
623 : ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: element_symbols, reference_basis_names, &
624 : initial_basis_names
625 :
626 : CHARACTER(LEN=default_string_length) :: element
627 : CHARACTER(LEN=default_string_length), &
628 1 : DIMENSION(:), POINTER :: basis_info
629 : INTEGER :: ifit, ikind, nkinds
630 : LOGICAL :: element_explicit, found
631 : TYPE(section_vals_type), POINTER :: kind_section, subsys_section
632 :
633 2 : subsys_section => section_vals_get_subs_vals(section, "FORCE_EVAL%SUBSYS")
634 1 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
635 1 : CALL section_vals_get(kind_section, n_repetition=nkinds)
636 :
637 0 : ALLOCATE (element_symbols(opt_bas%nkind), reference_basis_names(opt_bas%nkind), &
638 5 : initial_basis_names(opt_bas%nkind))
639 3 : DO ifit = 1, opt_bas%nkind
640 2 : found = .FALSE.
641 3 : DO ikind = 1, nkinds
642 : CALL section_vals_val_get(kind_section, "ELEMENT", explicit=element_explicit, &
643 3 : i_rep_section=ikind)
644 3 : IF (element_explicit) THEN
645 0 : CALL section_vals_val_get(kind_section, "ELEMENT", c_val=element, i_rep_section=ikind)
646 : ELSE
647 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
648 3 : c_val=element, i_rep_section=ikind)
649 : END IF
650 3 : CALL uppercase(element)
651 6 : IF (TRIM(element) == TRIM(opt_bas%kind_basis(ifit)%element)) THEN
652 2 : NULLIFY (basis_info)
653 : CALL section_vals_val_get(kind_section, "BASIS_SET", c_vals=basis_info, &
654 2 : i_rep_val=1, i_rep_section=ikind)
655 2 : element_symbols(ifit) = opt_bas%kind_basis(ifit)%element
656 2 : reference_basis_names(ifit) = basis_info(SIZE(basis_info))
657 2 : initial_basis_names(ifit) = opt_bas%kind_basis(ifit)%basis_name
658 : found = .TRUE.
659 : EXIT
660 : END IF
661 : END DO
662 1 : IF (.NOT. found) THEN
663 : CALL cp_abort(__LOCATION__, &
664 : "Could not find the reference basis set for fitted atom kind "// &
665 0 : TRIM(opt_bas%kind_basis(ifit)%element))
666 : END IF
667 : END DO
668 :
669 2 : END SUBROUTINE get_frontier_orbital_basis_sets
670 :
671 : ! **************************************************************************************************
672 : !> \brief Print the compact, aligned human-readable screening summary.
673 : !> \param unit_nr ...
674 : !> \param filename ...
675 : !> \param section ...
676 : !> \param candidates ...
677 : !> \param reference_gap ...
678 : !> \param reference_energy ...
679 : !> \param reference_condition_number ...
680 : !> \param initial_gap ...
681 : !> \param initial_energy ...
682 : !> \param initial_condition_number ...
683 : ! **************************************************************************************************
684 2 : SUBROUTINE print_frontier_orbital_summary(unit_nr, filename, section, candidates, &
685 : reference_gap, reference_energy, reference_condition_number, &
686 : initial_gap, initial_energy, initial_condition_number)
687 : INTEGER, INTENT(IN) :: unit_nr
688 : CHARACTER(LEN=*), INTENT(IN) :: filename
689 : TYPE(section_vals_type), POINTER :: section
690 : TYPE(frontier_orbital_screening_candidate_type), &
691 : DIMENSION(:), INTENT(IN) :: candidates
692 : REAL(KIND=dp), INTENT(IN) :: reference_gap, reference_energy, reference_condition_number, &
693 : initial_gap, initial_energy, initial_condition_number
694 :
695 : CHARACTER(LEN=default_string_length), &
696 2 : ALLOCATABLE, DIMENSION(:) :: element_symbols, initial_basis_names, &
697 2 : reference_basis_names
698 : INTEGER :: candidate, separator
699 :
700 2 : IF (unit_nr <= 0) RETURN
701 1 : separator = SCAN(TRIM(filename), "/", BACK=.TRUE.)
702 1 : WRITE (unit_nr, '(A)') ""
703 1 : WRITE (unit_nr, '(1X,A)') REPEAT("-", 79)
704 : WRITE (unit_nr, '(1X,A)') &
705 1 : "SCF calcs. with ref. basis set, initial small and optimized basis sets"
706 1 : WRITE (unit_nr, '(1X,A)') REPEAT("-", 79)
707 1 : WRITE (unit_nr, '(A)') ""
708 : CALL get_frontier_orbital_basis_sets(section, candidates(1)%opt_bas, element_symbols, &
709 1 : reference_basis_names, initial_basis_names)
710 : CALL print_frontier_orbital_basis_sets(unit_nr, element_symbols, reference_basis_names, &
711 1 : initial_basis_names)
712 1 : DEALLOCATE (element_symbols, reference_basis_names, initial_basis_names)
713 1 : WRITE (unit_nr, '(A)') ""
714 1 : CALL print_frontier_orbital_table_header(unit_nr)
715 1 : WRITE (unit_nr, '(A)') ""
716 : WRITE (unit_nr, '(1X,A,T14,F8.3,T24,F16.1,T46,ES11.3,T59,F17.12)') &
717 1 : "Reference", reference_gap*evolt, 0.0_dp, reference_condition_number, reference_energy
718 : WRITE (unit_nr, '(1X,A,T14,F8.3,T24,F16.1,T46,ES11.3,T59,F17.12)') &
719 1 : "Initial", initial_gap*evolt, ABS(initial_gap - reference_gap)*evolt*1000.0_dp, &
720 2 : initial_condition_number, initial_energy
721 2 : DO candidate = 1, SIZE(candidates)
722 : WRITE (unit_nr, '(1X,A,I0,T14,F8.3,T24,F16.1,T46,ES11.3,T59,F17.12)') &
723 1 : "Basis ", candidate, candidates(candidate)%gap*evolt, &
724 1 : ABS(candidates(candidate)%gap - reference_gap)*evolt*1000.0_dp, &
725 3 : candidates(candidate)%condition_number, candidates(candidate)%energy
726 : END DO
727 1 : WRITE (unit_nr, '(A)') ""
728 1 : WRITE (unit_nr, '(1X,A)') REPEAT("-", 79)
729 1 : WRITE (unit_nr, '(A)') ""
730 : WRITE (unit_nr, '(1X,A,A)') &
731 1 : "Basis set file containing all optimized basis sets: ", TRIM(filename(separator + 1:))
732 1 : WRITE (unit_nr, '(A)') ""
733 1 : WRITE (unit_nr, '(1X,A)') REPEAT("=", 79)
734 1 : WRITE (unit_nr, '(A)') ""
735 :
736 2 : END SUBROUTINE print_frontier_orbital_summary
737 :
738 : ! **************************************************************************************************
739 : !> \brief Join a screening working directory and filename.
740 : !> \param directory ...
741 : !> \param filename ...
742 : !> \return ...
743 : ! **************************************************************************************************
744 2 : PURE FUNCTION work_path(directory, filename) RESULT(path)
745 : CHARACTER(LEN=*), INTENT(IN) :: directory, filename
746 : CHARACTER(LEN=default_path_length) :: path
747 :
748 : INTEGER :: n
749 :
750 2 : n = LEN_TRIM(directory)
751 2 : IF (filename(1:1) == "/") THEN
752 0 : path = filename
753 2 : ELSE IF (n == 0 .OR. directory(n:n) == "/") THEN
754 0 : path = TRIM(directory)//filename
755 : ELSE
756 2 : path = TRIM(directory)//"/"//filename
757 : END IF
758 :
759 2 : END FUNCTION work_path
760 :
761 0 : END MODULE optbas_frontier_orbital_screening
|
