Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Optimizes exponents and contraction coefficients of the lri auxiliary
10 : !> basis sets using the UOBYQA minimizer
11 : !> lri : local resolution of the identity
12 : !> \par History
13 : !> created Dorothea Golze [05.2014]
14 : !> \authors Dorothea Golze
15 : ! **************************************************************************************************
16 : MODULE lri_optimize_ri_basis
17 :
18 : USE atomic_kind_types, ONLY: atomic_kind_type
19 : USE basis_set_types, ONLY: get_gto_basis_set,&
20 : gto_basis_set_type,&
21 : init_orb_basis_set
22 : USE cell_types, ONLY: cell_type
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_type
25 : USE cp_output_handling, ONLY: cp_p_file,&
26 : cp_print_key_finished_output,&
27 : cp_print_key_generate_filename,&
28 : cp_print_key_should_output,&
29 : cp_print_key_unit_nr
30 : USE dbcsr_api, ONLY: dbcsr_get_block_p,&
31 : dbcsr_p_type,&
32 : dbcsr_type
33 : USE generic_os_integrals, ONLY: int_overlap_aabb_os
34 : USE input_constants, ONLY: do_lri_opt_all,&
35 : do_lri_opt_coeff,&
36 : do_lri_opt_exps
37 : USE input_section_types, ONLY: section_vals_get,&
38 : section_vals_get_subs_vals,&
39 : section_vals_type,&
40 : section_vals_val_get
41 : USE kinds, ONLY: default_path_length,&
42 : dp
43 : USE lri_environment_init, ONLY: lri_basis_init
44 : USE lri_environment_methods, ONLY: calculate_avec_lri,&
45 : calculate_lri_integrals
46 : USE lri_environment_types, ONLY: allocate_lri_ints_rho,&
47 : deallocate_lri_ints_rho,&
48 : lri_density_type,&
49 : lri_environment_type,&
50 : lri_int_rho_type,&
51 : lri_int_type,&
52 : lri_list_type,&
53 : lri_rhoab_type
54 : USE lri_optimize_ri_basis_types, ONLY: create_lri_opt,&
55 : deallocate_lri_opt,&
56 : get_original_gcc,&
57 : lri_opt_type,&
58 : orthonormalize_gcc
59 : USE memory_utilities, ONLY: reallocate
60 : USE message_passing, ONLY: mp_para_env_type
61 : USE particle_types, ONLY: particle_type
62 : USE powell, ONLY: opt_state_type,&
63 : powell_optimize
64 : USE qs_environment_types, ONLY: get_qs_env,&
65 : qs_environment_type,&
66 : set_qs_env
67 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
68 : neighbor_list_iterate,&
69 : neighbor_list_iterator_create,&
70 : neighbor_list_iterator_p_type,&
71 : neighbor_list_iterator_release,&
72 : neighbor_list_set_p_type
73 : USE qs_rho_types, ONLY: qs_rho_get,&
74 : qs_rho_type
75 :
76 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
77 : #include "./base/base_uses.f90"
78 :
79 : IMPLICIT NONE
80 :
81 : PRIVATE
82 :
83 : ! **************************************************************************************************
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_optimize_ri_basis'
86 :
87 : PUBLIC :: optimize_lri_basis, &
88 : get_condition_number_of_overlap
89 :
90 : ! **************************************************************************************************
91 :
92 : CONTAINS
93 :
94 : ! **************************************************************************************************
95 : !> \brief optimizes the lri basis set
96 : !> \param qs_env qs environment
97 : ! **************************************************************************************************
98 6 : SUBROUTINE optimize_lri_basis(qs_env)
99 :
100 : TYPE(qs_environment_type), POINTER :: qs_env
101 :
102 : INTEGER :: iunit, nkind
103 6 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
104 : TYPE(cp_logger_type), POINTER :: logger
105 6 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
106 : TYPE(lri_density_type), POINTER :: lri_density
107 : TYPE(lri_environment_type), POINTER :: lri_env
108 : TYPE(lri_opt_type), POINTER :: lri_opt
109 : TYPE(mp_para_env_type), POINTER :: para_env
110 : TYPE(opt_state_type) :: opt_state
111 : TYPE(qs_rho_type), POINTER :: rho_struct
112 : TYPE(section_vals_type), POINTER :: dft_section, input, lri_optbas_section
113 :
114 6 : NULLIFY (atomic_kind_set, dft_section, lri_density, lri_env, &
115 6 : lri_opt, lri_optbas_section, rho_struct)
116 6 : NULLIFY (input, logger, para_env)
117 :
118 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, input=input, &
119 : lri_env=lri_env, lri_density=lri_density, nkind=nkind, &
120 6 : para_env=para_env, rho=rho_struct)
121 :
122 : ! density matrix
123 6 : CALL qs_rho_get(rho_struct, rho_ao_kp=pmatrix)
124 :
125 6 : logger => cp_get_default_logger()
126 6 : dft_section => section_vals_get_subs_vals(input, "DFT")
127 : lri_optbas_section => section_vals_get_subs_vals(input, &
128 6 : "DFT%QS%OPTIMIZE_LRI_BASIS")
129 : iunit = cp_print_key_unit_nr(logger, input, "PRINT%PROGRAM_RUN_INFO", &
130 6 : extension=".opt")
131 :
132 6 : IF (iunit > 0) THEN
133 3 : WRITE (iunit, '(/," POWELL| Start optimization procedure")')
134 : END IF
135 :
136 : ! *** initialization
137 6 : CALL create_lri_opt(lri_opt)
138 : CALL init_optimization(lri_env, lri_opt, lri_optbas_section, &
139 6 : opt_state, lri_opt%x, lri_opt%zet_init, nkind, iunit)
140 :
141 6 : CALL calculate_lri_overlap_aabb(lri_env, qs_env)
142 :
143 : ! *** ======================= START optimization =====================
144 6 : opt_state%state = 0
145 30 : DO
146 36 : IF (opt_state%state == 2) THEN
147 : CALL calc_lri_integrals_get_objective(lri_env, lri_density, qs_env, &
148 : lri_opt, opt_state, pmatrix, para_env, &
149 18 : nkind)
150 : ! lri_density has been re-initialized!
151 18 : CALL set_qs_env(qs_env, lri_density=lri_density)
152 : END IF
153 :
154 36 : IF (opt_state%state == -1) EXIT
155 :
156 30 : CALL powell_optimize(opt_state%nvar, lri_opt%x, opt_state)
157 30 : CALL update_exponents(lri_env, lri_opt, lri_opt%x, lri_opt%zet_init, nkind)
158 30 : CALL print_optimization_update(opt_state, lri_opt, iunit)
159 : END DO
160 : ! *** ======================= END optimization =======================
161 :
162 : ! *** get final optimized parameters
163 6 : opt_state%state = 8
164 6 : CALL powell_optimize(opt_state%nvar, lri_opt%x, opt_state)
165 6 : CALL update_exponents(lri_env, lri_opt, lri_opt%x, lri_opt%zet_init, nkind)
166 :
167 : CALL write_optimized_lri_basis(lri_env, dft_section, nkind, lri_opt, &
168 6 : atomic_kind_set)
169 :
170 6 : IF (iunit > 0) THEN
171 3 : WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') opt_state%nf
172 3 : WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') opt_state%fopt
173 3 : WRITE (iunit, '(/," Printed optimized lri basis set to file")')
174 : END IF
175 :
176 : CALL cp_print_key_finished_output(iunit, logger, input, &
177 6 : "PRINT%PROGRAM_RUN_INFO")
178 :
179 6 : CALL deallocate_lri_opt(lri_opt)
180 :
181 6 : END SUBROUTINE optimize_lri_basis
182 :
183 : ! **************************************************************************************************
184 : !> \brief calculates overlap integrals (aabb) of the orbital basis set,
185 : !> required for LRI basis set optimization
186 : !> \param lri_env ...
187 : !> \param qs_env ...
188 : ! **************************************************************************************************
189 6 : SUBROUTINE calculate_lri_overlap_aabb(lri_env, qs_env)
190 :
191 : TYPE(lri_environment_type), POINTER :: lri_env
192 : TYPE(qs_environment_type), POINTER :: qs_env
193 :
194 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_lri_overlap_aabb'
195 :
196 : INTEGER :: handle, iac, iatom, ikind, ilist, jatom, &
197 : jkind, jneighbor, nba, nbb, nkind, &
198 : nlist, nneighbor
199 : REAL(KIND=dp) :: dab
200 : REAL(KIND=dp), DIMENSION(3) :: rab
201 : TYPE(cell_type), POINTER :: cell
202 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb
203 : TYPE(lri_int_rho_type), POINTER :: lriir
204 : TYPE(lri_list_type), POINTER :: lri_ints_rho
205 : TYPE(neighbor_list_iterator_p_type), &
206 6 : DIMENSION(:), POINTER :: nl_iterator
207 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
208 6 : POINTER :: soo_list
209 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
210 :
211 6 : CALL timeset(routineN, handle)
212 6 : NULLIFY (cell, lriir, lri_ints_rho, nl_iterator, obasa, obasb, &
213 6 : particle_set, soo_list)
214 :
215 6 : IF (ASSOCIATED(lri_env%soo_list)) THEN
216 6 : soo_list => lri_env%soo_list
217 :
218 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, particle_set=particle_set, &
219 6 : cell=cell)
220 :
221 6 : IF (ASSOCIATED(lri_env%lri_ints_rho)) THEN
222 0 : CALL deallocate_lri_ints_rho(lri_env%lri_ints_rho)
223 : END IF
224 :
225 6 : CALL allocate_lri_ints_rho(lri_env, lri_env%lri_ints_rho, nkind)
226 6 : lri_ints_rho => lri_env%lri_ints_rho
227 :
228 6 : CALL neighbor_list_iterator_create(nl_iterator, soo_list)
229 15 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
230 :
231 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
232 : nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor, &
233 9 : iatom=iatom, jatom=jatom, r=rab)
234 :
235 9 : iac = ikind + nkind*(jkind - 1)
236 36 : dab = SQRT(SUM(rab*rab))
237 :
238 9 : obasa => lri_env%orb_basis(ikind)%gto_basis_set
239 9 : obasb => lri_env%orb_basis(jkind)%gto_basis_set
240 9 : IF (.NOT. ASSOCIATED(obasa)) CYCLE
241 9 : IF (.NOT. ASSOCIATED(obasb)) CYCLE
242 :
243 9 : lriir => lri_ints_rho%lri_atom(iac)%lri_node(ilist)%lri_int_rho(jneighbor)
244 :
245 9 : nba = obasa%nsgf
246 9 : nbb = obasb%nsgf
247 :
248 : ! calculate integrals (aa,bb)
249 : CALL int_overlap_aabb_os(lriir%soaabb, obasa, obasb, rab, lri_env%debug, &
250 9 : lriir%dmax_aabb)
251 :
252 : END DO
253 :
254 6 : CALL neighbor_list_iterator_release(nl_iterator)
255 :
256 : END IF
257 :
258 6 : CALL timestop(handle)
259 :
260 6 : END SUBROUTINE calculate_lri_overlap_aabb
261 :
262 : ! **************************************************************************************************
263 : !> \brief initialize optimization parameter
264 : !> \param lri_env lri environment
265 : !> \param lri_opt optimization environment
266 : !> \param lri_optbas_section ...
267 : !> \param opt_state state of the optimizer
268 : !> \param x parameters to be optimized, i.e. exponents and contraction coeffs
269 : !> of the lri basis set
270 : !> \param zet_init initial values of the exponents
271 : !> \param nkind number of atom kinds
272 : !> \param iunit output unit
273 : ! **************************************************************************************************
274 6 : SUBROUTINE init_optimization(lri_env, lri_opt, lri_optbas_section, opt_state, &
275 : x, zet_init, nkind, iunit)
276 :
277 : TYPE(lri_environment_type), POINTER :: lri_env
278 : TYPE(lri_opt_type), POINTER :: lri_opt
279 : TYPE(section_vals_type), POINTER :: lri_optbas_section
280 : TYPE(opt_state_type) :: opt_state
281 : REAL(KIND=dp), DIMENSION(:), POINTER :: x, zet_init
282 : INTEGER, INTENT(IN) :: nkind, iunit
283 :
284 : INTEGER :: ikind, iset, ishell, n, nset
285 6 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
286 6 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
287 6 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc_orig
288 : TYPE(gto_basis_set_type), POINTER :: fbas
289 :
290 6 : NULLIFY (fbas, gcc_orig, npgf, nshell, zet)
291 :
292 18 : ALLOCATE (lri_opt%ri_gcc_orig(nkind))
293 :
294 : ! *** get parameters
295 : CALL get_optimization_parameter(lri_opt, lri_optbas_section, &
296 6 : opt_state)
297 :
298 6 : opt_state%nvar = 0
299 6 : opt_state%nf = 0
300 6 : opt_state%iprint = 1
301 6 : opt_state%unit = iunit
302 :
303 : ! *** init exponents
304 6 : IF (lri_opt%opt_exps) THEN
305 : n = 0
306 12 : DO ikind = 1, nkind
307 6 : fbas => lri_env%ri_basis(ikind)%gto_basis_set
308 : CALL get_gto_basis_set(gto_basis_set=fbas, &
309 6 : npgf=npgf, nset=nset, zet=zet)
310 52 : DO iset = 1, nset
311 40 : IF (lri_opt%use_geometric_seq .AND. npgf(iset) > 2) THEN
312 4 : opt_state%nvar = opt_state%nvar + 2
313 4 : CALL reallocate(x, 1, opt_state%nvar)
314 32 : x(n + 1) = MAXVAL(zet(1:npgf(iset), iset))
315 32 : x(n + 2) = MINVAL(zet(1:npgf(iset), iset))
316 4 : n = n + 2
317 : ELSE
318 36 : opt_state%nvar = opt_state%nvar + npgf(iset)
319 36 : CALL reallocate(x, 1, opt_state%nvar)
320 108 : x(n + 1:n + npgf(iset)) = zet(1:npgf(iset), iset)
321 36 : n = n + npgf(iset)
322 : END IF
323 46 : lri_opt%nexp = lri_opt%nexp + npgf(iset)
324 : END DO
325 : END DO
326 :
327 : ! *** constraints on exponents
328 6 : IF (lri_opt%use_constraints) THEN
329 6 : ALLOCATE (zet_init(SIZE(x)))
330 34 : zet_init(:) = x
331 : ELSE
332 32 : x(:) = SQRT(x)
333 : END IF
334 : END IF
335 :
336 : ! *** get the original gcc without normalization factor
337 12 : DO ikind = 1, nkind
338 6 : fbas => lri_env%ri_basis(ikind)%gto_basis_set
339 : CALL get_original_gcc(lri_opt%ri_gcc_orig(ikind)%gcc_orig, fbas, &
340 12 : lri_opt)
341 : END DO
342 :
343 : ! *** init coefficients
344 6 : IF (lri_opt%opt_coeffs) THEN
345 4 : DO ikind = 1, nkind
346 2 : fbas => lri_env%ri_basis(ikind)%gto_basis_set
347 2 : gcc_orig => lri_opt%ri_gcc_orig(ikind)%gcc_orig
348 : CALL get_gto_basis_set(gto_basis_set=fbas, &
349 2 : npgf=npgf, nset=nset, nshell=nshell, zet=zet)
350 : ! *** Gram Schmidt orthonormalization
351 2 : CALL orthonormalize_gcc(gcc_orig, fbas, lri_opt)
352 2 : n = opt_state%nvar
353 18 : DO iset = 1, nset
354 32 : DO ishell = 1, nshell(iset)
355 18 : opt_state%nvar = opt_state%nvar + npgf(iset)
356 18 : CALL reallocate(x, 1, opt_state%nvar)
357 158 : x(n + 1:n + npgf(iset)) = gcc_orig(1:npgf(iset), ishell, iset)
358 18 : lri_opt%ncoeff = lri_opt%ncoeff + npgf(iset)
359 30 : n = n + npgf(iset)
360 : END DO
361 : END DO
362 : END DO
363 : END IF
364 :
365 6 : IF (iunit > 0) THEN
366 3 : WRITE (iunit, '(/," POWELL| Accuracy",T69,ES12.5)') opt_state%rhoend
367 3 : WRITE (iunit, '(" POWELL| Initial step size",T69,ES12.5)') opt_state%rhobeg
368 : WRITE (iunit, '(" POWELL| Maximum number of evaluations",T71,I10)') &
369 3 : opt_state%maxfun
370 : WRITE (iunit, '(" POWELL| Total number of parameters",T71,I10)') &
371 3 : opt_state%nvar
372 : END IF
373 :
374 6 : END SUBROUTINE init_optimization
375 :
376 : ! **************************************************************************************************
377 : !> \brief read input for optimization
378 : !> \param lri_opt optimization environment
379 : !> \param lri_optbas_section ...
380 : !> \param opt_state state of the optimizer
381 : ! **************************************************************************************************
382 12 : SUBROUTINE get_optimization_parameter(lri_opt, lri_optbas_section, &
383 : opt_state)
384 :
385 : TYPE(lri_opt_type), POINTER :: lri_opt
386 : TYPE(section_vals_type), POINTER :: lri_optbas_section
387 : TYPE(opt_state_type) :: opt_state
388 :
389 : INTEGER :: degree_freedom
390 : TYPE(section_vals_type), POINTER :: constrain_exp_section
391 :
392 6 : NULLIFY (constrain_exp_section)
393 :
394 : ! *** parameter for POWELL optimizer
395 : CALL section_vals_val_get(lri_optbas_section, "ACCURACY", &
396 6 : r_val=opt_state%rhoend)
397 : CALL section_vals_val_get(lri_optbas_section, "STEP_SIZE", &
398 6 : r_val=opt_state%rhobeg)
399 : CALL section_vals_val_get(lri_optbas_section, "MAX_FUN", &
400 6 : i_val=opt_state%maxfun)
401 :
402 : ! *** parameters which are optimized, i.e. exps or coeff or both
403 : CALL section_vals_val_get(lri_optbas_section, "DEGREES_OF_FREEDOM", &
404 6 : i_val=degree_freedom)
405 :
406 2 : SELECT CASE (degree_freedom)
407 : CASE (do_lri_opt_all)
408 2 : lri_opt%opt_coeffs = .TRUE.
409 2 : lri_opt%opt_exps = .TRUE.
410 : CASE (do_lri_opt_coeff)
411 0 : lri_opt%opt_coeffs = .TRUE.
412 : CASE (do_lri_opt_exps)
413 4 : lri_opt%opt_exps = .TRUE.
414 : CASE DEFAULT
415 6 : CPABORT("No initialization available?????")
416 : END SELECT
417 :
418 : ! *** restraint
419 : CALL section_vals_val_get(lri_optbas_section, "USE_CONDITION_NUMBER", &
420 6 : l_val=lri_opt%use_condition_number)
421 : CALL section_vals_val_get(lri_optbas_section, "CONDITION_WEIGHT", &
422 6 : r_val=lri_opt%cond_weight)
423 : CALL section_vals_val_get(lri_optbas_section, "GEOMETRIC_SEQUENCE", &
424 6 : l_val=lri_opt%use_geometric_seq)
425 :
426 : ! *** get constraint info
427 : constrain_exp_section => section_vals_get_subs_vals(lri_optbas_section, &
428 6 : "CONSTRAIN_EXPONENTS")
429 6 : CALL section_vals_get(constrain_exp_section, explicit=lri_opt%use_constraints)
430 :
431 6 : IF (lri_opt%use_constraints) THEN
432 : CALL section_vals_val_get(constrain_exp_section, "SCALE", &
433 2 : r_val=lri_opt%scale_exp)
434 : CALL section_vals_val_get(constrain_exp_section, "FERMI_EXP", &
435 2 : r_val=lri_opt%fermi_exp)
436 : END IF
437 :
438 6 : END SUBROUTINE get_optimization_parameter
439 :
440 : ! **************************************************************************************************
441 : !> \brief update exponents after optimization step
442 : !> \param lri_env lri environment
443 : !> \param lri_opt optimization environment
444 : !> \param x optimization parameters
445 : !> \param zet_init initial values of the exponents
446 : !> \param nkind number of atomic kinds
447 : ! **************************************************************************************************
448 36 : SUBROUTINE update_exponents(lri_env, lri_opt, x, zet_init, nkind)
449 :
450 : TYPE(lri_environment_type), POINTER :: lri_env
451 : TYPE(lri_opt_type), POINTER :: lri_opt
452 : REAL(KIND=dp), DIMENSION(:), POINTER :: x, zet_init
453 : INTEGER, INTENT(IN) :: nkind
454 :
455 : INTEGER :: ikind, iset, ishell, n, nset, nvar_exp
456 36 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
457 : REAL(KIND=dp) :: zet_max, zet_min
458 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet, zet_trans
459 36 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc_orig
460 : TYPE(gto_basis_set_type), POINTER :: fbas
461 :
462 36 : NULLIFY (fbas, gcc_orig, npgf, nshell, zet_trans, zet)
463 :
464 : ! nvar_exp: number of exponents that are variables
465 36 : nvar_exp = SIZE(x) - lri_opt%ncoeff
466 108 : ALLOCATE (zet_trans(nvar_exp))
467 :
468 : ! *** update exponents
469 36 : IF (lri_opt%opt_exps) THEN
470 36 : IF (lri_opt%use_constraints) THEN
471 14 : zet => x(1:nvar_exp)
472 14 : CALL transfer_exp(lri_opt, zet, zet_init, zet_trans, nvar_exp)
473 : ELSE
474 330 : zet_trans(:) = x(1:nvar_exp)**2.0_dp
475 : END IF
476 36 : n = 0
477 72 : DO ikind = 1, nkind
478 36 : fbas => lri_env%ri_basis(ikind)%gto_basis_set
479 36 : CALL get_gto_basis_set(gto_basis_set=fbas, npgf=npgf, nset=nset)
480 310 : DO iset = 1, nset
481 274 : IF (lri_opt%use_geometric_seq .AND. npgf(iset) > 2) THEN
482 112 : zet_max = MAXVAL(zet_trans(n + 1:n + 2))
483 112 : zet_min = MINVAL(zet_trans(n + 1:n + 2))
484 28 : zet => fbas%zet(1:npgf(iset), iset)
485 28 : CALL geometric_progression(zet, zet_max, zet_min, npgf(iset))
486 28 : n = n + 2
487 : ELSE
488 630 : fbas%zet(1:npgf(iset), iset) = zet_trans(n + 1:n + npgf(iset))
489 210 : n = n + npgf(iset)
490 : END IF
491 : END DO
492 : END DO
493 : END IF
494 :
495 : ! *** update coefficients
496 36 : IF (lri_opt%opt_coeffs) THEN
497 14 : n = nvar_exp
498 28 : DO ikind = 1, nkind
499 14 : fbas => lri_env%ri_basis(ikind)%gto_basis_set
500 14 : gcc_orig => lri_opt%ri_gcc_orig(ikind)%gcc_orig
501 : CALL get_gto_basis_set(gto_basis_set=fbas, &
502 14 : nshell=nshell, npgf=npgf, nset=nset)
503 98 : DO iset = 1, nset
504 224 : DO ishell = 1, nshell(iset)
505 1106 : gcc_orig(1:npgf(iset), ishell, iset) = x(n + 1:n + npgf(iset))
506 210 : n = n + npgf(iset)
507 : END DO
508 : END DO
509 : ! *** Gram Schmidt orthonormalization
510 42 : CALL orthonormalize_gcc(gcc_orig, fbas, lri_opt)
511 : END DO
512 : END IF
513 :
514 36 : DEALLOCATE (zet_trans)
515 36 : END SUBROUTINE update_exponents
516 :
517 : ! **************************************************************************************************
518 : !> \brief employ Fermi constraint, transfer exponents
519 : !> \param lri_opt optimization environment
520 : !> \param zet untransferred exponents
521 : !> \param zet_init initial value of the exponents
522 : !> \param zet_trans transferred exponents
523 : !> \param nvar number of optimized exponents
524 : ! **************************************************************************************************
525 14 : SUBROUTINE transfer_exp(lri_opt, zet, zet_init, zet_trans, nvar)
526 :
527 : TYPE(lri_opt_type), POINTER :: lri_opt
528 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet, zet_init, zet_trans
529 : INTEGER, INTENT(IN) :: nvar
530 :
531 : REAL(KIND=dp) :: a
532 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet_max, zet_min
533 :
534 70 : ALLOCATE (zet_max(nvar), zet_min(nvar))
535 :
536 238 : zet_min(:) = zet_init(:)*(1.0_dp - lri_opt%scale_exp)
537 238 : zet_max(:) = zet_init(:)*(1.0_dp + lri_opt%scale_exp)
538 :
539 14 : a = lri_opt%fermi_exp
540 :
541 238 : zet_trans = zet_min + (zet_max - zet_min)/(1 + EXP(-a*(zet - zet_init)))
542 :
543 14 : DEALLOCATE (zet_max, zet_min)
544 :
545 14 : END SUBROUTINE transfer_exp
546 :
547 : ! **************************************************************************************************
548 : !> \brief complete geometric sequence
549 : !> \param zet all exponents of the set
550 : !> \param zet_max maximal exponent of the set
551 : !> \param zet_min minimal exponent of the set
552 : !> \param nexp number of exponents of the set
553 : ! **************************************************************************************************
554 28 : SUBROUTINE geometric_progression(zet, zet_max, zet_min, nexp)
555 :
556 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet
557 : REAL(KIND=dp), INTENT(IN) :: zet_max, zet_min
558 : INTEGER, INTENT(IN) :: nexp
559 :
560 : INTEGER :: i, n
561 : REAL(KIND=dp) :: q
562 :
563 28 : n = nexp - 1
564 :
565 28 : q = (zet_min/zet_max)**(1._dp/REAL(n, dp))
566 :
567 196 : DO i = 1, nexp
568 196 : zet(i) = zet_max*q**(i - 1)
569 : END DO
570 :
571 28 : END SUBROUTINE geometric_progression
572 :
573 : ! **************************************************************************************************
574 : !> \brief calculates the lri integrals and coefficients with the new exponents
575 : !> of the lri basis sets and calculates the objective function
576 : !> \param lri_env lri environment
577 : !> \param lri_density ...
578 : !> \param qs_env ...
579 : !> \param lri_opt optimization environment
580 : !> \param opt_state state of the optimizer
581 : !> \param pmatrix density matrix
582 : !> \param para_env ...
583 : !> \param nkind number of atomic kinds
584 : ! **************************************************************************************************
585 18 : SUBROUTINE calc_lri_integrals_get_objective(lri_env, lri_density, qs_env, &
586 : lri_opt, opt_state, pmatrix, para_env, &
587 : nkind)
588 :
589 : TYPE(lri_environment_type), POINTER :: lri_env
590 : TYPE(lri_density_type), POINTER :: lri_density
591 : TYPE(qs_environment_type), POINTER :: qs_env
592 : TYPE(lri_opt_type), POINTER :: lri_opt
593 : TYPE(opt_state_type) :: opt_state
594 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
595 : TYPE(mp_para_env_type), POINTER :: para_env
596 : INTEGER, INTENT(IN) :: nkind
597 :
598 : INTEGER :: ikind, nset
599 18 : INTEGER, DIMENSION(:), POINTER :: npgf
600 18 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
601 : TYPE(gto_basis_set_type), POINTER :: fbas
602 :
603 18 : NULLIFY (fbas, npgf)
604 :
605 : !*** build new transformation matrices sphi with new exponents
606 18 : lri_env%store_integrals = .TRUE.
607 36 : DO ikind = 1, nkind
608 18 : fbas => lri_env%ri_basis(ikind)%gto_basis_set
609 18 : CALL get_gto_basis_set(gto_basis_set=fbas, npgf=npgf, nset=nset)
610 : !build new sphi
611 3696 : fbas%gcc = lri_opt%ri_gcc_orig(ikind)%gcc_orig
612 54 : CALL init_orb_basis_set(fbas)
613 : END DO
614 18 : CALL lri_basis_init(lri_env)
615 18 : CALL calculate_lri_integrals(lri_env, qs_env)
616 18 : CALL calculate_avec_lri(lri_env, lri_density, pmatrix, cell_to_index)
617 18 : IF (lri_opt%use_condition_number) THEN
618 8 : CALL get_condition_number_of_overlap(lri_env)
619 : END IF
620 : CALL calculate_objective(lri_env, lri_density, lri_opt, pmatrix, para_env, &
621 18 : opt_state%f)
622 :
623 18 : END SUBROUTINE calc_lri_integrals_get_objective
624 :
625 : ! **************************************************************************************************
626 : !> \brief calculates the objective function defined as integral of the square
627 : !> of rhoexact - rhofit, i.e. integral[(rhoexact-rhofit)**2]
628 : !> rhoexact is the exact pair density and rhofit the lri pair density
629 : !> \param lri_env lri environment
630 : !> \param lri_density ...
631 : !> \param lri_opt optimization environment
632 : !> \param pmatrix density matrix
633 : !> \param para_env ...
634 : !> \param fobj objective function
635 : ! **************************************************************************************************
636 18 : SUBROUTINE calculate_objective(lri_env, lri_density, lri_opt, pmatrix, para_env, &
637 : fobj)
638 :
639 : TYPE(lri_environment_type), POINTER :: lri_env
640 : TYPE(lri_density_type), POINTER :: lri_density
641 : TYPE(lri_opt_type), POINTER :: lri_opt
642 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
643 : TYPE(mp_para_env_type), POINTER :: para_env
644 : REAL(KIND=dp), INTENT(OUT) :: fobj
645 :
646 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_objective'
647 :
648 : INTEGER :: handle, iac, iatom, ikind, ilist, isgfa, ispin, jatom, jkind, jneighbor, jsgfa, &
649 : ksgfb, lsgfb, mepos, nba, nbb, nfa, nfb, nkind, nlist, nn, nneighbor, nspin, nthread
650 : LOGICAL :: found, trans
651 : REAL(KIND=dp) :: obj_ab, rhoexact_sq, rhofit_sq, rhomix
652 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pbij
653 : TYPE(dbcsr_type), POINTER :: pmat
654 : TYPE(lri_int_rho_type), POINTER :: lriir
655 : TYPE(lri_int_type), POINTER :: lrii
656 : TYPE(lri_list_type), POINTER :: lri_rho
657 : TYPE(lri_rhoab_type), POINTER :: lrho
658 : TYPE(neighbor_list_iterator_p_type), &
659 18 : DIMENSION(:), POINTER :: nl_iterator
660 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
661 18 : POINTER :: soo_list
662 :
663 18 : CALL timeset(routineN, handle)
664 18 : NULLIFY (lrii, lriir, lri_rho, lrho, nl_iterator, pmat, soo_list)
665 :
666 18 : IF (ASSOCIATED(lri_env%soo_list)) THEN
667 18 : soo_list => lri_env%soo_list
668 :
669 18 : nkind = lri_env%lri_ints%nkind
670 18 : nspin = SIZE(pmatrix, 1)
671 18 : CPASSERT(SIZE(pmatrix, 2) == 1)
672 : nthread = 1
673 18 : !$ nthread = omp_get_max_threads()
674 :
675 18 : fobj = 0._dp
676 18 : lri_opt%rho_diff = 0._dp
677 :
678 54 : DO ispin = 1, nspin
679 :
680 36 : pmat => pmatrix(ispin, 1)%matrix
681 36 : lri_rho => lri_density%lri_rhos(ispin)%lri_list
682 :
683 36 : CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
684 : !$OMP PARALLEL DEFAULT(NONE)&
685 : !$OMP SHARED (nthread,nl_iterator,pmat,nkind,fobj,lri_env,lri_opt,lri_rho)&
686 : !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,nneighbor,jneighbor,&
687 : !$OMP iac,lrii,lriir,lrho,nfa,nfb,nba,nbb,nn,rhoexact_sq,rhomix,rhofit_sq,&
688 36 : !$OMP obj_ab,pbij,trans,found,isgfa,jsgfa,ksgfb,lsgfb)
689 :
690 : mepos = 0
691 : !$ mepos = omp_get_thread_num()
692 :
693 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
694 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, &
695 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor)
696 :
697 : iac = ikind + nkind*(jkind - 1)
698 :
699 : IF (.NOT. ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) CYCLE
700 :
701 : lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
702 : lriir => lri_env%lri_ints_rho%lri_atom(iac)%lri_node(ilist)%lri_int_rho(jneighbor)
703 : lrho => lri_rho%lri_atom(iac)%lri_node(ilist)%lri_rhoab(jneighbor)
704 : nfa = lrii%nfa
705 : nfb = lrii%nfb
706 : nba = lrii%nba
707 : nbb = lrii%nbb
708 : nn = nfa + nfb
709 :
710 : rhoexact_sq = 0._dp
711 : rhomix = 0._dp
712 : rhofit_sq = 0._dp
713 : obj_ab = 0._dp
714 :
715 : NULLIFY (pbij)
716 : IF (iatom <= jatom) THEN
717 : CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pbij, found=found)
718 : trans = .FALSE.
719 : ELSE
720 : CALL dbcsr_get_block_p(matrix=pmat, row=jatom, col=iatom, block=pbij, found=found)
721 : trans = .TRUE.
722 : END IF
723 : CPASSERT(found)
724 :
725 : ! *** calculate integral of the square of exact density rhoexact_sq
726 : IF (trans) THEN
727 : DO isgfa = 1, nba
728 : DO jsgfa = 1, nba
729 : DO ksgfb = 1, nbb
730 : DO lsgfb = 1, nbb
731 : rhoexact_sq = rhoexact_sq + pbij(ksgfb, isgfa)*pbij(lsgfb, jsgfa) &
732 : *lriir%soaabb(isgfa, jsgfa, ksgfb, lsgfb)
733 : END DO
734 : END DO
735 : END DO
736 : END DO
737 : ELSE
738 : DO isgfa = 1, nba
739 : DO jsgfa = 1, nba
740 : DO ksgfb = 1, nbb
741 : DO lsgfb = 1, nbb
742 : rhoexact_sq = rhoexact_sq + pbij(isgfa, ksgfb)*pbij(jsgfa, lsgfb) &
743 : *lriir%soaabb(isgfa, jsgfa, ksgfb, lsgfb)
744 : END DO
745 : END DO
746 : END DO
747 : END DO
748 : END IF
749 :
750 : ! *** calculate integral of the square of the fitted density rhofit_sq
751 : DO isgfa = 1, nfa
752 : DO jsgfa = 1, nfa
753 : rhofit_sq = rhofit_sq + lrho%avec(isgfa)*lrho%avec(jsgfa) &
754 : *lri_env%bas_prop(ikind)%ri_ovlp(isgfa, jsgfa)
755 : END DO
756 : END DO
757 : IF (iatom /= jatom) THEN
758 : DO ksgfb = 1, nfb
759 : DO lsgfb = 1, nfb
760 : rhofit_sq = rhofit_sq + lrho%avec(nfa + ksgfb)*lrho%avec(nfa + lsgfb) &
761 : *lri_env%bas_prop(jkind)%ri_ovlp(ksgfb, lsgfb)
762 : END DO
763 : END DO
764 : DO isgfa = 1, nfa
765 : DO ksgfb = 1, nfb
766 : rhofit_sq = rhofit_sq + 2._dp*lrho%avec(isgfa)*lrho%avec(nfa + ksgfb) &
767 : *lrii%sab(isgfa, ksgfb)
768 : END DO
769 : END DO
770 : END IF
771 :
772 : ! *** and integral of the product of exact and fitted density rhomix
773 : IF (iatom == jatom) THEN
774 : rhomix = SUM(lrho%avec(1:nfa)*lrho%tvec(1:nfa))
775 : ELSE
776 : rhomix = SUM(lrho%avec(1:nn)*lrho%tvec(1:nn))
777 : END IF
778 :
779 : ! *** calculate contribution to the objective function for pair ab
780 : ! *** taking density matrix symmetry in account, double-count for off-diagonal blocks
781 : IF (iatom == jatom) THEN
782 : obj_ab = rhoexact_sq - 2._dp*rhomix + rhofit_sq
783 : ELSE
784 : obj_ab = 2.0_dp*(rhoexact_sq - 2._dp*rhomix + rhofit_sq)
785 : END IF
786 :
787 : !$OMP CRITICAL(addfun)
788 : IF (lri_opt%use_condition_number) THEN
789 : fobj = fobj + obj_ab + lri_opt%cond_weight*LOG(lrii%cond_num)
790 : lri_opt%rho_diff = lri_opt%rho_diff + obj_ab
791 : ELSE
792 : fobj = fobj + obj_ab
793 : END IF
794 : !$OMP END CRITICAL(addfun)
795 :
796 : END DO
797 : !$OMP END PARALLEL
798 :
799 54 : CALL neighbor_list_iterator_release(nl_iterator)
800 :
801 : END DO
802 18 : CALL para_env%sum(fobj)
803 :
804 : END IF
805 :
806 18 : CALL timestop(handle)
807 :
808 18 : END SUBROUTINE calculate_objective
809 :
810 : ! **************************************************************************************************
811 : !> \brief get condition number of overlap matrix
812 : !> \param lri_env lri environment
813 : ! **************************************************************************************************
814 8 : SUBROUTINE get_condition_number_of_overlap(lri_env)
815 :
816 : TYPE(lri_environment_type), POINTER :: lri_env
817 :
818 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_condition_number_of_overlap'
819 :
820 : INTEGER :: handle, iac, iatom, ikind, ilist, info, &
821 : jatom, jkind, jneighbor, lwork, mepos, &
822 : nfa, nfb, nkind, nlist, nn, nneighbor, &
823 : nthread
824 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag, off_diag, tau
825 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
826 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: smat
827 : TYPE(lri_int_type), POINTER :: lrii
828 : TYPE(neighbor_list_iterator_p_type), &
829 8 : DIMENSION(:), POINTER :: nl_iterator
830 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
831 8 : POINTER :: soo_list
832 :
833 8 : CALL timeset(routineN, handle)
834 8 : NULLIFY (lrii, nl_iterator, smat, soo_list)
835 :
836 8 : soo_list => lri_env%soo_list
837 :
838 8 : nkind = lri_env%lri_ints%nkind
839 : nthread = 1
840 8 : !$ nthread = omp_get_max_threads()
841 :
842 8 : CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
843 : !$OMP PARALLEL DEFAULT(NONE)&
844 : !$OMP SHARED (nthread,nl_iterator,nkind,lri_env)&
845 : !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,nneighbor,jneighbor,&
846 8 : !$OMP diag,off_diag,smat,tau,work,iac,lrii,nfa,nfb,nn,info,lwork)
847 :
848 : mepos = 0
849 : !$ mepos = omp_get_thread_num()
850 :
851 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
852 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, &
853 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, inode=jneighbor)
854 :
855 : iac = ikind + nkind*(jkind - 1)
856 : IF (.NOT. ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) CYCLE
857 : lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
858 :
859 : nfa = lrii%nfa
860 : nfb = lrii%nfb
861 : nn = nfa + nfb
862 :
863 : ! build the overlap matrix
864 : IF (iatom == jatom) THEN
865 : ALLOCATE (smat(nfa, nfa))
866 : ELSE
867 : ALLOCATE (smat(nn, nn))
868 : END IF
869 : smat(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
870 : IF (iatom /= jatom) THEN
871 : nn = nfa + nfb
872 : smat(1:nfa, nfa + 1:nn) = lrii%sab(1:nfa, 1:nfb)
873 : smat(nfa + 1:nn, 1:nfa) = TRANSPOSE(lrii%sab(1:nfa, 1:nfb))
874 : smat(nfa + 1:nn, nfa + 1:nn) = lri_env%bas_prop(jkind)%ri_ovlp(1:nfb, 1:nfb)
875 : END IF
876 :
877 : IF (iatom == jatom) nn = nfa
878 : ALLOCATE (diag(nn), off_diag(nn - 1), tau(nn - 1), work(1))
879 : diag = 0.0_dp
880 : off_diag = 0.0_dp
881 : tau = 0.0_dp
882 : work = 0.0_dp
883 : lwork = -1
884 : ! get lwork
885 : CALL DSYTRD('U', nn, smat, nn, diag, off_diag, tau, work, lwork, info)
886 : lwork = INT(work(1))
887 : CALL reallocate(work, 1, lwork)
888 : ! get the eigenvalues
889 : CALL DSYTRD('U', nn, smat, nn, diag, off_diag, tau, work, lwork, info)
890 : CALL DSTERF(nn, diag, off_diag, info)
891 :
892 : lrii%cond_num = MAXVAL(ABS(diag))/MINVAL(ABS(diag))
893 :
894 : DEALLOCATE (diag, off_diag, smat, tau, work)
895 : END DO
896 : !$OMP END PARALLEL
897 :
898 8 : CALL neighbor_list_iterator_release(nl_iterator)
899 :
900 8 : CALL timestop(handle)
901 :
902 16 : END SUBROUTINE get_condition_number_of_overlap
903 :
904 : ! **************************************************************************************************
905 : !> \brief print recent information on optimization
906 : !> \param opt_state state of the optimizer
907 : !> \param lri_opt optimization environment
908 : !> \param iunit ...
909 : ! **************************************************************************************************
910 30 : SUBROUTINE print_optimization_update(opt_state, lri_opt, iunit)
911 :
912 : TYPE(opt_state_type) :: opt_state
913 : TYPE(lri_opt_type), POINTER :: lri_opt
914 : INTEGER, INTENT(IN) :: iunit
915 :
916 : INTEGER :: n10
917 :
918 30 : n10 = MAX(opt_state%maxfun/100, 1)
919 :
920 30 : IF (opt_state%nf == 2 .AND. opt_state%state == 2 .AND. iunit > 0) THEN
921 2 : WRITE (iunit, '(/," POWELL| Initial value of function",T61,F20.10)') opt_state%f
922 : END IF
923 30 : IF (MOD(opt_state%nf, n10) == 0 .AND. opt_state%nf > 1 .AND. iunit > 0) THEN
924 : WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
925 10 : INT(REAL(opt_state%nf, dp)/REAL(opt_state%maxfun, dp)*100._dp), opt_state%fopt
926 : END IF
927 30 : IF (lri_opt%use_condition_number) THEN
928 12 : IF (MOD(opt_state%nf, n10) == 0 .AND. opt_state%nf > 1 .AND. iunit > 0) THEN
929 : WRITE (iunit, '(" POWELL| Recent value of function without condition nr.",T61,F20.10)') &
930 5 : lri_opt%rho_diff
931 : END IF
932 : END IF
933 :
934 30 : END SUBROUTINE print_optimization_update
935 :
936 : ! **************************************************************************************************
937 : !> \brief write optimized LRI basis set to file
938 : !> \param lri_env ...
939 : !> \param dft_section ...
940 : !> \param nkind ...
941 : !> \param lri_opt ...
942 : !> \param atomic_kind_set ...
943 : ! **************************************************************************************************
944 6 : SUBROUTINE write_optimized_lri_basis(lri_env, dft_section, nkind, lri_opt, &
945 : atomic_kind_set)
946 :
947 : TYPE(lri_environment_type), POINTER :: lri_env
948 : TYPE(section_vals_type), POINTER :: dft_section
949 : INTEGER, INTENT(IN) :: nkind
950 : TYPE(lri_opt_type), POINTER :: lri_opt
951 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
952 :
953 : CHARACTER(LEN=default_path_length) :: filename
954 : INTEGER :: cc_l, ikind, ipgf, iset, ishell, nset, &
955 : output_file
956 6 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
957 6 : INTEGER, DIMENSION(:, :), POINTER :: l
958 6 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
959 6 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc_orig
960 : TYPE(cp_logger_type), POINTER :: logger
961 : TYPE(gto_basis_set_type), POINTER :: fbas
962 : TYPE(section_vals_type), POINTER :: print_key
963 :
964 6 : NULLIFY (fbas, gcc_orig, l, lmax, lmin, logger, npgf, nshell, print_key, zet)
965 :
966 : !*** do the printing
967 : print_key => section_vals_get_subs_vals(dft_section, &
968 12 : "PRINT%OPTIMIZE_LRI_BASIS")
969 6 : logger => cp_get_default_logger()
970 6 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
971 : dft_section, "PRINT%OPTIMIZE_LRI_BASIS"), &
972 : cp_p_file)) THEN
973 : output_file = cp_print_key_unit_nr(logger, dft_section, &
974 : "PRINT%OPTIMIZE_LRI_BASIS", &
975 : extension=".opt", &
976 : file_status="REPLACE", &
977 : file_action="WRITE", &
978 6 : file_form="FORMATTED")
979 :
980 6 : IF (output_file > 0) THEN
981 :
982 : filename = cp_print_key_generate_filename(logger, &
983 : print_key, extension=".opt", &
984 3 : my_local=.TRUE.)
985 :
986 6 : DO ikind = 1, nkind
987 3 : fbas => lri_env%ri_basis(ikind)%gto_basis_set
988 3 : gcc_orig => lri_opt%ri_gcc_orig(ikind)%gcc_orig
989 : CALL get_gto_basis_set(gto_basis_set=fbas, &
990 : l=l, lmax=lmax, lmin=lmin, &
991 : npgf=npgf, nshell=nshell, &
992 3 : nset=nset, zet=zet)
993 3 : WRITE (output_file, '(T1,A2,T5,A)') TRIM(atomic_kind_set(ikind)%name), &
994 6 : TRIM(fbas%name)
995 3 : WRITE (output_file, '(T1,I4)') nset
996 29 : DO iset = 1, nset
997 20 : WRITE (output_file, '(4(1X,I0))', advance='no') 2, lmin(iset), &
998 40 : lmax(iset), npgf(iset)
999 20 : cc_l = 1
1000 65 : DO ishell = 1, nshell(iset)
1001 65 : IF (ishell /= nshell(iset)) THEN
1002 25 : IF (l(ishell, iset) == l(ishell + 1, iset)) THEN
1003 3 : cc_l = cc_l + 1
1004 : ELSE
1005 22 : WRITE (output_file, '(1X,I0)', advance='no') cc_l
1006 22 : cc_l = 1
1007 : END IF
1008 : ELSE
1009 20 : WRITE (output_file, '(1X,I0)') cc_l
1010 : END IF
1011 : END DO
1012 53 : DO ipgf = 1, npgf(iset)
1013 30 : WRITE (output_file, '(F18.12)', advance='no') zet(ipgf, iset)
1014 121 : DO ishell = 1, nshell(iset)
1015 101 : IF (ishell == nshell(iset)) THEN
1016 30 : WRITE (output_file, '(T5,F18.12)') gcc_orig(ipgf, ishell, iset)
1017 : ELSE
1018 41 : WRITE (output_file, '(T5,F18.12)', advance='no') gcc_orig(ipgf, ishell, iset)
1019 : END IF
1020 : END DO
1021 : END DO
1022 : END DO
1023 : END DO
1024 :
1025 : END IF
1026 :
1027 : CALL cp_print_key_finished_output(output_file, logger, dft_section, &
1028 6 : "PRINT%OPTIMIZE_LRI_BASIS")
1029 : END IF
1030 :
1031 6 : END SUBROUTINE write_optimized_lri_basis
1032 :
1033 : END MODULE lri_optimize_ri_basis
|
