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