Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 22 : 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()) &
200 57 : CALL powell_optimize(opt_bas%powell_param%nvar, opt_bas%x_opt, opt_bas%powell_param)
201 114 : CALL para_env_top%bcast(opt_bas%powell_param%state)
202 114 : CALL para_env_top%bcast(opt_bas%x_opt)
203 114 : CALL update_free_vars(opt_bas)
204 114 : write_basis = MOD(iopt, opt_bas%write_frequency) == 0
205 : CALL update_derived_basis_sets(opt_bas, write_basis, opt_bas%output_basis_file, &
206 114 : para_env_top)
207 118 : IF (opt_bas%powell_param%state == -1) EXIT
208 : END DO
209 :
210 : ! ======= Update the basis set and print the final basis =======
211 4 : IF (para_env_top%is_source()) THEN
212 2 : opt_bas%powell_param%state = 8
213 2 : CALL powell_optimize(opt_bas%powell_param%nvar, opt_bas%x_opt, opt_bas%powell_param)
214 : END IF
215 :
216 4 : CALL para_env_top%bcast(opt_bas%x_opt)
217 4 : CALL update_free_vars(opt_bas)
218 : CALL update_derived_basis_sets(opt_bas, .TRUE., opt_bas%output_basis_file, &
219 4 : para_env_top)
220 :
221 : ! ====== get rid of the f_env again =====
222 :
223 13 : DO icalc = SIZE(opt_bas%comp_group(mp_id)%member_list), 1, -1
224 9 : CALL f_env_get_from_id(f_env_id(icalc), f_env)
225 13 : CALL destroy_force_env(f_env_id(icalc), stat)
226 : END DO
227 4 : DEALLOCATE (f_env_id); DEALLOCATE (tot_time)
228 4 : CALL cp_fm_release(matrix_s_inv)
229 4 : CALL mp_para_env_release(para_env)
230 4 : CALL timestop(handle)
231 :
232 8 : END SUBROUTINE driver_optimization_para_low
233 :
234 : ! **************************************************************************************************
235 : !> \brief compute all ingredients for powell optimizer. Rho_diff,
236 : !> condition number, energy,... for all ttraining sets in
237 : !> the computational group
238 : !> \param opt_bas ...
239 : !> \param f_env_id ...
240 : !> \param matrix_S_inv ...
241 : !> \param tot_time ...
242 : !> \param para_env_top ...
243 : !> \param para_env ...
244 : !> \param iopt ...
245 : ! **************************************************************************************************
246 :
247 114 : SUBROUTINE compute_residuum_vectors(opt_bas, f_env_id, matrix_S_inv, tot_time, &
248 : para_env_top, para_env, iopt)
249 : TYPE(basis_optimization_type) :: opt_bas
250 : INTEGER, ALLOCATABLE, DIMENSION(:) :: f_env_id
251 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: matrix_S_inv
252 : REAL(KIND=dp), DIMENSION(:) :: tot_time
253 : TYPE(mp_para_env_type), POINTER :: para_env_top, para_env
254 : INTEGER :: iopt
255 :
256 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_residuum_vectors'
257 :
258 : CHARACTER(len=8) :: basis_type
259 : INTEGER :: bas_id, handle, icalc, icomb, ispin, &
260 : mp_id, my_id, nao, ncalc, nelectron, &
261 : nmo, nspins, set_id
262 : REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f
263 114 : REAL(KIND=dp), DIMENSION(:), POINTER :: cond_vec, energy, f_vec, my_time, &
264 114 : start_time
265 114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gdata
266 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
267 : TYPE(cp_fm_type), POINTER :: mo_coeff
268 114 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s_aux, matrix_s_aux_orb
269 : TYPE(f_env_type), POINTER :: f_env
270 : TYPE(force_env_type), POINTER :: force_env
271 114 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_aux
272 114 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
273 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
274 114 : POINTER :: sab_aux, sab_aux_orb
275 : TYPE(qs_environment_type), POINTER :: qs_env
276 114 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
277 : TYPE(qs_ks_env_type), POINTER :: ks_env
278 :
279 114 : CALL timeset(routineN, handle)
280 :
281 114 : basis_type = "AUX_OPT"
282 : !
283 114 : ncalc = opt_bas%ncombinations*opt_bas%ntraining_sets
284 342 : ALLOCATE (gdata(ncalc, 4))
285 114 : f_vec => gdata(:, 1)
286 114 : my_time => gdata(:, 2)
287 114 : cond_vec => gdata(:, 3)
288 114 : energy => gdata(:, 4)
289 : !
290 1986 : f_vec = 0.0_dp; cond_vec = 0.0_dp; my_time = 0.0_dp; energy = 0.0_dp
291 114 : mp_id = opt_bas%opt_id
292 342 : ALLOCATE (start_time(SIZE(opt_bas%comp_group(mp_id)%member_list)))
293 : !
294 348 : DO icalc = 1, SIZE(opt_bas%comp_group(mp_id)%member_list)
295 234 : my_id = opt_bas%comp_group(mp_id)%member_list(icalc) + 1
296 : ! setup timings
297 234 : start_time(icalc) = m_walltime()
298 :
299 234 : NULLIFY (matrix_s_aux_orb, matrix_s_aux)
300 234 : CALL get_set_and_basis_id(opt_bas%comp_group(mp_id)%member_list(icalc), opt_bas, set_id, bas_id)
301 234 : CALL f_env_get_from_id(f_env_id(icalc), f_env)
302 234 : force_env => f_env%force_env
303 234 : CALL force_env_get(force_env, qs_env=qs_env)
304 234 : CALL get_qs_env(qs_env, ks_env=ks_env)
305 234 : CALL update_basis_set(opt_bas, bas_id, basis_type, qs_env)
306 234 : NULLIFY (sab_aux, sab_aux_orb)
307 234 : CALL optbas_build_neighborlist(qs_env, sab_aux, sab_aux_orb, basis_type)
308 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s_aux, &
309 : basis_type_a=basis_type, &
310 : basis_type_b=basis_type, &
311 234 : sab_nl=sab_aux)
312 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s_aux_orb, &
313 : basis_type_a=basis_type, &
314 : basis_type_b="ORB", &
315 234 : sab_nl=sab_aux_orb)
316 234 : CALL release_neighbor_list_sets(sab_aux)
317 234 : CALL release_neighbor_list_sets(sab_aux_orb)
318 234 : CALL get_qs_env(qs_env, mos=mos, matrix_ks=matrix_ks)
319 :
320 234 : nspins = SIZE(mos)
321 936 : ALLOCATE (mos_aux(nspins))
322 234 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
323 234 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao, basis_type=basis_type)
324 468 : DO ispin = 1, nspins
325 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, maxocc=maxocc, nelectron=nelectron, &
326 234 : n_el_f=n_el_f, nmo=nmo, flexible_electron_count=flexible_electron_count)
327 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nmo, &
328 : context=mo_coeff%matrix_struct%context, &
329 234 : para_env=mo_coeff%matrix_struct%para_env)
330 : CALL allocate_mo_set(mos_aux(ispin), nao, nmo, nelectron, &
331 234 : n_el_f, maxocc, flexible_electron_count)
332 234 : CALL init_mo_set(mo_set=mos_aux(ispin), fm_struct=fm_struct, name="MO_AUX")
333 702 : CALL cp_fm_struct_release(fm_struct)
334 : END DO
335 :
336 234 : CALL fit_mo_coeffs(matrix_s_aux, matrix_s_aux_orb, mos, mos_aux)
337 : CALL evaluate_optvals(mos, mos_aux, matrix_ks, matrix_s_aux_orb(1)%matrix, &
338 : matrix_s_aux(1)%matrix, matrix_S_inv(icalc), &
339 234 : f_vec(my_id), energy(my_id), cond_vec(my_id))
340 :
341 468 : DO ispin = 1, nspins
342 468 : CALL deallocate_mo_set(mos_aux(ispin))
343 : END DO
344 234 : DEALLOCATE (mos_aux)
345 234 : IF (ASSOCIATED(matrix_s_aux)) CALL dbcsr_deallocate_matrix_set(matrix_s_aux)
346 234 : IF (ASSOCIATED(matrix_s_aux_orb)) CALL dbcsr_deallocate_matrix_set(matrix_s_aux_orb)
347 :
348 582 : my_time(my_id) = m_walltime() - start_time(icalc)
349 : END DO
350 :
351 114 : DEALLOCATE (start_time)
352 :
353 114 : IF (.NOT. para_env%is_source()) THEN
354 0 : f_vec = 0.0_dp; cond_vec = 0.0_dp; my_time = 0.0_dp; energy = 0.0_dp
355 : END IF
356 : ! collect date from all subgroup ionodes on the main ionode
357 4770 : CALL para_env_top%sum(gdata)
358 :
359 114 : opt_bas%powell_param%f = 0.0_dp
360 114 : IF (para_env_top%is_source()) THEN
361 291 : DO icalc = 1, SIZE(f_vec)
362 234 : icomb = MOD(icalc - 1, opt_bas%ncombinations)
363 : opt_bas%powell_param%f = opt_bas%powell_param%f + &
364 234 : (f_vec(icalc) + energy(icalc))*opt_bas%fval_weight(icomb)
365 234 : IF (opt_bas%use_condition_number) &
366 : opt_bas%powell_param%f = opt_bas%powell_param%f + &
367 291 : LOG(cond_vec(icalc))*opt_bas%condition_weight(icomb)
368 : END DO
369 : ELSE
370 993 : f_vec = 0.0_dp; cond_vec = 0.0_dp; my_time = 0.0_dp; energy = 0.0_dp
371 : END IF
372 114 : CALL para_env_top%bcast(opt_bas%powell_param%f)
373 :
374 : ! output info if required
375 114 : CALL output_opt_info(f_vec, cond_vec, my_time, tot_time, opt_bas, iopt, para_env_top)
376 114 : DEALLOCATE (gdata)
377 :
378 114 : CALL para_env_top%sync()
379 :
380 114 : CALL timestop(handle)
381 :
382 228 : END SUBROUTINE compute_residuum_vectors
383 :
384 : ! **************************************************************************************************
385 : !> \brief create the force_envs for every input in the computational group
386 : !> \param opt_bas ...
387 : !> \param f_env_id ...
388 : !> \param input_declaration ...
389 : !> \param matrix_s_inv ...
390 : !> \param para_env ...
391 : !> \param mpi_comm_opt ...
392 : ! **************************************************************************************************
393 :
394 17 : SUBROUTINE init_training_force_envs(opt_bas, f_env_id, input_declaration, matrix_s_inv, para_env, mpi_comm_opt)
395 :
396 : TYPE(basis_optimization_type) :: opt_bas
397 : INTEGER, ALLOCATABLE, DIMENSION(:) :: f_env_id
398 : TYPE(section_type), POINTER :: input_declaration
399 : TYPE(cp_fm_type), DIMENSION(:), INTENT(OUT) :: matrix_S_inv
400 : TYPE(mp_para_env_type), POINTER :: para_env
401 : TYPE(mp_comm_type) :: mpi_comm_opt
402 :
403 : CHARACTER(len=*), PARAMETER :: routineN = 'init_training_force_envs'
404 :
405 : CHARACTER(len=default_path_length) :: main_dir
406 : INTEGER :: bas_id, handle, icalc, ierr, mp_id, &
407 : set_id, stat
408 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
409 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
410 : TYPE(f_env_type), POINTER :: f_env
411 : TYPE(force_env_type), POINTER :: force_env
412 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
413 4 : POINTER :: sab_orb
414 : TYPE(qs_environment_type), POINTER :: qs_env
415 : TYPE(qs_ks_env_type), POINTER :: ks_env
416 : TYPE(section_vals_type), POINTER :: input_file
417 :
418 4 : CALL timeset(routineN, handle)
419 :
420 4 : NULLIFY (matrix_s, blacs_env, ks_env)
421 :
422 4 : mp_id = opt_bas%opt_id
423 4 : CALL m_getcwd(main_dir)
424 :
425 : ! ======= Create f_env for all calculations in MPI group =======
426 13 : DO icalc = 1, SIZE(opt_bas%comp_group(mp_id)%member_list)
427 9 : NULLIFY (input_file)
428 : ! parse the input of the training sets
429 9 : CALL get_set_and_basis_id(opt_bas%comp_group(mp_id)%member_list(icalc), opt_bas, set_id, bas_id)
430 9 : CALL m_chdir(TRIM(opt_bas%training_dir(set_id)), ierr)
431 9 : IF (ierr /= 0) THEN
432 : CALL cp_abort(__LOCATION__, &
433 0 : "Could not change to directory <"//TRIM(opt_bas%training_dir(set_id))//">")
434 : END IF
435 : input_file => read_input(input_declaration, &
436 : opt_bas%training_input(set_id), &
437 : initial_variables=empty_initial_variables, &
438 9 : para_env=para_env)
439 :
440 9 : CALL modify_input_settings(opt_bas, bas_id, input_file)
441 : CALL create_force_env(f_env_id(icalc), &
442 : input_declaration=input_declaration, &
443 : input_path=opt_bas%training_input(set_id), &
444 : input=input_file, &
445 : output_path="scrap_information", &
446 : mpi_comm=mpi_comm_opt, &
447 9 : ierr=stat)
448 :
449 : ! some weirdness with the default stacks defaults have to be addded to get the
450 : ! correct default program name this causes trouble with the timer stack if kept
451 9 : CALL f_env_add_defaults(f_env_id(icalc), f_env)
452 9 : force_env => f_env%force_env
453 9 : CALL force_env_get(force_env, qs_env=qs_env)
454 9 : CALL allocate_mo_sets(qs_env)
455 9 : CALL f_env_rm_defaults(f_env, stat)
456 9 : CALL get_qs_env(qs_env, ks_env=ks_env)
457 : CALL build_qs_neighbor_lists(qs_env, para_env, molecular=.FALSE., &
458 9 : force_env_section=qs_env%input)
459 : CALL get_ks_env(ks_env, &
460 : matrix_s=matrix_s, &
461 9 : sab_orb=sab_orb)
462 : CALL build_overlap_matrix(ks_env, matrix_s=matrix_s, &
463 : matrix_name="OVERLAP", &
464 : basis_type_a="ORB", &
465 : basis_type_b="ORB", &
466 9 : sab_nl=sab_orb)
467 9 : CALL set_ks_env(ks_env, matrix_s=matrix_s)
468 9 : CALL get_qs_env(qs_env, matrix_s=matrix_s, blacs_env=blacs_env)
469 : CALL calculate_overlap_inverse(matrix_s(1)%matrix, matrix_s_inv(icalc), &
470 9 : para_env, blacs_env)
471 9 : CALL calculate_ks_matrix(qs_env)
472 :
473 9 : CALL section_vals_release(input_file)
474 :
475 9 : CALL qs_env_part_release(qs_env)
476 :
477 22 : CALL m_chdir(TRIM(ADJUSTL(main_dir)), ierr)
478 : END DO
479 :
480 4 : CALL timestop(handle)
481 :
482 4 : END SUBROUTINE init_training_force_envs
483 :
484 : ! **************************************************************************************************
485 : !> \brief variable update from the powell vector for all sets
486 : !> \param opt_bas ...
487 : !> \author Florian Schiffmann
488 : ! **************************************************************************************************
489 :
490 118 : SUBROUTINE update_free_vars(opt_bas)
491 : TYPE(basis_optimization_type) :: opt_bas
492 :
493 : CHARACTER(len=*), PARAMETER :: routineN = 'update_free_vars'
494 :
495 : INTEGER :: handle, ikind, iset, ix
496 :
497 118 : CALL timeset(routineN, handle)
498 118 : ix = 0
499 354 : DO ikind = 1, opt_bas%nkind
500 590 : DO iset = 1, opt_bas%kind_basis(ikind)%flex_basis(0)%nsets
501 472 : CALL update_subset_freevars(opt_bas%kind_basis(ikind)%flex_basis(0)%subset(iset), ix, opt_bas%x_opt)
502 : END DO
503 : END DO
504 118 : CALL timestop(handle)
505 :
506 118 : END SUBROUTINE update_free_vars
507 :
508 : ! **************************************************************************************************
509 : !> \brief low level update for the basis sets. Exponents are transformed according to constraint
510 : !> \param subset ...
511 : !> \param ix ...
512 : !> \param x ...
513 : !> \author Florian Schiffmann
514 : ! **************************************************************************************************
515 :
516 236 : SUBROUTINE update_subset_freevars(subset, ix, x)
517 : TYPE(subset_type) :: subset
518 : INTEGER :: ix
519 : REAL(KIND=dp), DIMENSION(:) :: x
520 :
521 : CHARACTER(len=*), PARAMETER :: routineN = 'update_subset_freevars'
522 :
523 : INTEGER :: handle, icon1, icon2, icont, iexp, il, &
524 : istart
525 : REAL(KIND=dp) :: fermi_f, gs_scale
526 :
527 236 : CALL timeset(routineN, handle)
528 1888 : DO iexp = 1, subset%nexp
529 1652 : IF (subset%opt_exps(iexp)) THEN
530 0 : ix = ix + 1
531 0 : subset%exps(iexp) = ABS(x(ix))
532 0 : IF (subset%exp_has_const(iexp)) THEN
533 : !use a fermi function to keep exponents in a given range around their initial value
534 0 : fermi_f = 1.0_dp/(EXP((x(ix) - 1.0_dp)/0.5_dp) + 1.0_dp)
535 : subset%exps(iexp) = (2.0_dp*fermi_f - 1.0_dp)*subset%exp_const(iexp)%var_fac*subset%exp_const(iexp)%init + &
536 0 : subset%exp_const(iexp)%init
537 : ELSE
538 :
539 : END IF
540 : END IF
541 10974 : DO icont = 1, subset%ncon_tot
542 10738 : IF (subset%opt_coeff(iexp, icont)) THEN
543 9086 : ix = ix + 1
544 9086 : subset%coeff(iexp, icont) = x(ix)
545 : END IF
546 : END DO
547 : END DO
548 :
549 : ! orthonormalize contraction coefficients using gram schmidt
550 236 : istart = 1
551 826 : DO il = 1, subset%nl
552 1298 : DO icon1 = istart, istart + subset%l(il) - 2
553 2360 : DO icon2 = icon1 + 1, istart + subset%l(il) - 1
554 : gs_scale = DOT_PRODUCT(subset%coeff(:, icon2), subset%coeff(:, icon1))/ &
555 15930 : DOT_PRODUCT(subset%coeff(:, icon1), subset%coeff(:, icon1))
556 9204 : subset%coeff(:, icon2) = subset%coeff(:, icon2) - gs_scale*subset%coeff(:, icon1)
557 : END DO
558 : END DO
559 826 : istart = istart + subset%l(il)
560 : END DO
561 :
562 1534 : DO icon1 = 1, subset%ncon_tot
563 : subset%coeff(:, icon1) = subset%coeff(:, icon1)/ &
564 19706 : SQRT(DOT_PRODUCT(subset%coeff(:, icon1), subset%coeff(:, icon1)))
565 : END DO
566 236 : CALL timestop(handle)
567 :
568 236 : END SUBROUTINE update_subset_freevars
569 :
570 : ! **************************************************************************************************
571 : !> \brief variable initialization for the powell vector for all sets
572 : !> \param opt_bas ...
573 : !> \author Florian Schiffmann
574 : ! **************************************************************************************************
575 :
576 4 : SUBROUTINE init_free_vars(opt_bas)
577 : TYPE(basis_optimization_type) :: opt_bas
578 :
579 : CHARACTER(len=*), PARAMETER :: routineN = 'init_free_vars'
580 :
581 : INTEGER :: handle, ikind, iset, ix
582 :
583 4 : CALL timeset(routineN, handle)
584 4 : ix = 0
585 12 : DO ikind = 1, opt_bas%nkind
586 20 : DO iset = 1, opt_bas%kind_basis(ikind)%flex_basis(0)%nsets
587 16 : CALL init_subset_freevars(opt_bas%kind_basis(ikind)%flex_basis(0)%subset(iset), ix, opt_bas%x_opt)
588 : END DO
589 : END DO
590 4 : CALL timestop(handle)
591 :
592 4 : END SUBROUTINE init_free_vars
593 :
594 : ! **************************************************************************************************
595 : !> \brief variable initialization for the powell vector from low level informations
596 : !> constraint exponents will be mapped on a fermi function
597 : !> \param subset ...
598 : !> \param ix ...
599 : !> \param x ...
600 : !> \author Florian Schiffmann
601 : ! **************************************************************************************************
602 :
603 8 : SUBROUTINE init_subset_freevars(subset, ix, x)
604 : TYPE(subset_type) :: subset
605 : INTEGER :: ix
606 : REAL(KIND=dp), DIMENSION(:) :: x
607 :
608 : CHARACTER(len=*), PARAMETER :: routineN = 'init_subset_freevars'
609 :
610 : INTEGER :: handle, icont, iexp
611 : REAL(KIND=dp) :: fract
612 :
613 8 : CALL timeset(routineN, handle)
614 :
615 64 : DO iexp = 1, subset%nexp
616 56 : IF (subset%opt_exps(iexp)) THEN
617 0 : ix = ix + 1
618 0 : x(ix) = subset%exps(iexp)
619 0 : IF (subset%exp_has_const(iexp)) THEN
620 0 : IF (subset%exp_const(iexp)%const_type == 0) THEN
621 : fract = 1.0_dp + (subset%exps(iexp) - subset%exp_const(iexp)%init)/ &
622 0 : (subset%exp_const(iexp)%init*subset%exp_const(iexp)%var_fac)
623 0 : x(ix) = 0.5_dp*LOG((2.0_dp/fract - 1.0_dp)) + 1.0_dp
624 : END IF
625 0 : IF (subset%exp_const(iexp)%const_type == 1) THEN
626 0 : x(ix) = 1.0_dp
627 : END IF
628 : END IF
629 : END IF
630 372 : DO icont = 1, subset%ncon_tot
631 364 : IF (subset%opt_coeff(iexp, icont)) THEN
632 308 : ix = ix + 1
633 308 : x(ix) = subset%coeff(iexp, icont)
634 : END IF
635 : END DO
636 : END DO
637 8 : CALL timestop(handle)
638 :
639 8 : END SUBROUTINE init_subset_freevars
640 :
641 : ! **************************************************************************************************
642 : !> \brief commuticates all info to the master and assembles the output
643 : !> \param f_vec ...
644 : !> \param cond_vec ...
645 : !> \param my_time ...
646 : !> \param tot_time ...
647 : !> \param opt_bas ...
648 : !> \param iopt ...
649 : !> \param para_env_top ...
650 : !> \author Florian Schiffmann
651 : ! **************************************************************************************************
652 :
653 114 : SUBROUTINE output_opt_info(f_vec, cond_vec, my_time, tot_time, opt_bas, iopt, para_env_top)
654 : REAL(KIND=dp), DIMENSION(:) :: f_vec, cond_vec, my_time, tot_time
655 : TYPE(basis_optimization_type) :: opt_bas
656 : INTEGER :: iopt
657 : TYPE(mp_para_env_type), POINTER :: para_env_top
658 :
659 : CHARACTER(len=*), PARAMETER :: routineN = 'output_opt_info'
660 :
661 : INTEGER :: handle, ibasis, icalc, iset, unit_nr
662 : TYPE(cp_logger_type), POINTER :: logger
663 :
664 114 : CALL timeset(routineN, handle)
665 114 : logger => cp_get_default_logger()
666 :
667 582 : tot_time = tot_time + my_time
668 :
669 114 : unit_nr = -1
670 114 : IF (para_env_top%is_source() .AND. (MOD(iopt, opt_bas%write_frequency) == 0 .OR. iopt == opt_bas%powell_param%maxfun)) &
671 5 : unit_nr = cp_logger_get_default_unit_nr(logger)
672 :
673 5 : IF (unit_nr .GT. 0) THEN
674 5 : WRITE (unit_nr, '(1X,A,I8)') "BASOPT| Information at iteration number:", iopt
675 5 : WRITE (unit_nr, '(1X,A)') "BASOPT| Training set | Combination | Rho difference | Condition num. | Time"
676 5 : WRITE (unit_nr, '(1X,A)') "BASOPT| -----------------------------------------------------------------------"
677 5 : icalc = 0
678 12 : DO iset = 1, opt_bas%ntraining_sets
679 33 : DO ibasis = 1, opt_bas%ncombinations
680 21 : icalc = icalc + 1
681 : WRITE (unit_nr, '(1X,A,2(5X,I3,5X,A),2(1X,E14.8,1X,A),1X,F8.1)') &
682 28 : 'BASOPT| ', iset, "|", ibasis, "|", f_vec(icalc), "|", cond_vec(icalc), "|", tot_time(icalc)
683 : END DO
684 : END DO
685 5 : WRITE (unit_nr, '(1X,A)') "BASOPT| -----------------------------------------------------------------------"
686 5 : WRITE (unit_nr, '(1X,A,E14.8)') "BASOPT| Total residuum value: ", opt_bas%powell_param%f
687 5 : WRITE (unit_nr, '(A)') ""
688 : END IF
689 114 : CALL timestop(handle)
690 114 : END SUBROUTINE output_opt_info
691 :
692 : END MODULE optimize_basis
693 :
|
