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