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_utils
8 : USE cp_files, ONLY: close_file,&
9 : open_file
10 : USE cp_log_handling, ONLY: cp_get_default_logger,&
11 : cp_logger_get_default_unit_nr,&
12 : cp_logger_type,&
13 : cp_to_string
14 : USE cp_parser_methods, ONLY: parser_get_object,&
15 : parser_search_string
16 : USE cp_parser_types, ONLY: cp_parser_type,&
17 : parser_create,&
18 : parser_release
19 : USE input_constants, ONLY: do_opt_all,&
20 : do_opt_coeff,&
21 : do_opt_exps,&
22 : do_opt_none
23 : USE input_section_types, ONLY: section_vals_get,&
24 : section_vals_get_subs_vals,&
25 : section_vals_type,&
26 : section_vals_val_get
27 : USE kinds, ONLY: default_path_length,&
28 : default_string_length,&
29 : dp
30 : USE machine, ONLY: default_output_unit,&
31 : m_getcwd
32 : USE message_passing, ONLY: mp_para_env_type
33 : USE optimize_basis_types, ONLY: basis_optimization_type,&
34 : derived_basis_info,&
35 : flex_basis_type,&
36 : method_mo_fit_occ,&
37 : method_mo_fit_occ_virtual,&
38 : subset_type
39 : USE powell, ONLY: opt_state_type
40 : USE string_utilities, ONLY: uppercase
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 : PRIVATE
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_basis_utils'
47 :
48 : PUBLIC :: optimize_basis_init_read_input, get_set_and_basis_id, &
49 : update_derived_basis_sets, write_basis
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief initialize all parts of the optimization type and read input settings
55 : !> \param opt_bas ...
56 : !> \param root_section ...
57 : !> \param para_env ...
58 : !> \param quiet_output ...
59 : !> \param embedded_training ...
60 : !> \author Florian Schiffmann
61 : ! **************************************************************************************************
62 :
63 8 : SUBROUTINE optimize_basis_init_read_input(opt_bas, root_section, para_env, quiet_output, embedded_training)
64 : TYPE(basis_optimization_type) :: opt_bas
65 : TYPE(section_vals_type), POINTER :: root_section
66 : TYPE(mp_para_env_type), POINTER :: para_env
67 : LOGICAL, INTENT(IN), OPTIONAL :: quiet_output, embedded_training
68 :
69 : CHARACTER(LEN=default_path_length) :: main_dir
70 : INTEGER :: ikind, iset, iweight, nrep
71 : LOGICAL :: frontier_orbitals_explicit, &
72 : use_embedded_training
73 : TYPE(section_vals_type), POINTER :: frontier_orbitals_section, kind_section, &
74 : optbas_section, powell_section, &
75 : train_section
76 :
77 8 : optbas_section => section_vals_get_subs_vals(root_section, "OPTIMIZE_BASIS")
78 8 : powell_section => section_vals_get_subs_vals(optbas_section, "OPTIMIZATION")
79 8 : train_section => section_vals_get_subs_vals(optbas_section, "TRAINING_FILES")
80 8 : kind_section => section_vals_get_subs_vals(optbas_section, "FIT_KIND")
81 8 : frontier_orbitals_section => section_vals_get_subs_vals(optbas_section, "FRONTIER_ORBITALS")
82 :
83 8 : IF (PRESENT(quiet_output)) opt_bas%quiet_output = quiet_output
84 8 : use_embedded_training = .FALSE.
85 8 : IF (PRESENT(embedded_training)) use_embedded_training = embedded_training
86 :
87 8 : CALL section_vals_val_get(optbas_section, "BASIS_TEMPLATE_FILE", c_val=opt_bas%template_basis_file)
88 8 : CALL section_vals_val_get(optbas_section, "BASIS_OUTPUT_FILE", c_val=opt_bas%output_basis_file)
89 8 : CALL m_getcwd(main_dir)
90 :
91 8 : CALL section_vals_val_get(optbas_section, "WRITE_FREQUENCY", i_val=opt_bas%write_frequency)
92 8 : CALL section_vals_val_get(optbas_section, "USE_CONDITION_NUMBER", l_val=opt_bas%use_condition_number)
93 :
94 8 : CALL generate_initial_basis(kind_section, opt_bas, para_env)
95 :
96 8 : CALL section_vals_get(frontier_orbitals_section, explicit=frontier_orbitals_explicit)
97 8 : IF (use_embedded_training) THEN
98 2 : opt_bas%method = method_mo_fit_occ_virtual
99 6 : ELSE IF (frontier_orbitals_explicit) THEN
100 2 : opt_bas%method = method_mo_fit_occ_virtual
101 : CALL section_vals_val_get(frontier_orbitals_section, "OCCUPIED_WEIGHT", &
102 2 : r_val=opt_bas%occupied_weight)
103 : CALL section_vals_val_get(frontier_orbitals_section, "VIRTUAL_WEIGHT", &
104 2 : r_val=opt_bas%virtual_weight)
105 : CALL section_vals_val_get(frontier_orbitals_section, "EMPTY_OVERLAP_WEIGHT", &
106 2 : r_val=opt_bas%empty_overlap_weight)
107 : CALL section_vals_val_get(frontier_orbitals_section, "GAP_WEIGHT", &
108 2 : r_val=opt_bas%gap_weight)
109 : CALL section_vals_val_get(frontier_orbitals_section, "COEFFICIENT_WEIGHT", &
110 2 : r_val=opt_bas%coefficient_weight)
111 : CALL section_vals_val_get(frontier_orbitals_section, "GAP_ENERGY_SCALE", &
112 2 : r_val=opt_bas%gap_energy_scale)
113 : CALL section_vals_val_get(frontier_orbitals_section, "VIRTUAL_ENERGY_CUTOFF", &
114 2 : r_val=opt_bas%virtual_energy_cutoff)
115 : CALL section_vals_val_get(frontier_orbitals_section, "VIRTUAL_ENERGY_SMOOTHING", &
116 2 : r_val=opt_bas%virtual_energy_smoothing)
117 : END IF
118 8 : IF (opt_bas%method == method_mo_fit_occ) THEN
119 4 : CALL section_vals_val_get(optbas_section, "BASIS_WORK_FILE", c_val=opt_bas%work_basis_file)
120 4 : opt_bas%work_basis_file = TRIM(ADJUSTL(main_dir))//"/"//TRIM(ADJUSTL(opt_bas%work_basis_file))
121 : END IF
122 8 : IF (frontier_orbitals_explicit) THEN
123 2 : IF (opt_bas%gap_energy_scale <= 0.0_dp) CPABORT("GAP_ENERGY_SCALE must be positive")
124 2 : IF (opt_bas%virtual_energy_smoothing <= 0.0_dp) THEN
125 0 : CPABORT("VIRTUAL_ENERGY_SMOOTHING must be positive")
126 : END IF
127 2 : IF (MIN(opt_bas%occupied_weight, opt_bas%virtual_weight, &
128 : opt_bas%empty_overlap_weight, opt_bas%gap_weight, &
129 : opt_bas%coefficient_weight) < 0.0_dp) THEN
130 0 : CPABORT("FRONTIER_ORBITALS weights must be non-negative")
131 : END IF
132 : END IF
133 8 : IF (opt_bas%method == method_mo_fit_occ_virtual) THEN
134 12 : DO ikind = 1, opt_bas%nkind
135 20 : DO iset = 1, opt_bas%kind_basis(ikind)%flex_basis(0)%nsets
136 68 : IF (ANY(opt_bas%kind_basis(ikind)%flex_basis(0)%subset(iset)%opt_exps)) THEN
137 0 : CPABORT("FRONTIER_ORBITALS requires fixed Gaussian exponents")
138 : END IF
139 : END DO
140 : END DO
141 : END IF
142 :
143 8 : CALL section_vals_get(train_section, n_repetition=opt_bas%ntraining_sets)
144 8 : IF (use_embedded_training) THEN
145 2 : IF (opt_bas%ntraining_sets /= 0) THEN
146 : CALL cp_abort(__LOCATION__, &
147 : "TRAINING_FILES is not allowed together with FRONTIER_ORBITAL_SCREENING "// &
148 0 : "because it is overwritten automatically")
149 : END IF
150 2 : opt_bas%ntraining_sets = 1
151 6 : ELSE IF (opt_bas%ntraining_sets == 0) THEN
152 0 : CPABORT("No training set was specified in the Input")
153 : END IF
154 :
155 24 : ALLOCATE (opt_bas%training_input(opt_bas%ntraining_sets))
156 16 : ALLOCATE (opt_bas%training_dir(opt_bas%ntraining_sets))
157 8 : IF (use_embedded_training) THEN
158 2 : opt_bas%training_input(1) = ""
159 2 : opt_bas%training_dir(1) = "."
160 : ELSE
161 14 : DO iset = 1, opt_bas%ntraining_sets
162 : CALL section_vals_val_get(train_section, "DIRECTORY", c_val=opt_bas%training_dir(iset), &
163 8 : i_rep_section=iset)
164 : CALL section_vals_val_get(train_section, "INPUT_FILE_NAME", c_val=opt_bas%training_input(iset), &
165 14 : i_rep_section=iset)
166 : END DO
167 : END IF
168 :
169 8 : CALL init_powell_var(opt_bas%powell_param, powell_section)
170 8 : opt_bas%powell_param%nvar = SIZE(opt_bas%x_opt)
171 :
172 : CALL generate_derived_basis_sets( &
173 8 : opt_bas, para_env, write_work_file=opt_bas%method == method_mo_fit_occ)
174 :
175 8 : CALL generate_basis_combinations(opt_bas, optbas_section)
176 :
177 8 : CALL section_vals_val_get(optbas_section, "RESIDUUM_WEIGHT", n_rep_val=nrep)
178 24 : ALLOCATE (opt_bas%fval_weight(0:opt_bas%ncombinations))
179 32 : opt_bas%fval_weight = 1.0_dp
180 24 : DO iweight = 1, nrep
181 : CALL section_vals_val_get(optbas_section, "RESIDUUM_WEIGHT", r_val=opt_bas%fval_weight(iweight - 1), &
182 24 : i_rep_val=iweight)
183 : END DO
184 :
185 8 : CALL section_vals_val_get(optbas_section, "CONDITION_WEIGHT", n_rep_val=nrep)
186 24 : ALLOCATE (opt_bas%condition_weight(0:opt_bas%ncombinations))
187 32 : opt_bas%condition_weight = 1.0_dp
188 24 : DO iweight = 1, nrep
189 : CALL section_vals_val_get(optbas_section, "CONDITION_WEIGHT", r_val=opt_bas%condition_weight(iweight - 1), &
190 24 : i_rep_val=iweight)
191 : END DO
192 :
193 8 : CALL generate_computation_groups(opt_bas, optbas_section, para_env)
194 :
195 8 : IF (.NOT. opt_bas%quiet_output) CALL print_opt_info(opt_bas)
196 :
197 24 : END SUBROUTINE optimize_basis_init_read_input
198 :
199 : ! **************************************************************************************************
200 : !> \brief ...
201 : !> \param opt_bas ...
202 : ! **************************************************************************************************
203 6 : SUBROUTINE print_opt_info(opt_bas)
204 : TYPE(basis_optimization_type) :: opt_bas
205 :
206 : INTEGER :: icomb, ikind, unit_nr
207 : TYPE(cp_logger_type), POINTER :: logger
208 :
209 6 : logger => cp_get_default_logger()
210 6 : unit_nr = -1
211 6 : IF (logger%para_env%is_source()) THEN
212 3 : unit_nr = cp_logger_get_default_unit_nr(logger)
213 : END IF
214 :
215 3 : IF (unit_nr > 0) THEN
216 3 : WRITE (unit_nr, '(1X,A,A)') "BASOPT| Total number of calculations ", &
217 6 : TRIM(cp_to_string(opt_bas%ncombinations*opt_bas%ntraining_sets))
218 3 : WRITE (unit_nr, '(A)') ""
219 10 : DO icomb = 1, opt_bas%ncombinations
220 7 : WRITE (unit_nr, '(1X,A,A)') "BASOPT| Content of basis combination ", TRIM(cp_to_string(icomb))
221 21 : DO ikind = 1, opt_bas%nkind
222 14 : WRITE (unit_nr, '(1X,A,A4,4X,A,3X,A)') "BASOPT| Element: ", TRIM(opt_bas%kind_basis(ikind)%element), &
223 35 : "Basis set: ", TRIM(opt_bas%kind_basis(ikind)%flex_basis(opt_bas%combination(icomb, ikind))%basis_name)
224 : END DO
225 10 : WRITE (unit_nr, '(A)') ""
226 : END DO
227 : END IF
228 6 : END SUBROUTINE print_opt_info
229 :
230 : ! **************************************************************************************************
231 : !> \brief Generation of the requested basis set combinations if multiple kinds
232 : !> are fitted at the same time (if not specified create all possible)
233 : !> \param opt_bas ...
234 : !> \param optbas_section ...
235 : !> \author Florian Schiffmann
236 : ! **************************************************************************************************
237 8 : SUBROUTINE generate_basis_combinations(opt_bas, optbas_section)
238 : TYPE(basis_optimization_type) :: opt_bas
239 : TYPE(section_vals_type), POINTER :: optbas_section
240 :
241 : INTEGER :: i, ikind, j, n_rep
242 8 : INTEGER, DIMENSION(:), POINTER :: i_vals, tmp_i, tmp_i2
243 : LOGICAL :: explicit, raise
244 :
245 : !setup the basis combinations to optimize
246 :
247 8 : CALL section_vals_val_get(optbas_section, "BASIS_COMBINATIONS", explicit=explicit, n_rep_val=n_rep)
248 8 : IF (.NOT. explicit) THEN
249 4 : opt_bas%ncombinations = 1
250 12 : ALLOCATE (tmp_i(opt_bas%nkind))
251 8 : ALLOCATE (tmp_i2(opt_bas%nkind))
252 12 : DO ikind = 1, opt_bas%nkind
253 8 : opt_bas%ncombinations = opt_bas%ncombinations*SIZE(opt_bas%kind_basis(ikind)%flex_basis)
254 12 : tmp_i(ikind) = opt_bas%kind_basis(ikind)%nbasis_deriv
255 : END DO
256 16 : ALLOCATE (opt_bas%combination(opt_bas%ncombinations, opt_bas%nkind))
257 12 : tmp_i2 = 0
258 8 : DO i = 1, opt_bas%ncombinations
259 12 : DO j = 1, opt_bas%nkind
260 12 : opt_bas%combination(i, j) = tmp_i2(j)
261 : END DO
262 4 : tmp_i2(opt_bas%nkind) = tmp_i2(opt_bas%nkind) + 1
263 4 : raise = .FALSE.
264 16 : DO j = opt_bas%nkind, 1, -1
265 8 : IF (raise) tmp_i2(j) = tmp_i2(j) + 1
266 12 : IF (tmp_i2(j) > tmp_i(j)) THEN
267 8 : tmp_i2(j) = 0
268 8 : raise = .TRUE.
269 : END IF
270 : END DO
271 : END DO
272 4 : DEALLOCATE (tmp_i)
273 4 : DEALLOCATE (tmp_i2)
274 : ELSE
275 4 : opt_bas%ncombinations = n_rep
276 16 : ALLOCATE (opt_bas%combination(opt_bas%ncombinations, opt_bas%nkind))
277 16 : DO i = 1, n_rep
278 12 : CALL section_vals_val_get(optbas_section, "BASIS_COMBINATIONS", i_vals=i_vals, i_rep_val=i)
279 40 : opt_bas%combination(i, :) = i_vals(:)
280 : END DO
281 : END IF
282 :
283 8 : END SUBROUTINE generate_basis_combinations
284 :
285 : ! **************************************************************************************************
286 : !> \brief returns a mapping from the calculation id to the trainings set id and
287 : !> basis combination id
288 : !> \param calc_id ...
289 : !> \param opt_bas ...
290 : !> \param set_id ...
291 : !> \param bas_id ...
292 : !> \author Florian Schiffmann
293 : ! **************************************************************************************************
294 :
295 656 : SUBROUTINE get_set_and_basis_id(calc_id, opt_bas, set_id, bas_id)
296 :
297 : INTEGER :: calc_id
298 : TYPE(basis_optimization_type) :: opt_bas
299 : INTEGER :: set_id, bas_id
300 :
301 : INTEGER :: ncom, nset
302 :
303 656 : ncom = opt_bas%ncombinations
304 656 : nset = opt_bas%ntraining_sets
305 :
306 656 : set_id = (calc_id)/ncom + 1
307 656 : bas_id = MOD(calc_id, ncom) + 1
308 :
309 656 : END SUBROUTINE get_set_and_basis_id
310 :
311 : ! **************************************************************************************************
312 : !> \brief Pack calculations in groups. If less MPI tasks than systems are used
313 : !> multiple systems will be assigned to a single MPI task
314 : !> \param opt_bas ...
315 : !> \param optbas_section ...
316 : !> \param para_env ...
317 : !> \author Florian Schiffmann
318 : ! **************************************************************************************************
319 :
320 8 : SUBROUTINE generate_computation_groups(opt_bas, optbas_section, para_env)
321 : TYPE(basis_optimization_type) :: opt_bas
322 : TYPE(section_vals_type), POINTER :: optbas_section
323 : TYPE(mp_para_env_type), POINTER :: para_env
324 :
325 : INTEGER :: iadd1, iadd2, icount, igroup, isize, j, &
326 : ncalc, nproc, nptot
327 8 : INTEGER, DIMENSION(:), POINTER :: i_vals
328 : LOGICAL :: explicit
329 :
330 8 : nproc = para_env%num_pe
331 8 : ncalc = opt_bas%ncombinations*opt_bas%ntraining_sets
332 8 : CALL section_vals_val_get(optbas_section, "GROUP_PARTITION", explicit=explicit)
333 :
334 : ! No input information available, try to equally distribute
335 8 : IF (.NOT. explicit) THEN
336 8 : IF (nproc >= ncalc) THEN
337 4 : iadd1 = nproc/ncalc
338 4 : iadd2 = MOD(nproc, ncalc)
339 16 : ALLOCATE (opt_bas%comp_group(ncalc))
340 12 : ALLOCATE (opt_bas%group_partition(0:ncalc - 1))
341 8 : DO igroup = 0, ncalc - 1
342 4 : ALLOCATE (opt_bas%comp_group(igroup + 1)%member_list(1))
343 4 : opt_bas%comp_group(igroup + 1)%member_list(1) = igroup
344 4 : opt_bas%group_partition(igroup) = iadd1
345 8 : IF (igroup < iadd2) opt_bas%group_partition(igroup) = opt_bas%group_partition(igroup) + 1
346 : END DO
347 : ELSE
348 4 : iadd1 = ncalc/nproc
349 4 : iadd2 = MOD(ncalc, nproc)
350 20 : ALLOCATE (opt_bas%comp_group(nproc))
351 12 : ALLOCATE (opt_bas%group_partition(0:nproc - 1))
352 4 : icount = 0
353 12 : DO igroup = 0, nproc - 1
354 8 : opt_bas%group_partition(igroup) = 1
355 8 : isize = iadd1
356 8 : IF (igroup < iadd2) isize = isize + 1
357 24 : ALLOCATE (opt_bas%comp_group(igroup + 1)%member_list(isize))
358 30 : DO j = 1, isize
359 18 : opt_bas%comp_group(igroup + 1)%member_list(j) = icount
360 26 : icount = icount + 1
361 : END DO
362 : END DO
363 : END IF
364 : ELSE
365 :
366 : ! Group partition from input. see if all systems can be assigned. If not add to existing group
367 0 : CALL section_vals_val_get(optbas_section, "GROUP_PARTITION", i_vals=i_vals)
368 0 : isize = SIZE(i_vals)
369 0 : nptot = SUM(i_vals)
370 0 : IF (nptot /= nproc) THEN
371 : CALL cp_abort(__LOCATION__, &
372 : "Number of processors in group distribution does not match number of MPI tasks."// &
373 0 : " Please change input.")
374 : END IF
375 0 : IF (.NOT. isize <= ncalc) THEN
376 : CALL cp_abort(__LOCATION__, &
377 : "Number of Groups larger than number of calculations"// &
378 0 : " Please change input.")
379 : END IF
380 0 : CPASSERT(nptot == nproc)
381 0 : ALLOCATE (opt_bas%comp_group(isize))
382 0 : ALLOCATE (opt_bas%group_partition(0:isize - 1))
383 0 : IF (isize < ncalc) THEN
384 0 : iadd1 = ncalc/isize
385 0 : iadd2 = MOD(ncalc, isize)
386 0 : icount = 0
387 0 : DO igroup = 0, isize - 1
388 0 : opt_bas%group_partition(igroup) = i_vals(igroup + 1)
389 0 : isize = iadd1
390 0 : IF (igroup < iadd2) isize = isize + 1
391 0 : ALLOCATE (opt_bas%comp_group(igroup + 1)%member_list(isize))
392 0 : DO j = 1, isize
393 0 : opt_bas%comp_group(igroup + 1)%member_list(j) = icount
394 0 : icount = icount + 1
395 : END DO
396 : END DO
397 : ELSE
398 0 : DO igroup = 0, isize - 1
399 0 : opt_bas%group_partition(igroup) = i_vals(igroup + 1)
400 0 : ALLOCATE (opt_bas%comp_group(igroup + 1)%member_list(1))
401 0 : opt_bas%comp_group(igroup + 1)%member_list(1) = igroup
402 : END DO
403 : END IF
404 : END IF
405 :
406 8 : END SUBROUTINE generate_computation_groups
407 :
408 : ! **************************************************************************************************
409 : !> \brief Regenerate the basis sets from reference 0 after an update from the
410 : !> optimizer to reference was performed, and print basis file if required
411 : !> \param opt_bas ...
412 : !> \param write_it ...
413 : !> \param output_file ...
414 : !> \param para_env ...
415 : !> \author Florian Schiffmann
416 : ! **************************************************************************************************
417 526 : SUBROUTINE update_derived_basis_sets(opt_bas, write_it, output_file, para_env)
418 : TYPE(basis_optimization_type) :: opt_bas
419 : LOGICAL :: write_it
420 : CHARACTER(LEN=default_path_length) :: output_file
421 : TYPE(mp_para_env_type), POINTER :: para_env
422 :
423 : INTEGER :: ibasis, ikind, unit_nr
424 :
425 1578 : DO ikind = 1, opt_bas%nkind
426 2050 : DO ibasis = 1, opt_bas%kind_basis(ikind)%nbasis_deriv
427 : CALL update_used_parts(opt_bas%kind_basis(ikind)%deriv_info(ibasis), &
428 : opt_bas%kind_basis(ikind)%flex_basis(0), &
429 1524 : opt_bas%kind_basis(ikind)%flex_basis(ibasis))
430 : END DO
431 : END DO
432 :
433 526 : IF (write_it) THEN
434 14 : IF (para_env%is_source()) THEN
435 : CALL open_file(file_name=output_file, file_status="UNKNOWN", &
436 7 : file_action="WRITE", unit_number=unit_nr)
437 : ELSE
438 7 : unit_nr = -999
439 : END IF
440 42 : DO ikind = 1, opt_bas%nkind
441 118 : DO ibasis = 0, opt_bas%kind_basis(ikind)%nbasis_deriv
442 : CALL write_basis(opt_bas%kind_basis(ikind)%flex_basis(ibasis), opt_bas%kind_basis(ikind)%element, &
443 104 : unit_nr)
444 : END DO
445 : END DO
446 14 : IF (para_env%is_source()) CALL close_file(unit_number=unit_nr)
447 : END IF
448 :
449 526 : END SUBROUTINE update_derived_basis_sets
450 :
451 : ! **************************************************************************************************
452 : !> \brief Update the the information in a given basis set
453 : !> \param info_new ...
454 : !> \param basis ...
455 : !> \param basis_new ...
456 : !> \author Florian Schiffmann
457 : ! **************************************************************************************************
458 :
459 472 : SUBROUTINE update_used_parts(info_new, basis, basis_new)
460 : TYPE(derived_basis_info) :: info_new
461 : TYPE(flex_basis_type) :: basis, basis_new
462 :
463 : INTEGER :: icont, iset, jcont, jset
464 :
465 472 : jset = 0
466 944 : DO iset = 1, basis%nsets
467 944 : IF (info_new%in_use_set(iset)) THEN
468 472 : jset = jset + 1
469 3776 : basis_new%subset(jset)%exps(:) = basis%subset(iset)%exps
470 472 : jcont = 0
471 3068 : DO icont = 1, basis%subset(iset)%ncon_tot
472 3068 : IF (info_new%use_contr(iset)%in_use(icont)) THEN
473 1298 : jcont = jcont + 1
474 10384 : basis_new%subset(jset)%coeff(:, jcont) = basis%subset(iset)%coeff(:, icont)
475 : END IF
476 : END DO
477 : END IF
478 : END DO
479 :
480 472 : END SUBROUTINE update_used_parts
481 :
482 : ! **************************************************************************************************
483 : !> \brief Initial generation of the basis set from the file and DERIVED_SET
484 : !> \param opt_bas ...
485 : !> \param para_env ...
486 : !> \param write_work_file ...
487 : !> \author Florian Schiffmann
488 : ! **************************************************************************************************
489 :
490 8 : SUBROUTINE generate_derived_basis_sets(opt_bas, para_env, write_work_file)
491 : TYPE(basis_optimization_type) :: opt_bas
492 : TYPE(mp_para_env_type), POINTER :: para_env
493 : LOGICAL, INTENT(IN) :: write_work_file
494 :
495 : INTEGER :: ibasis, ikind, iref, jbasis, unit_nr
496 :
497 24 : DO ikind = 1, opt_bas%nkind
498 16 : CALL init_deriv_info_ref(opt_bas%kind_basis(ikind)%deriv_info(0), opt_bas%kind_basis(ikind)%flex_basis(0))
499 16 : opt_bas%kind_basis(ikind)%deriv_info(0)%basis_name = TRIM(ADJUSTL(opt_bas%kind_basis(ikind)%basis_name))
500 : ! initialize the reference set used as template for the rest
501 40 : DO ibasis = 1, opt_bas%kind_basis(ikind)%nbasis_deriv
502 16 : iref = opt_bas%kind_basis(ikind)%deriv_info(ibasis)%reference_set
503 80 : DO jbasis = 0, opt_bas%kind_basis(ikind)%nbasis_deriv
504 48 : IF (iref == jbasis) CALL setup_used_parts_init_basis(opt_bas%kind_basis(ikind)%deriv_info(ibasis), &
505 : opt_bas%kind_basis(ikind)%deriv_info(iref), &
506 : opt_bas%kind_basis(ikind)%flex_basis(0), &
507 32 : opt_bas%kind_basis(ikind)%flex_basis(ibasis))
508 : END DO
509 : END DO
510 : END DO
511 :
512 8 : IF (write_work_file .AND. para_env%is_source()) THEN
513 : CALL open_file(file_name=opt_bas%work_basis_file, file_status="UNKNOWN", &
514 2 : file_action="WRITE", unit_number=unit_nr)
515 : ELSE
516 6 : unit_nr = -999
517 : END IF
518 24 : DO ikind = 1, opt_bas%nkind
519 56 : DO ibasis = 0, opt_bas%kind_basis(ikind)%nbasis_deriv
520 32 : IF (LEN_TRIM(opt_bas%kind_basis(ikind)%deriv_info(ibasis)%basis_name) > 0) THEN
521 : opt_bas%kind_basis(ikind)%flex_basis(ibasis)%basis_name = &
522 24 : TRIM(ADJUSTL(opt_bas%kind_basis(ikind)%deriv_info(ibasis)%basis_name))
523 : ELSE
524 : opt_bas%kind_basis(ikind)%flex_basis(ibasis)%basis_name = &
525 8 : TRIM(ADJUSTL(opt_bas%kind_basis(ikind)%basis_name))//"-DERIVED_SET-"//ADJUSTL(cp_to_string(ibasis))
526 : END IF
527 48 : IF (write_work_file) THEN
528 : CALL write_basis(opt_bas%kind_basis(ikind)%flex_basis(ibasis), &
529 24 : opt_bas%kind_basis(ikind)%element, unit_nr)
530 : END IF
531 : END DO
532 : END DO
533 8 : IF (write_work_file .AND. para_env%is_source()) CALL close_file(unit_number=unit_nr)
534 :
535 8 : END SUBROUTINE generate_derived_basis_sets
536 :
537 : ! **************************************************************************************************
538 : !> \brief Write a basis set file which can be used from CP2K
539 : !> \param basis ...
540 : !> \param element ...
541 : !> \param unit_nr ...
542 : !> \author Florian Schiffmann
543 : ! **************************************************************************************************
544 :
545 104 : SUBROUTINE write_basis(basis, element, unit_nr)
546 : TYPE(flex_basis_type) :: basis
547 : CHARACTER(LEN=default_string_length) :: element
548 : INTEGER :: unit_nr
549 :
550 : INTEGER :: iexp, iset
551 :
552 104 : IF (unit_nr > 0) THEN
553 52 : WRITE (UNIT=unit_nr, FMT="(A)") TRIM(ADJUSTL(element))//" "//TRIM(ADJUSTL(basis%basis_name))
554 52 : WRITE (UNIT=unit_nr, FMT="(1X,I0)") basis%nsets
555 104 : DO iset = 1, basis%nsets
556 52 : WRITE (UNIT=unit_nr, FMT="(30(1X,I0))") basis%subset(iset)%n, basis%subset(iset)%lmin, basis%subset(iset)%lmax, &
557 218 : basis%subset(iset)%nexp, basis%subset(iset)%l
558 466 : DO iexp = 1, basis%subset(iset)%nexp
559 : WRITE (UNIT=unit_nr, FMT="(T2,F24.14,30(1X,ES24.14))") &
560 1762 : basis%subset(iset)%exps(iexp), basis%subset(iset)%coeff(iexp, :)
561 : END DO
562 : END DO
563 : END IF
564 :
565 104 : END SUBROUTINE write_basis
566 :
567 : ! **************************************************************************************************
568 : !> \brief Initialize the derived basis sets and the vectors containing their
569 : !> assembly information ehich is used for regeneration of the sets.
570 : !> \param info_new ...
571 : !> \param info_ref ...
572 : !> \param basis ...
573 : !> \param basis_new ...
574 : !> \author Florian Schiffmann
575 : ! **************************************************************************************************
576 :
577 16 : SUBROUTINE setup_used_parts_init_basis(info_new, info_ref, basis, basis_new)
578 : TYPE(derived_basis_info) :: info_new, info_ref
579 : TYPE(flex_basis_type) :: basis, basis_new
580 :
581 : INTEGER :: i, jset, lind, nsets
582 :
583 : ! copy the reference information on the new set
584 :
585 48 : ALLOCATE (info_new%in_use_set(SIZE(info_ref%in_use_set)))
586 32 : info_new%in_use_set(:) = info_ref%in_use_set
587 64 : ALLOCATE (info_new%use_contr(SIZE(info_ref%in_use_set)))
588 32 : DO i = 1, SIZE(info_ref%in_use_set)
589 48 : ALLOCATE (info_new%use_contr(i)%in_use(SIZE(info_ref%use_contr(i)%in_use)))
590 120 : info_new%use_contr(i)%in_use(:) = info_ref%use_contr(i)%in_use
591 : END DO
592 16 : DO i = 1, info_new%nsets
593 16 : info_new%in_use_set(info_new%remove_set(i)) = .FALSE.
594 : END DO
595 48 : DO i = 1, info_new%ncontr
596 : lind = convert_l_contr_to_entry(basis%subset(info_new%remove_contr(i, 1))%lmin, &
597 : basis%subset(info_new%remove_contr(i, 1))%l, &
598 32 : info_new%remove_contr(i, 3), info_new%remove_contr(i, 2))
599 :
600 48 : info_new%use_contr(info_new%remove_contr(i, 1))%in_use(lind) = .FALSE.
601 : END DO
602 :
603 16 : nsets = 0
604 32 : DO i = 1, basis%nsets
605 32 : IF (info_new%in_use_set(i)) nsets = nsets + 1
606 : END DO
607 16 : basis_new%nsets = nsets
608 64 : ALLOCATE (basis_new%subset(nsets))
609 16 : jset = 0
610 32 : DO i = 1, basis%nsets
611 32 : IF (info_new%in_use_set(i)) THEN
612 16 : jset = jset + 1
613 16 : CALL create_new_subset(basis%subset(i), basis_new%subset(jset), info_new%use_contr(jset)%in_use)
614 : END IF
615 : END DO
616 :
617 16 : END SUBROUTINE setup_used_parts_init_basis
618 :
619 : ! **************************************************************************************************
620 : !> \brief Fill the low level information of the derived basis set.
621 : !> \param subset ...
622 : !> \param subset_new ...
623 : !> \param in_use ...
624 : !> \author Florian Schiffmann
625 : ! **************************************************************************************************
626 :
627 16 : SUBROUTINE create_new_subset(subset, subset_new, in_use)
628 : TYPE(subset_type) :: subset, subset_new
629 : LOGICAL, DIMENSION(:) :: in_use
630 :
631 : INTEGER :: icon, iind, il
632 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_l
633 :
634 48 : ALLOCATE (tmp_l(SIZE(subset%l)))
635 56 : tmp_l(:) = subset%l
636 16 : subset_new%lmin = subset%lmin
637 16 : subset_new%lmax = subset%lmin - 1
638 16 : subset_new%nexp = subset%nexp
639 16 : subset_new%n = subset%n
640 56 : DO il = 1, SIZE(subset%l)
641 128 : DO icon = 1, subset%l(il)
642 88 : iind = convert_l_contr_to_entry(subset%lmin, subset%l, icon, subset%lmin + il - 1)
643 128 : IF (.NOT. in_use(iind)) tmp_l(il) = tmp_l(il) - 1
644 : END DO
645 56 : IF (tmp_l(il) > 0) subset_new%lmax = subset_new%lmax + 1
646 : END DO
647 16 : subset_new%nl = subset_new%lmax - subset_new%lmin + 1
648 56 : subset_new%ncon_tot = SUM(tmp_l)
649 48 : ALLOCATE (subset_new%l(subset_new%nl))
650 64 : ALLOCATE (subset_new%coeff(subset_new%nexp, subset_new%ncon_tot))
651 48 : ALLOCATE (subset_new%exps(subset_new%nexp))
652 128 : subset_new%exps(:) = subset%exps
653 48 : DO il = 1, SIZE(subset%l)
654 40 : IF (tmp_l(il) == 0) EXIT
655 48 : subset_new%l(il) = tmp_l(il)
656 : END DO
657 16 : DEALLOCATE (tmp_l)
658 16 : iind = 0
659 104 : DO icon = 1, subset%ncon_tot
660 104 : IF (in_use(icon)) THEN
661 44 : iind = iind + 1
662 352 : subset_new%coeff(:, iind) = subset%coeff(:, icon)
663 : END IF
664 : END DO
665 :
666 16 : END SUBROUTINE create_new_subset
667 :
668 : ! **************************************************************************************************
669 : !> \brief for completeness generate the derived info for set 0(reference from file)
670 : !> \param info ...
671 : !> \param basis ...
672 : !> \author Florian Schiffmann
673 : ! **************************************************************************************************
674 :
675 16 : SUBROUTINE init_deriv_info_ref(info, basis)
676 : TYPE(derived_basis_info) :: info
677 : TYPE(flex_basis_type) :: basis
678 :
679 : INTEGER :: i
680 :
681 48 : ALLOCATE (info%in_use_set(basis%nsets))
682 32 : info%in_use_set = .TRUE.
683 64 : ALLOCATE (info%use_contr(basis%nsets))
684 32 : DO i = 1, basis%nsets
685 48 : ALLOCATE (info%use_contr(i)%in_use(basis%subset(i)%ncon_tot))
686 112 : info%use_contr(i)%in_use = .TRUE.
687 : END DO
688 :
689 16 : END SUBROUTINE init_deriv_info_ref
690 :
691 : ! **************************************************************************************************
692 : !> \brief get the general information for the basis sets. E.g. what to optimize,
693 : !> Basis set name, constraints upon optimization and read the reference basis
694 : !> \param kind_section ...
695 : !> \param opt_bas ...
696 : !> \param para_env ...
697 : !> \author Florian Schiffmann
698 : ! **************************************************************************************************
699 :
700 8 : SUBROUTINE generate_initial_basis(kind_section, opt_bas, para_env)
701 : TYPE(section_vals_type), POINTER :: kind_section
702 : TYPE(basis_optimization_type) :: opt_bas
703 : TYPE(mp_para_env_type), POINTER :: para_env
704 :
705 : INTEGER :: ikind, variable_counter
706 : LOGICAL :: explicit
707 : TYPE(section_vals_type), POINTER :: set_section
708 :
709 8 : CALL section_vals_get(kind_section, n_repetition=opt_bas%nkind)
710 40 : ALLOCATE (opt_bas%kind_basis(opt_bas%nkind))
711 :
712 : ! counter to get the number of free variables in optimization
713 8 : variable_counter = 0
714 24 : DO ikind = 1, opt_bas%nkind
715 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", c_val=opt_bas%kind_basis(ikind)%element, &
716 16 : i_rep_section=ikind)
717 : CALL section_vals_val_get(kind_section, "BASIS_SET", c_val=opt_bas%kind_basis(ikind)%basis_name, &
718 16 : i_rep_section=ikind)
719 : set_section => section_vals_get_subs_vals(kind_section, "DERIVED_BASIS_SETS", &
720 16 : i_rep_section=ikind)
721 16 : CALL section_vals_get(set_section, n_repetition=opt_bas%kind_basis(ikind)%nbasis_deriv, explicit=explicit)
722 16 : IF (.NOT. explicit) opt_bas%kind_basis(ikind)%nbasis_deriv = 0
723 80 : ALLOCATE (opt_bas%kind_basis(ikind)%flex_basis(0:opt_bas%kind_basis(ikind)%nbasis_deriv))
724 80 : ALLOCATE (opt_bas%kind_basis(ikind)%deriv_info(0:opt_bas%kind_basis(ikind)%nbasis_deriv))
725 :
726 : CALL fill_basis_template(kind_section, opt_bas%kind_basis(ikind)%flex_basis(0), opt_bas%template_basis_file, &
727 16 : opt_bas%kind_basis(ikind)%element, opt_bas%kind_basis(ikind)%basis_name, para_env, ikind)
728 :
729 16 : CALL setup_exp_constraints(kind_section, opt_bas%kind_basis(ikind)%flex_basis(0))
730 :
731 16 : CALL parse_derived_basis(kind_section, opt_bas%kind_basis(ikind)%deriv_info, ikind)
732 :
733 40 : variable_counter = variable_counter + opt_bas%kind_basis(ikind)%flex_basis(0)%nopt
734 : END DO
735 :
736 24 : ALLOCATE (opt_bas%x_opt(variable_counter))
737 :
738 8 : variable_counter = 0
739 24 : DO ikind = 1, opt_bas%nkind
740 24 : CALL assign_x_to_basis(opt_bas%x_opt, opt_bas%kind_basis(ikind)%flex_basis(0), variable_counter)
741 : END DO
742 :
743 8 : CPASSERT(variable_counter == SIZE(opt_bas%x_opt))
744 :
745 8 : END SUBROUTINE generate_initial_basis
746 :
747 : ! **************************************************************************************************
748 : !> \brief get low level information about how to construc new basis sets from reference
749 : !> \param kind_section ...
750 : !> \param deriv_info ...
751 : !> \param ikind ...
752 : !> \author Florian Schiffmann
753 : ! **************************************************************************************************
754 :
755 16 : SUBROUTINE parse_derived_basis(kind_section, deriv_info, ikind)
756 : TYPE(section_vals_type), POINTER :: kind_section
757 : TYPE(derived_basis_info), DIMENSION(:) :: deriv_info
758 : INTEGER :: ikind
759 :
760 : INTEGER :: i_rep, iset, jset, n_rep, nsets
761 16 : INTEGER, DIMENSION(:), POINTER :: i_vals
762 : LOGICAL :: explicit
763 : TYPE(section_vals_type), POINTER :: set1_section
764 :
765 16 : nsets = SIZE(deriv_info) - 1
766 : set1_section => section_vals_get_subs_vals(kind_section, "DERIVED_BASIS_SETS", &
767 32 : i_rep_section=ikind)
768 32 : DO jset = 1, nsets
769 : ! stracnge but as derive info is allcated from 0 to n the count over here has to be shifted
770 16 : iset = jset + 1
771 : CALL section_vals_val_get(set1_section, "BASIS_SET_NAME", c_val=deriv_info(iset)%basis_name, &
772 16 : i_rep_section=jset)
773 16 : CALL section_vals_val_get(set1_section, "REFERENCE_SET", i_vals=i_vals, i_rep_section=jset)
774 16 : deriv_info(iset)%reference_set = i_vals(1)
775 : CALL section_vals_val_get(set1_section, "REMOVE_CONTRACTION", explicit=explicit, n_rep_val=n_rep, &
776 16 : i_rep_section=jset)
777 16 : deriv_info(iset)%ncontr = n_rep
778 16 : IF (explicit) THEN
779 48 : ALLOCATE (deriv_info(iset)%remove_contr(n_rep, 3))
780 48 : DO i_rep = 1, n_rep
781 : CALL section_vals_val_get(set1_section, "REMOVE_CONTRACTION", i_rep_val=i_rep, i_vals=i_vals, &
782 32 : i_rep_section=jset)
783 144 : deriv_info(iset)%remove_contr(i_rep, :) = i_vals(:)
784 : END DO
785 : END IF
786 : CALL section_vals_val_get(set1_section, "REMOVE_SET", explicit=explicit, n_rep_val=n_rep, &
787 16 : i_rep_section=jset)
788 16 : deriv_info(iset)%nsets = n_rep
789 48 : IF (explicit) THEN
790 0 : ALLOCATE (deriv_info(iset)%remove_set(n_rep))
791 0 : DO i_rep = 1, n_rep
792 : CALL section_vals_val_get(set1_section, "REMOVE_SET", i_rep_val=i_rep, i_vals=i_vals, &
793 0 : i_rep_section=jset)
794 0 : deriv_info(iset)%remove_set(i_rep) = i_vals(1)
795 : END DO
796 : END IF
797 : END DO
798 :
799 16 : END SUBROUTINE parse_derived_basis
800 :
801 : ! **************************************************************************************************
802 : !> \brief get low level information about constraint on exponents
803 : !> \param kind1_section ...
804 : !> \param flex_basis ...
805 : !> \author Florian Schiffmann
806 : ! **************************************************************************************************
807 :
808 32 : SUBROUTINE setup_exp_constraints(kind1_section, flex_basis)
809 : TYPE(section_vals_type), POINTER :: kind1_section
810 : TYPE(flex_basis_type) :: flex_basis
811 :
812 : INTEGER :: ipgf, irep, iset, nrep
813 16 : INTEGER, DIMENSION(:), POINTER :: def_exp
814 : LOGICAL :: is_bound, is_varlim
815 : TYPE(section_vals_type), POINTER :: const_section
816 :
817 32 : const_section => section_vals_get_subs_vals(kind1_section, "CONSTRAIN_EXPONENTS")
818 16 : CALL section_vals_get(const_section, n_repetition=nrep)
819 16 : DO irep = 1, nrep
820 0 : CALL section_vals_val_get(const_section, "USE_EXP", i_vals=def_exp, i_rep_section=irep)
821 0 : CALL section_vals_val_get(const_section, "BOUNDARIES", explicit=is_bound, i_rep_section=irep)
822 0 : CALL section_vals_val_get(const_section, "MAX_VAR_FRACTION", explicit=is_varlim, i_rep_section=irep)
823 0 : IF (is_bound .AND. is_varlim) THEN
824 : CALL cp_abort(__LOCATION__, "Exponent has two constraints. "// &
825 0 : "This is not possible at the moment. Please change input.")
826 : END IF
827 0 : IF (.NOT. is_bound .AND. .NOT. is_varlim) THEN
828 : CALL cp_abort(__LOCATION__, "Exponent is declared to be constraint but none is given"// &
829 0 : " Please change input.")
830 : END IF
831 16 : IF (def_exp(1) == -1) THEN
832 0 : DO iset = 1, flex_basis%nsets
833 0 : IF (def_exp(2) == -1) THEN
834 0 : DO ipgf = 1, flex_basis%subset(iset)%nexp
835 0 : CALL set_constraint(flex_basis, iset, ipgf, const_section, is_bound, is_varlim, irep)
836 : END DO
837 : ELSE
838 0 : IF (def_exp(2) <= flex_basis%subset(iset)%nexp) THEN
839 : CALL cp_abort(__LOCATION__, &
840 : "Exponent declared in constraint is larger than number of exponents in the set"// &
841 0 : " Please change input.")
842 : END IF
843 0 : CALL set_constraint(flex_basis, iset, def_exp(2), const_section, is_bound, is_varlim, irep)
844 : END IF
845 : END DO
846 : ELSE
847 0 : IF (.NOT. def_exp(1) <= flex_basis%nsets) THEN
848 : CALL cp_abort(__LOCATION__, &
849 : "Set number of constraint is larger than number of sets in the template basis set."// &
850 0 : " Please change input.")
851 : END IF
852 0 : IF (def_exp(2) == -1) THEN
853 0 : DO ipgf = 1, flex_basis%subset(iset)%nexp
854 0 : CALL set_constraint(flex_basis, def_exp(1), ipgf, const_section, is_bound, is_varlim, irep)
855 : END DO
856 : ELSE
857 0 : IF (.NOT. def_exp(2) <= flex_basis%subset(def_exp(1))%nexp) THEN
858 : CALL cp_abort(__LOCATION__, &
859 : "Exponent declared in constraint is larger than number of exponents in the set"// &
860 0 : " Please change input.")
861 : END IF
862 0 : CALL set_constraint(flex_basis, def_exp(1), def_exp(2), const_section, is_bound, is_varlim, irep)
863 : END IF
864 : END IF
865 : END DO
866 :
867 16 : END SUBROUTINE setup_exp_constraints
868 :
869 : ! **************************************************************************************************
870 : !> \brief put the constraint information in type and process if requires
871 : !> BOUNDARIES constraint gets transformed into MAX_VAR_FRACTION constraint.
872 : !> \param flex_basis ...
873 : !> \param iset ...
874 : !> \param ipgf ...
875 : !> \param const_section ...
876 : !> \param is_bound ...
877 : !> \param is_varlim ...
878 : !> \param irep ...
879 : !> \author Florian Schiffmann
880 : ! **************************************************************************************************
881 :
882 0 : SUBROUTINE set_constraint(flex_basis, iset, ipgf, const_section, is_bound, is_varlim, irep)
883 : TYPE(flex_basis_type) :: flex_basis
884 : INTEGER :: iset, ipgf
885 : TYPE(section_vals_type), POINTER :: const_section
886 : LOGICAL :: is_bound, is_varlim
887 : INTEGER :: irep
888 :
889 : REAL(KIND=dp) :: r_val
890 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: r_vals
891 :
892 0 : IF (flex_basis%subset(iset)%exp_has_const(ipgf)) THEN
893 : CALL cp_abort(__LOCATION__, &
894 : "Multiple constraints due to collision in CONSTRAIN_EXPONENTS."// &
895 0 : " Please change input.")
896 : END IF
897 0 : flex_basis%subset(iset)%exp_has_const(ipgf) = .TRUE.
898 0 : IF (is_bound) THEN
899 0 : flex_basis%subset(iset)%exp_const(ipgf)%const_type = 0
900 0 : CALL section_vals_val_get(const_section, "BOUNDARIES", r_vals=r_vals, i_rep_section=irep)
901 0 : flex_basis%subset(iset)%exp_const(ipgf)%llim = MINVAL(r_vals)
902 0 : flex_basis%subset(iset)%exp_const(ipgf)%ulim = MAXVAL(r_vals)
903 0 : r_val = flex_basis%subset(iset)%exps(ipgf)
904 0 : IF (flex_basis%subset(iset)%exps(ipgf) > MAXVAL(r_vals) .OR. flex_basis%subset(iset)%exps(ipgf) < MINVAL(r_vals)) THEN
905 : CALL cp_abort(__LOCATION__, &
906 : "Exponent "//cp_to_string(r_val)// &
907 : " declared in constraint is out of bounds of constraint"//cp_to_string(MINVAL(r_vals))// &
908 : " to"//cp_to_string(MAXVAL(r_vals))// &
909 0 : " Please change input.")
910 : END IF
911 0 : flex_basis%subset(iset)%exp_const(ipgf)%init = SUM(r_vals)/2.0_dp
912 0 : flex_basis%subset(iset)%exp_const(ipgf)%var_fac = MAXVAL(r_vals)/flex_basis%subset(iset)%exp_const(ipgf)%init - 1.0_dp
913 : END IF
914 0 : IF (is_varlim) THEN
915 0 : flex_basis%subset(iset)%exp_const(ipgf)%const_type = 1
916 0 : CALL section_vals_val_get(const_section, "MAX_VAR_FRACTION", r_vals=r_vals, i_rep_section=irep)
917 0 : flex_basis%subset(iset)%exp_const(ipgf)%var_fac = r_vals(1)
918 0 : flex_basis%subset(iset)%exp_const(ipgf)%init = flex_basis%subset(iset)%exps(ipgf)
919 : END IF
920 :
921 0 : END SUBROUTINE set_constraint
922 :
923 : ! **************************************************************************************************
924 : !> \brief Initialize the optimization vector with the values from the refernece sets
925 : !> \param x ...
926 : !> \param basis ...
927 : !> \param x_ind ...
928 : !> \author Florian Schiffmann
929 : ! **************************************************************************************************
930 :
931 16 : SUBROUTINE assign_x_to_basis(x, basis, x_ind)
932 : REAL(KIND=dp), DIMENSION(:) :: x
933 : TYPE(flex_basis_type) :: basis
934 : INTEGER :: x_ind
935 :
936 : INTEGER :: icont, ipgf, iset
937 :
938 32 : DO iset = 1, basis%nsets
939 140 : DO ipgf = 1, basis%subset(iset)%nexp
940 108 : IF (basis%subset(iset)%opt_exps(ipgf)) THEN
941 0 : x_ind = x_ind + 1
942 0 : basis%subset(iset)%exp_x_ind(ipgf) = x_ind
943 0 : x(x_ind) = basis%subset(iset)%exps(ipgf)
944 : END IF
945 664 : DO icont = 1, basis%subset(iset)%ncon_tot
946 648 : IF (basis%subset(iset)%opt_coeff(ipgf, icont)) THEN
947 380 : x_ind = x_ind + 1
948 380 : basis%subset(iset)%coeff_x_ind(ipgf, icont) = x_ind
949 380 : x(x_ind) = basis%subset(iset)%coeff(ipgf, icont)
950 : END IF
951 : END DO
952 : END DO
953 : END DO
954 :
955 16 : END SUBROUTINE assign_x_to_basis
956 :
957 : ! **************************************************************************************************
958 : !> \brief Fill the reference set and get the free varialbles from input
959 : !> \param kind1_section ...
960 : !> \param flex_basis ...
961 : !> \param template_basis_file ...
962 : !> \param element ...
963 : !> \param basis_name ...
964 : !> \param para_env ...
965 : !> \param ikind ...
966 : !> \author Florian Schiffmann
967 : ! **************************************************************************************************
968 :
969 96 : SUBROUTINE fill_basis_template(kind1_section, flex_basis, template_basis_file, element, basis_name, para_env, ikind)
970 : TYPE(section_vals_type), POINTER :: kind1_section
971 : TYPE(flex_basis_type) :: flex_basis
972 : CHARACTER(LEN=default_path_length) :: template_basis_file
973 : CHARACTER(LEN=default_string_length) :: element, basis_name
974 : TYPE(mp_para_env_type), POINTER :: para_env
975 : INTEGER :: ikind
976 :
977 : INTEGER :: icont, idof, ipgf, irep, iset, nrep
978 16 : INTEGER, DIMENSION(:), POINTER :: switch
979 :
980 16 : CALL parse_basis(flex_basis, template_basis_file, element, basis_name, para_env)
981 :
982 : ! get the optimizable parameters. Many way to modify them but in the end only logical matrix
983 : ! is either set or values get flipped according to the input
984 : CALL section_vals_val_get(kind1_section, "INITIAL_DEGREES_OF_FREEDOM", i_val=idof, &
985 16 : i_rep_section=ikind)
986 32 : DO iset = 1, flex_basis%nsets
987 16 : SELECT CASE (idof)
988 : CASE (do_opt_none)
989 : ! initialization in parse subset did the job
990 : CASE (do_opt_all)
991 0 : flex_basis%subset(iset)%opt_coeff = .TRUE.
992 0 : flex_basis%subset(iset)%opt_exps = .TRUE.
993 : CASE (do_opt_coeff)
994 360 : flex_basis%subset(iset)%opt_coeff = .TRUE.
995 : CASE (do_opt_exps)
996 0 : flex_basis%subset(iset)%opt_exps = .TRUE.
997 : CASE DEFAULT
998 16 : CPABORT("No initialization available?????")
999 : END SELECT
1000 : END DO
1001 :
1002 16 : CALL section_vals_val_get(kind1_section, "SWITCH_CONTRACTION_STATE", n_rep_val=nrep, i_rep_section=ikind)
1003 28 : DO irep = 1, nrep
1004 : CALL section_vals_val_get(kind1_section, "SWITCH_CONTRACTION_STATE", i_rep_val=irep, &
1005 12 : i_rep_section=ikind, i_vals=switch)
1006 12 : icont = convert_l_contr_to_entry(flex_basis%subset(switch(1))%lmin, flex_basis%subset(switch(1))%l, switch(3), switch(2))
1007 100 : DO ipgf = 1, flex_basis%subset(switch(1))%nexp
1008 84 : flex_basis%subset(switch(1))%opt_coeff(ipgf, icont) = .NOT. flex_basis%subset(switch(1))%opt_coeff(ipgf, icont)
1009 : END DO
1010 : END DO
1011 :
1012 16 : CALL section_vals_val_get(kind1_section, "SWITCH_COEFF_STATE", n_rep_val=nrep, i_rep_section=ikind)
1013 16 : DO irep = 1, nrep
1014 : CALL section_vals_val_get(kind1_section, "SWITCH_COEFF_STATE", i_rep_val=irep, &
1015 0 : i_rep_section=ikind, i_vals=switch)
1016 0 : icont = convert_l_contr_to_entry(flex_basis%subset(switch(1))%lmin, flex_basis%subset(switch(1))%l, switch(3), switch(2))
1017 : flex_basis%subset(switch(1))%opt_coeff(switch(4), icont) = &
1018 16 : .NOT. flex_basis%subset(switch(1))%opt_coeff(switch(4), icont)
1019 : END DO
1020 :
1021 16 : CALL section_vals_val_get(kind1_section, "SWITCH_EXP_STATE", n_rep_val=nrep, i_rep_section=ikind)
1022 16 : DO irep = 1, nrep
1023 : CALL section_vals_val_get(kind1_section, "SWITCH_EXP_STATE", i_rep_val=irep, &
1024 0 : i_rep_section=ikind, i_vals=switch)
1025 16 : flex_basis%subset(switch(1))%opt_exps(switch(2)) = .NOT. flex_basis%subset(switch(1))%opt_exps(switch(2))
1026 : END DO
1027 :
1028 16 : CALL section_vals_val_get(kind1_section, "SWITCH_SET_STATE", n_rep_val=nrep, i_rep_section=ikind)
1029 16 : DO irep = 1, nrep
1030 : CALL section_vals_val_get(kind1_section, "SWITCH_SET_STATE", i_rep_val=irep, &
1031 0 : i_rep_section=ikind, i_vals=switch)
1032 16 : DO ipgf = 1, flex_basis%subset(switch(2))%nexp
1033 0 : SELECT CASE (switch(1))
1034 : CASE (0) ! switch all states in the set
1035 0 : DO icont = 1, flex_basis%subset(switch(2))%ncon_tot
1036 : flex_basis%subset(switch(2))%opt_coeff(ipgf, icont) = &
1037 0 : .NOT. flex_basis%subset(switch(2))%opt_coeff(ipgf, icont)
1038 : END DO
1039 0 : flex_basis%subset(switch(2))%opt_exps(ipgf) = .NOT. flex_basis%subset(switch(2))%opt_exps(ipgf)
1040 : CASE (1) ! switch only exp
1041 0 : flex_basis%subset(switch(2))%opt_exps(ipgf) = .NOT. flex_basis%subset(switch(2))%opt_exps(ipgf)
1042 : CASE (2) ! switch only coeff
1043 0 : DO icont = 1, flex_basis%subset(switch(2))%ncon_tot
1044 : flex_basis%subset(switch(2))%opt_coeff(ipgf, icont) = &
1045 0 : .NOT. flex_basis%subset(switch(2))%opt_coeff(ipgf, icont)
1046 : END DO
1047 : CASE DEFAULT
1048 0 : CPABORT("Invalid option in SWITCH_SET_STATE, 1st value has to be 0, 1 or 2")
1049 : END SELECT
1050 : END DO
1051 : END DO
1052 :
1053 : ! perform a final modification. If basis set is uncontracted coefficient will never have to be optimized
1054 32 : DO irep = 1, flex_basis%nsets
1055 32 : IF (flex_basis%subset(irep)%nexp == 1) THEN
1056 0 : DO ipgf = 1, flex_basis%subset(irep)%nexp
1057 0 : flex_basis%subset(irep)%opt_coeff(ipgf, 1) = .FALSE.
1058 : END DO
1059 : END IF
1060 : END DO
1061 :
1062 : ! finally count the total number of free parameters
1063 16 : flex_basis%nopt = 0
1064 32 : DO irep = 1, flex_basis%nsets
1065 140 : DO ipgf = 1, flex_basis%subset(irep)%nexp
1066 648 : DO icont = 1, flex_basis%subset(irep)%ncon_tot
1067 648 : IF (flex_basis%subset(irep)%opt_coeff(ipgf, icont)) flex_basis%nopt = flex_basis%nopt + 1
1068 : END DO
1069 124 : IF (flex_basis%subset(irep)%opt_exps(ipgf)) flex_basis%nopt = flex_basis%nopt + 1
1070 : END DO
1071 : END DO
1072 :
1073 16 : END SUBROUTINE fill_basis_template
1074 :
1075 : ! **************************************************************************************************
1076 : !> \brief Helper function to parse input. Converts l and index position of
1077 : !> a contraction to index in the contraction array of the set using lmin and nl
1078 : !> \param lmin ...
1079 : !> \param nl ...
1080 : !> \param icontr ...
1081 : !> \param l ...
1082 : !> \return ...
1083 : !> \author Florian Schiffmann
1084 : ! **************************************************************************************************
1085 :
1086 132 : FUNCTION convert_l_contr_to_entry(lmin, nl, icontr, l) RESULT(ientry)
1087 : INTEGER :: lmin
1088 : INTEGER, DIMENSION(:) :: nl
1089 : INTEGER :: icontr, l, ientry
1090 :
1091 : INTEGER :: i, icon2l, iwork
1092 :
1093 132 : iwork = l - lmin
1094 132 : icon2l = 0
1095 212 : DO i = 1, iwork
1096 212 : icon2l = icon2l + nl(i)
1097 : END DO
1098 132 : ientry = icon2l + icontr
1099 :
1100 132 : END FUNCTION convert_l_contr_to_entry
1101 :
1102 : ! **************************************************************************************************
1103 : !> \brief Read the reference basis sets from the template basis file
1104 : !> \param flex_basis ...
1105 : !> \param template_basis_file ...
1106 : !> \param element ...
1107 : !> \param basis_name ...
1108 : !> \param para_env ...
1109 : !> \author Florian Schiffmann
1110 : ! **************************************************************************************************
1111 :
1112 32 : SUBROUTINE parse_basis(flex_basis, template_basis_file, element, basis_name, para_env)
1113 : TYPE(flex_basis_type) :: flex_basis
1114 : CHARACTER(LEN=default_path_length) :: template_basis_file
1115 : CHARACTER(LEN=default_string_length) :: element, basis_name
1116 : TYPE(mp_para_env_type), POINTER :: para_env
1117 :
1118 : CHARACTER(LEN=240) :: line
1119 : CHARACTER(LEN=242) :: line2
1120 : CHARACTER(LEN=LEN(basis_name)+2) :: basis_name2
1121 : CHARACTER(LEN=LEN(element)+2) :: element2
1122 : INTEGER :: iset, strlen1, strlen2
1123 : LOGICAL :: basis_found, found, match
1124 : TYPE(cp_parser_type) :: parser
1125 :
1126 16 : basis_found = .FALSE.
1127 16 : CALL uppercase(element)
1128 16 : CALL uppercase(basis_name)
1129 16 : CALL parser_create(parser, template_basis_file, para_env=para_env)
1130 :
1131 : search_loop: DO
1132 32 : CALL parser_search_string(parser, TRIM(basis_name), .TRUE., found, line)
1133 32 : IF (found) THEN
1134 32 : match = .FALSE.
1135 32 : CALL uppercase(line)
1136 : ! Check both the element symbol and the basis set name
1137 32 : line2 = " "//line//" "
1138 32 : element2 = " "//TRIM(element)//" "
1139 32 : basis_name2 = " "//TRIM(basis_name)//" "
1140 32 : strlen1 = LEN_TRIM(element2) + 1
1141 32 : strlen2 = LEN_TRIM(basis_name2) + 1
1142 32 : IF ((INDEX(line2, element2(:strlen1)) > 0) .AND. &
1143 16 : (INDEX(line2, basis_name2(:strlen2)) > 0)) match = .TRUE.
1144 32 : IF (match) THEN
1145 16 : CALL parser_get_object(parser, flex_basis%nsets, newline=.TRUE.)
1146 64 : ALLOCATE (flex_basis%subset(flex_basis%nsets))
1147 32 : DO iset = 1, flex_basis%nsets
1148 32 : CALL parse_subset(parser, flex_basis%subset(iset))
1149 : END DO
1150 : basis_found = .TRUE.
1151 : EXIT search_loop
1152 : END IF
1153 : ELSE
1154 : EXIT search_loop
1155 : END IF
1156 : END DO search_loop
1157 16 : CALL parser_release(parser)
1158 :
1159 16 : IF (.NOT. basis_found) CALL cp_abort(__LOCATION__, &
1160 : "The requested basis set <"//TRIM(basis_name)// &
1161 : "> for element <"//TRIM(element)//"> was not "// &
1162 0 : "found in the template basis set file ")
1163 :
1164 48 : END SUBROUTINE parse_basis
1165 :
1166 : ! **************************************************************************************************
1167 : !> \brief Read the subset information from the template basis file
1168 : !> \param parser ...
1169 : !> \param subset ...
1170 : !> \author Florian Schiffmann
1171 : ! **************************************************************************************************
1172 16 : SUBROUTINE parse_subset(parser, subset)
1173 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
1174 : TYPE(subset_type) :: subset
1175 :
1176 : CHARACTER(len=20*default_string_length) :: line_att
1177 : INTEGER :: icon1, icon2, il, ipgf, ishell, istart
1178 : REAL(KIND=dp) :: gs_scale
1179 : REAL(KIND=dp), POINTER :: r_val
1180 :
1181 : line_att = ""
1182 16 : CALL parser_get_object(parser, subset%n, newline=.TRUE.)
1183 16 : CALL parser_get_object(parser, subset%lmin)
1184 16 : CALL parser_get_object(parser, subset%lmax)
1185 16 : CALL parser_get_object(parser, subset%nexp)
1186 16 : subset%nl = subset%lmax - subset%lmin + 1
1187 16 : ALLOCATE (r_val)
1188 48 : ALLOCATE (subset%l(subset%nl))
1189 48 : ALLOCATE (subset%exps(subset%nexp))
1190 48 : ALLOCATE (subset%exp_has_const(subset%nexp))
1191 124 : subset%exp_has_const = .FALSE.
1192 32 : ALLOCATE (subset%opt_exps(subset%nexp))
1193 124 : subset%opt_exps = .FALSE.
1194 156 : ALLOCATE (subset%exp_const(subset%nexp))
1195 32 : ALLOCATE (subset%exp_x_ind(subset%nexp))
1196 56 : DO ishell = 1, subset%nl
1197 56 : CALL parser_get_object(parser, subset%l(ishell))
1198 : END DO
1199 56 : subset%ncon_tot = SUM(subset%l)
1200 64 : ALLOCATE (subset%coeff(subset%nexp, subset%ncon_tot))
1201 64 : ALLOCATE (subset%opt_coeff(subset%nexp, subset%ncon_tot))
1202 636 : subset%opt_coeff = .FALSE.
1203 48 : ALLOCATE (subset%coeff_x_ind(subset%nexp, subset%ncon_tot))
1204 124 : DO ipgf = 1, subset%nexp
1205 108 : CALL parser_get_object(parser, r_val, newline=.TRUE.)
1206 108 : subset%exps(ipgf) = r_val
1207 664 : DO ishell = 1, subset%ncon_tot
1208 540 : CALL parser_get_object(parser, r_val)
1209 648 : subset%coeff(ipgf, ishell) = r_val
1210 : END DO
1211 : END DO
1212 :
1213 : ! orthonormalize contraction coefficients using gram schmidt
1214 16 : istart = 1
1215 56 : DO il = 1, subset%nl
1216 80 : DO icon1 = istart, istart + subset%l(il) - 2
1217 136 : DO icon2 = icon1 + 1, istart + subset%l(il) - 1
1218 : gs_scale = DOT_PRODUCT(subset%coeff(:, icon2), subset%coeff(:, icon1))/ &
1219 824 : DOT_PRODUCT(subset%coeff(:, icon1), subset%coeff(:, icon1))
1220 480 : subset%coeff(:, icon2) = subset%coeff(:, icon2) - gs_scale*subset%coeff(:, icon1)
1221 : END DO
1222 : END DO
1223 56 : istart = istart + subset%l(il)
1224 : END DO
1225 :
1226 : ! just to get an understandable basis normalize coefficients
1227 96 : DO icon1 = 1, subset%ncon_tot
1228 1176 : subset%coeff(:, icon1) = subset%coeff(:, icon1)/NORM2(subset%coeff(:, icon1))
1229 : END DO
1230 16 : DEALLOCATE (r_val)
1231 :
1232 16 : END SUBROUTINE parse_subset
1233 :
1234 : ! **************************************************************************************************
1235 : !> \brief Initialize the variables for the powell optimizer
1236 : !> \param p_param ...
1237 : !> \param powell_section ...
1238 : !> \author Florian Schiffmann
1239 : ! **************************************************************************************************
1240 :
1241 8 : SUBROUTINE init_powell_var(p_param, powell_section)
1242 : TYPE(opt_state_type), INTENT(INOUT) :: p_param
1243 : TYPE(section_vals_type), POINTER :: powell_section
1244 :
1245 8 : p_param%state = 0
1246 8 : p_param%nvar = 0
1247 8 : p_param%iprint = 0
1248 8 : p_param%unit = default_output_unit
1249 8 : CALL section_vals_val_get(powell_section, "ACCURACY", r_val=p_param%rhoend)
1250 8 : CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=p_param%rhobeg)
1251 8 : CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=p_param%maxfun)
1252 :
1253 8 : END SUBROUTINE init_powell_var
1254 :
1255 : END MODULE optimize_basis_utils
|
