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 : MODULE optimize_basis
8 : USE cp_blacs_env, ONLY: cp_blacs_env_type
9 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
10 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
11 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
12 : cp_fm_struct_release,&
13 : cp_fm_struct_type
14 : USE cp_fm_types, ONLY: cp_fm_release,&
15 : cp_fm_type
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_get_default_unit_nr,&
18 : cp_logger_type
19 : USE f77_interface, ONLY: create_force_env,&
20 : destroy_force_env,&
21 : f_env_add_defaults,&
22 : f_env_get_from_id,&
23 : f_env_rm_defaults,&
24 : f_env_type
25 : USE force_env_types, ONLY: force_env_get,&
26 : force_env_type
27 : USE input_cp2k_read, ONLY: empty_initial_variables,&
28 : read_input
29 : USE input_section_types, ONLY: section_type,&
30 : section_vals_release,&
31 : section_vals_type
32 : USE kinds, ONLY: default_path_length,&
33 : dp
34 : USE machine, ONLY: m_chdir,&
35 : m_getcwd,&
36 : m_walltime
37 : USE message_passing, ONLY: mp_comm_type,&
38 : mp_para_env_release,&
39 : mp_para_env_type
40 : USE optbas_fenv_manipulation, ONLY: allocate_mo_sets,&
41 : calculate_ks_matrix,&
42 : calculate_overlap_inverse,&
43 : modify_input_settings,&
44 : update_basis_set
45 : USE optbas_opt_utils, ONLY: evaluate_optvals,&
46 : fit_mo_coeffs,&
47 : optbas_build_neighborlist
48 : USE optimize_basis_types, ONLY: basis_optimization_type,&
49 : deallocate_basis_optimization_type,&
50 : subset_type
51 : USE optimize_basis_utils, ONLY: get_set_and_basis_id,&
52 : optimize_basis_init_read_input,&
53 : update_derived_basis_sets
54 : USE powell, ONLY: powell_optimize
55 : USE qs_environment_types, ONLY: get_qs_env,&
56 : qs_env_part_release,&
57 : qs_environment_type
58 : USE qs_kind_types, ONLY: get_qs_kind_set,&
59 : qs_kind_type
60 : USE qs_ks_types, ONLY: get_ks_env,&
61 : qs_ks_env_type,&
62 : set_ks_env
63 : USE qs_mo_types, ONLY: allocate_mo_set,&
64 : deallocate_mo_set,&
65 : get_mo_set,&
66 : init_mo_set,&
67 : mo_set_type
68 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
69 : release_neighbor_list_sets
70 : USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists
71 : USE qs_overlap, ONLY: build_overlap_matrix
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 : PRIVATE
76 :
77 : PUBLIC :: run_optimize_basis
78 :
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_basis'
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief main entry point for methods aimed at optimizing basis sets
85 : !> \param input_declaration ...
86 : !> \param root_section ...
87 : !> \param para_env ...
88 : !> \author Florian Schiffmann
89 : ! **************************************************************************************************
90 4 : SUBROUTINE run_optimize_basis(input_declaration, root_section, para_env)
91 : TYPE(section_type), POINTER :: input_declaration
92 : TYPE(section_vals_type), POINTER :: root_section
93 : TYPE(mp_para_env_type), POINTER :: para_env
94 :
95 : CHARACTER(len=*), PARAMETER :: routineN = 'run_optimize_basis'
96 :
97 : INTEGER :: handle
98 4 : TYPE(basis_optimization_type) :: opt_bas
99 :
100 4 : CALL timeset(routineN, handle)
101 :
102 4 : CALL optimize_basis_init_read_input(opt_bas, root_section, para_env)
103 :
104 4 : CALL driver_para_opt_basis(opt_bas, input_declaration, para_env)
105 :
106 4 : CALL deallocate_basis_optimization_type(opt_bas)
107 4 : CALL timestop(handle)
108 :
109 4 : END SUBROUTINE run_optimize_basis
110 :
111 : ! **************************************************************************************************
112 : !> \brief driver routine for the parallel part of the method
113 : !> \param opt_bas ...
114 : !> \param input_declaration ...
115 : !> \param para_env ...
116 : !> \author Florian Schiffmann
117 : ! **************************************************************************************************
118 :
119 4 : SUBROUTINE driver_para_opt_basis(opt_bas, input_declaration, para_env)
120 : TYPE(basis_optimization_type) :: opt_bas
121 : TYPE(section_type), POINTER :: input_declaration
122 : TYPE(mp_para_env_type), POINTER :: para_env
123 :
124 : CHARACTER(len=*), PARAMETER :: routineN = 'driver_para_opt_basis'
125 :
126 : INTEGER :: handle, n_groups_created
127 : TYPE(mp_comm_type) :: opt_group
128 : INTEGER, DIMENSION(:), POINTER :: group_distribution_p
129 8 : INTEGER, DIMENSION(0:para_env%num_pe-1), TARGET :: group_distribution
130 :
131 4 : CALL timeset(routineN, handle)
132 4 : group_distribution_p => group_distribution
133 : CALL opt_group%from_split(para_env, n_groups_created, group_distribution_p, &
134 4 : n_subgroups=SIZE(opt_bas%group_partition), group_partition=opt_bas%group_partition)
135 4 : opt_bas%opt_id = group_distribution(para_env%mepos) + 1
136 4 : opt_bas%n_groups_created = n_groups_created
137 12 : ALLOCATE (opt_bas%sub_sources(0:para_env%num_pe - 1))
138 :
139 4 : CALL driver_optimization_para_low(opt_bas, input_declaration, para_env, opt_group)
140 :
141 4 : CALL opt_group%free()
142 4 : CALL timestop(handle)
143 :
144 4 : END SUBROUTINE driver_para_opt_basis
145 :
146 : ! **************************************************************************************************
147 : !> \brief low level optimization routine includes initialization of the subsytems
148 : !> powell optimizer and deallocation of the various force envs
149 : !> \param opt_bas ...
150 : !> \param input_declaration ...
151 : !> \param para_env_top ...
152 : !> \param mpi_comm_opt ...
153 : !> \author Florian Schiffmann
154 : ! **************************************************************************************************
155 :
156 4 : SUBROUTINE driver_optimization_para_low(opt_bas, input_declaration, para_env_top, mpi_comm_opt)
157 : TYPE(basis_optimization_type) :: opt_bas
158 : TYPE(section_type), POINTER :: input_declaration
159 : TYPE(mp_para_env_type), POINTER :: para_env_top
160 : TYPE(mp_comm_type), INTENT(IN) :: mpi_comm_opt
161 :
162 : CHARACTER(len=*), PARAMETER :: routineN = 'driver_optimization_para_low'
163 :
164 : INTEGER :: handle, icalc, iopt, is, mp_id, stat
165 : INTEGER, ALLOCATABLE, DIMENSION(:) :: f_env_id
166 : LOGICAL :: write_basis
167 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tot_time
168 4 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: matrix_S_inv
169 : TYPE(f_env_type), POINTER :: f_env
170 : TYPE(mp_para_env_type), POINTER :: para_env
171 :
172 4 : NULLIFY (f_env)
173 :
174 4 : CALL timeset(routineN, handle)
175 :
176 : ! ====== initialize the f_env and precompute some matrices =====
177 4 : mp_id = opt_bas%opt_id
178 4 : NULLIFY (para_env, f_env)
179 12 : ALLOCATE (f_env_id(SIZE(opt_bas%comp_group(mp_id)%member_list)))
180 12 : ALLOCATE (tot_time(opt_bas%ncombinations*opt_bas%ntraining_sets))
181 21 : ALLOCATE (matrix_s_inv(SIZE(opt_bas%comp_group(mp_id)%member_list)))
182 :
183 4 : ALLOCATE (para_env)
184 4 : para_env = mpi_comm_opt
185 :
186 4 : is = -1
187 4 : IF (para_env%is_source()) is = para_env_top%mepos
188 4 : CALL para_env_top%allgather(is, opt_bas%sub_sources)
189 :
190 4 : CALL init_training_force_envs(opt_bas, f_env_id, input_declaration, matrix_s_inv, para_env, mpi_comm_opt)
191 :
192 4 : CALL init_free_vars(opt_bas)
193 4 : tot_time = 0.0_dp
194 :
195 : ! ======= The real optimization loop =======
196 118 : DO iopt = 0, opt_bas%powell_param%maxfun
197 : CALL compute_residuum_vectors(opt_bas, f_env_id, matrix_S_inv, tot_time, &
198 114 : para_env_top, para_env, iopt)
199 114 : IF (para_env_top%is_source()) THEN
200 57 : CALL powell_optimize(opt_bas%powell_param%nvar, opt_bas%x_opt, opt_bas%powell_param)
201 : END IF
202 114 : CALL para_env_top%bcast(opt_bas%powell_param%state)
203 114 : CALL para_env_top%bcast(opt_bas%x_opt)
204 114 : CALL update_free_vars(opt_bas)
205 114 : write_basis = MOD(iopt, opt_bas%write_frequency) == 0
206 : CALL update_derived_basis_sets(opt_bas, write_basis, opt_bas%output_basis_file, &
207 114 : para_env_top)
208 118 : IF (opt_bas%powell_param%state == -1) EXIT
209 : END DO
210 :
211 : ! ======= Update the basis set and print the final basis =======
212 4 : IF (para_env_top%is_source()) THEN
213 2 : opt_bas%powell_param%state = 8
214 2 : CALL powell_optimize(opt_bas%powell_param%nvar, opt_bas%x_opt, opt_bas%powell_param)
215 : END IF
216 :
217 4 : CALL para_env_top%bcast(opt_bas%x_opt)
218 4 : CALL update_free_vars(opt_bas)
219 : CALL update_derived_basis_sets(opt_bas, .TRUE., opt_bas%output_basis_file, &
220 4 : para_env_top)
221 :
222 : ! ====== get rid of the f_env again =====
223 :
224 13 : DO icalc = SIZE(opt_bas%comp_group(mp_id)%member_list), 1, -1
225 9 : CALL f_env_get_from_id(f_env_id(icalc), f_env)
226 13 : CALL destroy_force_env(f_env_id(icalc), stat)
227 : END DO
228 4 : DEALLOCATE (f_env_id); DEALLOCATE (tot_time)
229 4 : CALL cp_fm_release(matrix_s_inv)
230 4 : CALL mp_para_env_release(para_env)
231 4 : CALL timestop(handle)
232 :
233 8 : END SUBROUTINE driver_optimization_para_low
234 :
235 : ! **************************************************************************************************
236 : !> \brief compute all ingredients for powell optimizer. Rho_diff,
237 : !> condition number, energy,... for all ttraining sets in
238 : !> the computational group
239 : !> \param opt_bas ...
240 : !> \param f_env_id ...
241 : !> \param matrix_S_inv ...
242 : !> \param tot_time ...
243 : !> \param para_env_top ...
244 : !> \param para_env ...
245 : !> \param iopt ...
246 : ! **************************************************************************************************
247 :
248 114 : SUBROUTINE compute_residuum_vectors(opt_bas, f_env_id, matrix_S_inv, tot_time, &
249 : para_env_top, para_env, iopt)
250 : TYPE(basis_optimization_type) :: opt_bas
251 : INTEGER, ALLOCATABLE, DIMENSION(:) :: f_env_id
252 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: matrix_S_inv
253 : REAL(KIND=dp), DIMENSION(:) :: tot_time
254 : TYPE(mp_para_env_type), POINTER :: para_env_top, para_env
255 : INTEGER :: iopt
256 :
257 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_residuum_vectors'
258 :
259 : CHARACTER(len=8) :: basis_type
260 : INTEGER :: bas_id, handle, icalc, icomb, ispin, &
261 : mp_id, my_id, nao, ncalc, nelectron, &
262 : nmo, nspins, set_id
263 : REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f
264 114 : REAL(KIND=dp), DIMENSION(:), POINTER :: cond_vec, energy, f_vec, my_time, &
265 114 : start_time
266 114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gdata
267 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
268 : TYPE(cp_fm_type), POINTER :: mo_coeff
269 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s_aux, matrix_s_aux_orb
270 : TYPE(f_env_type), POINTER :: f_env
271 : TYPE(force_env_type), POINTER :: force_env
272 114 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_aux
273 114 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
274 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
275 114 : POINTER :: sab_aux, sab_aux_orb
276 : TYPE(qs_environment_type), POINTER :: qs_env
277 114 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
278 : TYPE(qs_ks_env_type), POINTER :: ks_env
279 :
280 114 : CALL timeset(routineN, handle)
281 :
282 114 : basis_type = "AUX_OPT"
283 : !
284 114 : ncalc = opt_bas%ncombinations*opt_bas%ntraining_sets
285 342 : ALLOCATE (gdata(ncalc, 4))
286 114 : f_vec => gdata(:, 1)
287 114 : my_time => gdata(:, 2)
288 114 : cond_vec => gdata(:, 3)
289 114 : energy => gdata(:, 4)
290 : !
291 1986 : f_vec = 0.0_dp; cond_vec = 0.0_dp; my_time = 0.0_dp; energy = 0.0_dp
292 114 : mp_id = opt_bas%opt_id
293 342 : ALLOCATE (start_time(SIZE(opt_bas%comp_group(mp_id)%member_list)))
294 : !
295 348 : DO icalc = 1, SIZE(opt_bas%comp_group(mp_id)%member_list)
296 234 : my_id = opt_bas%comp_group(mp_id)%member_list(icalc) + 1
297 : ! setup timings
298 234 : start_time(icalc) = m_walltime()
299 :
300 234 : NULLIFY (matrix_s_aux_orb, matrix_s_aux)
301 234 : CALL get_set_and_basis_id(opt_bas%comp_group(mp_id)%member_list(icalc), opt_bas, set_id, bas_id)
302 234 : CALL f_env_get_from_id(f_env_id(icalc), f_env)
303 234 : force_env => f_env%force_env
304 234 : CALL force_env_get(force_env, qs_env=qs_env)
305 234 : CALL get_qs_env(qs_env, ks_env=ks_env)
306 234 : CALL update_basis_set(opt_bas, bas_id, basis_type, qs_env)
307 234 : NULLIFY (sab_aux, sab_aux_orb)
308 234 : CALL optbas_build_neighborlist(qs_env, sab_aux, sab_aux_orb, basis_type)
309 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s_aux, &
310 : basis_type_a=basis_type, &
311 : basis_type_b=basis_type, &
312 234 : sab_nl=sab_aux)
313 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s_aux_orb, &
314 : basis_type_a=basis_type, &
315 : basis_type_b="ORB", &
316 234 : sab_nl=sab_aux_orb)
317 234 : CALL release_neighbor_list_sets(sab_aux)
318 234 : CALL release_neighbor_list_sets(sab_aux_orb)
319 234 : CALL get_qs_env(qs_env, mos=mos, matrix_ks=matrix_ks)
320 :
321 234 : nspins = SIZE(mos)
322 936 : ALLOCATE (mos_aux(nspins))
323 234 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
324 234 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao, basis_type=basis_type)
325 468 : DO ispin = 1, nspins
326 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, maxocc=maxocc, nelectron=nelectron, &
327 234 : n_el_f=n_el_f, nmo=nmo, flexible_electron_count=flexible_electron_count)
328 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo, &
329 : context=mo_coeff%matrix_struct%context, &
330 234 : para_env=mo_coeff%matrix_struct%para_env)
331 : CALL allocate_mo_set(mos_aux(ispin), nao, nmo, nelectron, &
332 234 : n_el_f, maxocc, flexible_electron_count)
333 234 : CALL init_mo_set(mo_set=mos_aux(ispin), fm_struct=fm_struct, name="MO_AUX")
334 702 : CALL cp_fm_struct_release(fm_struct)
335 : END DO
336 :
337 234 : CALL fit_mo_coeffs(matrix_s_aux, matrix_s_aux_orb, mos, mos_aux)
338 : CALL evaluate_optvals(mos, mos_aux, matrix_ks, matrix_s_aux_orb(1)%matrix, &
339 : matrix_s_aux(1)%matrix, matrix_S_inv(icalc), &
340 234 : f_vec(my_id), energy(my_id), cond_vec(my_id))
341 :
342 468 : DO ispin = 1, nspins
343 468 : CALL deallocate_mo_set(mos_aux(ispin))
344 : END DO
345 234 : DEALLOCATE (mos_aux)
346 234 : IF (ASSOCIATED(matrix_s_aux)) CALL dbcsr_deallocate_matrix_set(matrix_s_aux)
347 234 : IF (ASSOCIATED(matrix_s_aux_orb)) CALL dbcsr_deallocate_matrix_set(matrix_s_aux_orb)
348 :
349 582 : my_time(my_id) = m_walltime() - start_time(icalc)
350 : END DO
351 :
352 114 : DEALLOCATE (start_time)
353 :
354 114 : IF (.NOT. para_env%is_source()) THEN
355 0 : f_vec = 0.0_dp; cond_vec = 0.0_dp; my_time = 0.0_dp; energy = 0.0_dp
356 : END IF
357 : ! collect date from all subgroup ionodes on the main ionode
358 4770 : CALL para_env_top%sum(gdata)
359 :
360 114 : opt_bas%powell_param%f = 0.0_dp
361 114 : IF (para_env_top%is_source()) THEN
362 291 : DO icalc = 1, SIZE(f_vec)
363 234 : icomb = MOD(icalc - 1, opt_bas%ncombinations)
364 : opt_bas%powell_param%f = opt_bas%powell_param%f + &
365 234 : (f_vec(icalc) + energy(icalc))*opt_bas%fval_weight(icomb)
366 291 : IF (opt_bas%use_condition_number) THEN
367 : opt_bas%powell_param%f = opt_bas%powell_param%f + &
368 234 : LOG(cond_vec(icalc))*opt_bas%condition_weight(icomb)
369 : END IF
370 : END DO
371 : ELSE
372 993 : f_vec = 0.0_dp; cond_vec = 0.0_dp; my_time = 0.0_dp; energy = 0.0_dp
373 : END IF
374 114 : CALL para_env_top%bcast(opt_bas%powell_param%f)
375 :
376 : ! output info if required
377 114 : CALL output_opt_info(f_vec, cond_vec, my_time, tot_time, opt_bas, iopt, para_env_top)
378 114 : DEALLOCATE (gdata)
379 :
380 114 : CALL para_env_top%sync()
381 :
382 114 : CALL timestop(handle)
383 :
384 228 : END SUBROUTINE compute_residuum_vectors
385 :
386 : ! **************************************************************************************************
387 : !> \brief create the force_envs for every input in the computational group
388 : !> \param opt_bas ...
389 : !> \param f_env_id ...
390 : !> \param input_declaration ...
391 : !> \param matrix_s_inv ...
392 : !> \param para_env ...
393 : !> \param mpi_comm_opt ...
394 : ! **************************************************************************************************
395 :
396 17 : SUBROUTINE init_training_force_envs(opt_bas, f_env_id, input_declaration, matrix_s_inv, para_env, mpi_comm_opt)
397 :
398 : TYPE(basis_optimization_type) :: opt_bas
399 : INTEGER, ALLOCATABLE, DIMENSION(:) :: f_env_id
400 : TYPE(section_type), POINTER :: input_declaration
401 : TYPE(cp_fm_type), DIMENSION(:), INTENT(OUT) :: matrix_S_inv
402 : TYPE(mp_para_env_type), POINTER :: para_env
403 : TYPE(mp_comm_type) :: mpi_comm_opt
404 :
405 : CHARACTER(len=*), PARAMETER :: routineN = 'init_training_force_envs'
406 :
407 : CHARACTER(len=default_path_length) :: main_dir
408 : INTEGER :: bas_id, handle, icalc, ierr, mp_id, &
409 : set_id, stat
410 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
411 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
412 : TYPE(f_env_type), POINTER :: f_env
413 : TYPE(force_env_type), POINTER :: force_env
414 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
415 4 : POINTER :: sab_orb
416 : TYPE(qs_environment_type), POINTER :: qs_env
417 : TYPE(qs_ks_env_type), POINTER :: ks_env
418 : TYPE(section_vals_type), POINTER :: input_file
419 :
420 4 : CALL timeset(routineN, handle)
421 :
422 4 : NULLIFY (matrix_s, blacs_env, ks_env)
423 :
424 4 : mp_id = opt_bas%opt_id
425 4 : CALL m_getcwd(main_dir)
426 :
427 : ! ======= Create f_env for all calculations in MPI group =======
428 13 : DO icalc = 1, SIZE(opt_bas%comp_group(mp_id)%member_list)
429 9 : NULLIFY (input_file)
430 : ! parse the input of the training sets
431 9 : CALL get_set_and_basis_id(opt_bas%comp_group(mp_id)%member_list(icalc), opt_bas, set_id, bas_id)
432 9 : CALL m_chdir(TRIM(opt_bas%training_dir(set_id)), ierr)
433 9 : IF (ierr /= 0) THEN
434 : CALL cp_abort(__LOCATION__, &
435 0 : "Could not change to directory <"//TRIM(opt_bas%training_dir(set_id))//">")
436 : END IF
437 : input_file => read_input(input_declaration, &
438 : opt_bas%training_input(set_id), &
439 : initial_variables=empty_initial_variables, &
440 9 : para_env=para_env)
441 :
442 9 : CALL modify_input_settings(opt_bas, bas_id, input_file)
443 : CALL create_force_env(f_env_id(icalc), &
444 : input_declaration=input_declaration, &
445 : input_path=opt_bas%training_input(set_id), &
446 : input=input_file, &
447 : output_path="scrap_information", &
448 : mpi_comm=mpi_comm_opt, &
449 9 : ierr=stat)
450 :
451 : ! some weirdness with the default stacks defaults have to be addded to get the
452 : ! correct default program name this causes trouble with the timer stack if kept
453 9 : CALL f_env_add_defaults(f_env_id(icalc), f_env)
454 9 : force_env => f_env%force_env
455 9 : CALL force_env_get(force_env, qs_env=qs_env)
456 9 : CALL allocate_mo_sets(qs_env)
457 9 : CALL f_env_rm_defaults(f_env, stat)
458 9 : CALL get_qs_env(qs_env, ks_env=ks_env)
459 : CALL build_qs_neighbor_lists(qs_env, para_env, molecular=.FALSE., &
460 9 : force_env_section=qs_env%input)
461 : CALL get_ks_env(ks_env, &
462 : matrix_s=matrix_s, &
463 9 : sab_orb=sab_orb)
464 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s, &
465 : matrix_name="OVERLAP", &
466 : basis_type_a="ORB", &
467 : basis_type_b="ORB", &
468 9 : sab_nl=sab_orb)
469 9 : CALL set_ks_env(ks_env, matrix_s=matrix_s)
470 9 : CALL get_qs_env(qs_env, matrix_s=matrix_s, blacs_env=blacs_env)
471 : CALL calculate_overlap_inverse(matrix_s(1)%matrix, matrix_s_inv(icalc), &
472 9 : para_env, blacs_env)
473 9 : CALL calculate_ks_matrix(qs_env)
474 :
475 9 : CALL section_vals_release(input_file)
476 :
477 9 : CALL qs_env_part_release(qs_env)
478 :
479 22 : CALL m_chdir(TRIM(ADJUSTL(main_dir)), ierr)
480 : END DO
481 :
482 4 : CALL timestop(handle)
483 :
484 4 : END SUBROUTINE init_training_force_envs
485 :
486 : ! **************************************************************************************************
487 : !> \brief variable update from the powell vector for all sets
488 : !> \param opt_bas ...
489 : !> \author Florian Schiffmann
490 : ! **************************************************************************************************
491 :
492 118 : SUBROUTINE update_free_vars(opt_bas)
493 : TYPE(basis_optimization_type) :: opt_bas
494 :
495 : CHARACTER(len=*), PARAMETER :: routineN = 'update_free_vars'
496 :
497 : INTEGER :: handle, ikind, iset, ix
498 :
499 118 : CALL timeset(routineN, handle)
500 118 : ix = 0
501 354 : DO ikind = 1, opt_bas%nkind
502 590 : DO iset = 1, opt_bas%kind_basis(ikind)%flex_basis(0)%nsets
503 472 : CALL update_subset_freevars(opt_bas%kind_basis(ikind)%flex_basis(0)%subset(iset), ix, opt_bas%x_opt)
504 : END DO
505 : END DO
506 118 : CALL timestop(handle)
507 :
508 118 : END SUBROUTINE update_free_vars
509 :
510 : ! **************************************************************************************************
511 : !> \brief low level update for the basis sets. Exponents are transformed according to constraint
512 : !> \param subset ...
513 : !> \param ix ...
514 : !> \param x ...
515 : !> \author Florian Schiffmann
516 : ! **************************************************************************************************
517 :
518 236 : SUBROUTINE update_subset_freevars(subset, ix, x)
519 : TYPE(subset_type) :: subset
520 : INTEGER :: ix
521 : REAL(KIND=dp), DIMENSION(:) :: x
522 :
523 : CHARACTER(len=*), PARAMETER :: routineN = 'update_subset_freevars'
524 :
525 : INTEGER :: handle, icon1, icon2, icont, iexp, il, &
526 : istart
527 : REAL(KIND=dp) :: fermi_f, gs_scale
528 :
529 236 : CALL timeset(routineN, handle)
530 1888 : DO iexp = 1, subset%nexp
531 1652 : IF (subset%opt_exps(iexp)) THEN
532 0 : ix = ix + 1
533 0 : subset%exps(iexp) = ABS(x(ix))
534 0 : IF (subset%exp_has_const(iexp)) THEN
535 : !use a fermi function to keep exponents in a given range around their initial value
536 0 : fermi_f = 1.0_dp/(EXP((x(ix) - 1.0_dp)/0.5_dp) + 1.0_dp)
537 : subset%exps(iexp) = (2.0_dp*fermi_f - 1.0_dp)*subset%exp_const(iexp)%var_fac*subset%exp_const(iexp)%init + &
538 0 : subset%exp_const(iexp)%init
539 : ELSE
540 :
541 : END IF
542 : END IF
543 10974 : DO icont = 1, subset%ncon_tot
544 10738 : IF (subset%opt_coeff(iexp, icont)) THEN
545 9086 : ix = ix + 1
546 9086 : subset%coeff(iexp, icont) = x(ix)
547 : END IF
548 : END DO
549 : END DO
550 :
551 : ! orthonormalize contraction coefficients using gram schmidt
552 236 : istart = 1
553 826 : DO il = 1, subset%nl
554 1298 : DO icon1 = istart, istart + subset%l(il) - 2
555 2360 : DO icon2 = icon1 + 1, istart + subset%l(il) - 1
556 : gs_scale = DOT_PRODUCT(subset%coeff(:, icon2), subset%coeff(:, icon1))/ &
557 15930 : DOT_PRODUCT(subset%coeff(:, icon1), subset%coeff(:, icon1))
558 9204 : subset%coeff(:, icon2) = subset%coeff(:, icon2) - gs_scale*subset%coeff(:, icon1)
559 : END DO
560 : END DO
561 826 : istart = istart + subset%l(il)
562 : END DO
563 :
564 1534 : DO icon1 = 1, subset%ncon_tot
565 19706 : subset%coeff(:, icon1) = subset%coeff(:, icon1)/NORM2(subset%coeff(:, icon1))
566 : END DO
567 236 : CALL timestop(handle)
568 :
569 236 : END SUBROUTINE update_subset_freevars
570 :
571 : ! **************************************************************************************************
572 : !> \brief variable initialization for the powell vector for all sets
573 : !> \param opt_bas ...
574 : !> \author Florian Schiffmann
575 : ! **************************************************************************************************
576 :
577 4 : SUBROUTINE init_free_vars(opt_bas)
578 : TYPE(basis_optimization_type) :: opt_bas
579 :
580 : CHARACTER(len=*), PARAMETER :: routineN = 'init_free_vars'
581 :
582 : INTEGER :: handle, ikind, iset, ix
583 :
584 4 : CALL timeset(routineN, handle)
585 4 : ix = 0
586 12 : DO ikind = 1, opt_bas%nkind
587 20 : DO iset = 1, opt_bas%kind_basis(ikind)%flex_basis(0)%nsets
588 16 : CALL init_subset_freevars(opt_bas%kind_basis(ikind)%flex_basis(0)%subset(iset), ix, opt_bas%x_opt)
589 : END DO
590 : END DO
591 4 : CALL timestop(handle)
592 :
593 4 : END SUBROUTINE init_free_vars
594 :
595 : ! **************************************************************************************************
596 : !> \brief variable initialization for the powell vector from low level informations
597 : !> constraint exponents will be mapped on a fermi function
598 : !> \param subset ...
599 : !> \param ix ...
600 : !> \param x ...
601 : !> \author Florian Schiffmann
602 : ! **************************************************************************************************
603 :
604 8 : SUBROUTINE init_subset_freevars(subset, ix, x)
605 : TYPE(subset_type) :: subset
606 : INTEGER :: ix
607 : REAL(KIND=dp), DIMENSION(:) :: x
608 :
609 : CHARACTER(len=*), PARAMETER :: routineN = 'init_subset_freevars'
610 :
611 : INTEGER :: handle, icont, iexp
612 : REAL(KIND=dp) :: fract
613 :
614 8 : CALL timeset(routineN, handle)
615 :
616 64 : DO iexp = 1, subset%nexp
617 56 : IF (subset%opt_exps(iexp)) THEN
618 0 : ix = ix + 1
619 0 : x(ix) = subset%exps(iexp)
620 0 : IF (subset%exp_has_const(iexp)) THEN
621 0 : IF (subset%exp_const(iexp)%const_type == 0) THEN
622 : fract = 1.0_dp + (subset%exps(iexp) - subset%exp_const(iexp)%init)/ &
623 0 : (subset%exp_const(iexp)%init*subset%exp_const(iexp)%var_fac)
624 0 : x(ix) = 0.5_dp*LOG((2.0_dp/fract - 1.0_dp)) + 1.0_dp
625 : END IF
626 0 : IF (subset%exp_const(iexp)%const_type == 1) THEN
627 0 : x(ix) = 1.0_dp
628 : END IF
629 : END IF
630 : END IF
631 372 : DO icont = 1, subset%ncon_tot
632 364 : IF (subset%opt_coeff(iexp, icont)) THEN
633 308 : ix = ix + 1
634 308 : x(ix) = subset%coeff(iexp, icont)
635 : END IF
636 : END DO
637 : END DO
638 8 : CALL timestop(handle)
639 :
640 8 : END SUBROUTINE init_subset_freevars
641 :
642 : ! **************************************************************************************************
643 : !> \brief commuticates all info to the master and assembles the output
644 : !> \param f_vec ...
645 : !> \param cond_vec ...
646 : !> \param my_time ...
647 : !> \param tot_time ...
648 : !> \param opt_bas ...
649 : !> \param iopt ...
650 : !> \param para_env_top ...
651 : !> \author Florian Schiffmann
652 : ! **************************************************************************************************
653 :
654 114 : SUBROUTINE output_opt_info(f_vec, cond_vec, my_time, tot_time, opt_bas, iopt, para_env_top)
655 : REAL(KIND=dp), DIMENSION(:) :: f_vec, cond_vec, my_time, tot_time
656 : TYPE(basis_optimization_type) :: opt_bas
657 : INTEGER :: iopt
658 : TYPE(mp_para_env_type), POINTER :: para_env_top
659 :
660 : CHARACTER(len=*), PARAMETER :: routineN = 'output_opt_info'
661 :
662 : INTEGER :: handle, ibasis, icalc, iset, unit_nr
663 : TYPE(cp_logger_type), POINTER :: logger
664 :
665 114 : CALL timeset(routineN, handle)
666 114 : logger => cp_get_default_logger()
667 :
668 582 : tot_time = tot_time + my_time
669 :
670 114 : unit_nr = -1
671 114 : IF (para_env_top%is_source() .AND. (MOD(iopt, opt_bas%write_frequency) == 0 .OR. iopt == opt_bas%powell_param%maxfun)) THEN
672 5 : unit_nr = cp_logger_get_default_unit_nr(logger)
673 : END IF
674 :
675 5 : IF (unit_nr > 0) THEN
676 5 : WRITE (unit_nr, '(1X,A,I8)') "BASOPT| Information at iteration number:", iopt
677 5 : WRITE (unit_nr, '(1X,A)') "BASOPT| Training set | Combination | Rho difference | Condition num. | Time"
678 5 : WRITE (unit_nr, '(1X,A)') "BASOPT| -----------------------------------------------------------------------"
679 5 : icalc = 0
680 12 : DO iset = 1, opt_bas%ntraining_sets
681 33 : DO ibasis = 1, opt_bas%ncombinations
682 21 : icalc = icalc + 1
683 : WRITE (unit_nr, '(1X,A,2(5X,I3,5X,A),2(1X,E14.8,1X,A),1X,F8.1)') &
684 28 : 'BASOPT| ', iset, "|", ibasis, "|", f_vec(icalc), "|", cond_vec(icalc), "|", tot_time(icalc)
685 : END DO
686 : END DO
687 5 : WRITE (unit_nr, '(1X,A)') "BASOPT| -----------------------------------------------------------------------"
688 5 : WRITE (unit_nr, '(1X,A,E14.8)') "BASOPT| Total residuum value: ", opt_bas%powell_param%f
689 5 : WRITE (unit_nr, '(A)') ""
690 : END IF
691 114 : CALL timestop(handle)
692 114 : END SUBROUTINE output_opt_info
693 :
694 : END MODULE optimize_basis
695 :
|
