(git:ccc2433)
mp2_optimize_ri_basis.F
Go to the documentation of this file.
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 Routines to optimize the RI-MP2 basis. Only exponents of
10 !> non-contracted auxiliary basis basis are optimized. The derivative
11 !> of the MP2 energy with respect to the exponents of the basis
12 !> are calculated numerically.
13 !> \par History
14 !> 08.2013 created [Mauro Del Ben]
15 !> \author Mauro Del Ben
16 ! **************************************************************************************************
17 MODULE mp2_optimize_ri_basis
18  USE atomic_kind_types, ONLY: atomic_kind_type
19  USE basis_set_container_types, ONLY: add_basis_set_to_container,&
20  remove_basis_from_container
21  USE basis_set_types, ONLY: allocate_gto_basis_set,&
22  gto_basis_set_type
23  USE cp_log_handling, ONLY: cp_add_default_logger,&
24  cp_get_default_logger,&
25  cp_logger_create,&
26  cp_logger_get_default_unit_nr,&
27  cp_logger_release,&
28  cp_logger_set,&
29  cp_logger_type,&
30  cp_rm_default_logger
31  USE hfx_types, ONLY: hfx_basis_info_type,&
32  hfx_basis_type
33  USE kinds, ONLY: dp
34  USE libint_wrapper, ONLY: build_eri_size
35  USE machine, ONLY: default_output_unit,&
36  m_flush
37  USE memory_utilities, ONLY: reallocate
38  USE message_passing, ONLY: mp_para_env_release,&
39  mp_para_env_type
40  USE mp2_direct_method, ONLY: mp2_canonical_direct_single_batch
41  USE mp2_ri_libint, ONLY: libint_ri_mp2,&
42  read_ri_basis_set,&
43  release_ri_basis_set
44  USE mp2_types, ONLY: mp2_biel_type,&
45  mp2_type
46  USE orbital_pointers, ONLY: indco,&
47  init_orbital_pointers,&
48  nco,&
49  ncoset,&
50  nso
51  USE orbital_symbols, ONLY: cgf_symbol,&
52  sgf_symbol
53  USE qs_environment_types, ONLY: get_qs_env,&
54  qs_environment_type
55  USE qs_kind_types, ONLY: get_qs_kind,&
56  qs_kind_type
57  USE util, ONLY: sort
58 #include "./base/base_uses.f90"
59 
60  IMPLICIT NONE
61 
62  PRIVATE
63 
64  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_optimize_ri_basis'
65 
66  PUBLIC :: optimize_ri_basis_main
67 
68 CONTAINS
69 
70 ! **************************************************************************************************
71 !> \brief optimize RI-MP2 basis set
72 !> \param Emp2 ...
73 !> \param Emp2_Cou ...
74 !> \param Emp2_ex ...
75 !> \param Emp2_S ...
76 !> \param Emp2_T ...
77 !> \param dimen ...
78 !> \param natom ...
79 !> \param homo ...
80 !> \param mp2_biel ...
81 !> \param mp2_env ...
82 !> \param C ...
83 !> \param Auto ...
84 !> \param kind_of ...
85 !> \param qs_env ...
86 !> \param para_env ...
87 !> \param unit_nr ...
88 !> \param homo_beta ...
89 !> \param C_beta ...
90 !> \param Auto_beta ...
91 !> \author Mauro Del Ben
92 ! **************************************************************************************************
93  SUBROUTINE optimize_ri_basis_main(Emp2, Emp2_Cou, Emp2_ex, Emp2_S, Emp2_T, dimen, natom, homo, &
94  mp2_biel, mp2_env, C, Auto, kind_of, &
95  qs_env, para_env, &
96  unit_nr, homo_beta, C_beta, Auto_beta)
97 
98  REAL(kind=dp), INTENT(OUT) :: emp2, emp2_cou, emp2_ex, emp2_s, emp2_t
99  INTEGER, INTENT(IN) :: dimen, natom, homo
100  TYPE(mp2_biel_type), INTENT(IN) :: mp2_biel
101  TYPE(mp2_type), INTENT(INOUT) :: mp2_env
102  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: c
103  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: auto
104  INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
105  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
106  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
107  INTEGER, INTENT(IN) :: unit_nr
108  INTEGER, INTENT(IN), OPTIONAL :: homo_beta
109  REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
110  OPTIONAL :: c_beta
111  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: auto_beta
112 
113  CHARACTER(len=*), PARAMETER :: routinen = 'optimize_ri_basis_main'
114 
115  INTEGER :: color_sub, dimen_ri, elements_ij_proc, handle, i, iiter, ikind, ipgf, iset, &
116  ishell, j, local_unit_nr, max_l_am, max_num_iter, ndof, nkind, number_groups, virtual, &
117  virtual_beta
118  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ij_list_proc, index_table_ri
119  LOGICAL :: hess_up, open_shell_case, reset_boundary
120  LOGICAL, ALLOCATABLE, DIMENSION(:) :: basis_was_assoc
121  REAL(kind=dp) :: di, di_new, dri, dri_new, emp2_aa, emp2_aa_cou, emp2_aa_ex, emp2_ab, &
122  emp2_ab_cou, emp2_ab_ex, emp2_bb, emp2_bb_cou, emp2_bb_ex, emp2_ri, emp2_ri_new, &
123  eps_di_rel, eps_dri, eps_step, expon, fac, fad, fae, reset_edge, sumdg, sumxi
124  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: deriv, dg, g, hdg, lower_b, max_dev, &
125  max_rel_dev, p, pnew, xi
126  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: exp_limits, hessin
127  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: integ_mp2, integ_mp2_aa, integ_mp2_ab, &
128  integ_mp2_bb
129  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
130  TYPE(cp_logger_type), POINTER :: logger, logger_sub
131  TYPE(gto_basis_set_type), POINTER :: ri_aux_basis
132  TYPE(hfx_basis_info_type) :: ri_basis_info
133  TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_s0, ri_basis_parameter
134  TYPE(mp_para_env_type), POINTER :: para_env_sub
135  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
136 
137  CALL timeset(routinen, handle)
138  logger => cp_get_default_logger()
139 
140  open_shell_case = .false.
141  IF (PRESENT(homo_beta) .AND. PRESENT(c_beta) .AND. PRESENT(auto_beta)) open_shell_case = .true.
142 
143  virtual = dimen - homo
144 
145  eps_dri = mp2_env%ri_opt_param%DRI
146  eps_di_rel = mp2_env%ri_opt_param%DI_rel
147  eps_step = mp2_env%ri_opt_param%eps_step
148  max_num_iter = mp2_env%ri_opt_param%max_num_iter
149 
150  ! calculate the ERI's over molecular integrals
151  emp2 = 0.0_dp
152  emp2_cou = 0.0_dp
153  emp2_ex = 0.0_dp
154  emp2_s = 0.0_dp
155  emp2_t = 0.0_dp
156  IF (open_shell_case) THEN
157  ! open shell case
158  virtual_beta = dimen - homo_beta
159 
160  ! alpha-aplha case
161  emp2_aa = 0.0_dp
162  emp2_aa_cou = 0.0_dp
163  emp2_aa_ex = 0.0_dp
164  CALL calc_elem_ij_proc(homo, homo, para_env, elements_ij_proc, ij_list_proc)
165  CALL mp2_canonical_direct_single_batch(emp2_aa, emp2_aa_cou, emp2_aa_ex, mp2_env, qs_env, para_env, &
166  mp2_biel, dimen, c, auto, 0, homo, homo, &
167  elements_ij_proc, ij_list_proc, homo, 0, &
168  integ_mp2=integ_mp2_aa)
169  CALL para_env%sum(emp2_aa_cou)
170  CALL para_env%sum(emp2_aa_ex)
171  CALL para_env%sum(emp2_aa)
172  DEALLOCATE (ij_list_proc)
173 
174  ! beta-beta case
175  emp2_bb = 0.0_dp
176  emp2_bb_cou = 0.0_dp
177  emp2_bb_ex = 0.0_dp
178  CALL calc_elem_ij_proc(homo_beta, homo_beta, para_env, elements_ij_proc, ij_list_proc)
179  CALL mp2_canonical_direct_single_batch(emp2_bb, emp2_bb_cou, emp2_bb_ex, mp2_env, qs_env, para_env, &
180  mp2_biel, dimen, c_beta, auto_beta, 0, homo_beta, homo_beta, &
181  elements_ij_proc, ij_list_proc, homo_beta, 0, &
182  integ_mp2=integ_mp2_bb)
183  CALL para_env%sum(emp2_bb_cou)
184  CALL para_env%sum(emp2_bb_ex)
185  CALL para_env%sum(emp2_bb)
186  DEALLOCATE (ij_list_proc)
187 
188  ! aplha-beta case
189  emp2_ab = 0.0_dp
190  emp2_ab_cou = 0.0_dp
191  emp2_ab_ex = 0.0_dp
192  CALL calc_elem_ij_proc(homo, homo_beta, para_env, elements_ij_proc, ij_list_proc)
193  CALL mp2_canonical_direct_single_batch(emp2_ab, emp2_ab_cou, emp2_ab_ex, mp2_env, qs_env, para_env, &
194  mp2_biel, dimen, c, auto, 0, homo, homo, &
195  elements_ij_proc, ij_list_proc, homo_beta, 0, &
196  homo_beta, c_beta, auto_beta, integ_mp2=integ_mp2_ab)
197  CALL para_env%sum(emp2_ab_cou)
198  CALL para_env%sum(emp2_ab_ex)
199  CALL para_env%sum(emp2_ab)
200  DEALLOCATE (ij_list_proc)
201 
202  emp2 = emp2_aa + emp2_bb + emp2_ab*2.0_dp !+Emp2_BA
203  emp2_cou = emp2_aa_cou + emp2_bb_cou + emp2_ab_cou*2.0_dp !+Emp2_BA
204  emp2_ex = emp2_aa_ex + emp2_bb_ex + emp2_ab_ex*2.0_dp !+Emp2_BA
205 
206  emp2_s = emp2_ab*2.0_dp
207  emp2_t = emp2_aa + emp2_bb
208 
209  ! Replicate the MO-ERI's over all processes
210  CALL para_env%sum(integ_mp2_aa)
211  CALL para_env%sum(integ_mp2_bb)
212  CALL para_env%sum(integ_mp2_ab)
213 
214  ELSE
215  ! close shell case
216  CALL calc_elem_ij_proc(homo, homo, para_env, elements_ij_proc, ij_list_proc)
217  CALL mp2_canonical_direct_single_batch(emp2, emp2_cou, emp2_ex, mp2_env, qs_env, para_env, &
218  mp2_biel, dimen, c, auto, 0, homo, homo, &
219  elements_ij_proc, ij_list_proc, homo, 0, &
220  integ_mp2=integ_mp2)
221  CALL para_env%sum(emp2_cou)
222  CALL para_env%sum(emp2_ex)
223  CALL para_env%sum(emp2)
224  DEALLOCATE (ij_list_proc)
225 
226  ! Replicate the MO-ERI's over all processes
227  CALL para_env%sum(integ_mp2)
228 
229  END IF
230 
231  ! create the para_env_sub
232  number_groups = para_env%num_pe/mp2_env%mp2_num_proc
233  color_sub = para_env%mepos/mp2_env%mp2_num_proc
234  ALLOCATE (para_env_sub)
235  CALL para_env_sub%from_split(para_env, color_sub)
236 
237  IF (para_env%is_source()) THEN
238  local_unit_nr = cp_logger_get_default_unit_nr(logger, local=.false.)
239  ELSE
240  local_unit_nr = default_output_unit
241  END IF
242  NULLIFY (logger_sub)
243  CALL cp_logger_create(logger_sub, para_env=para_env_sub, &
244  default_global_unit_nr=local_unit_nr, close_global_unit_on_dealloc=.false.)
245  CALL cp_logger_set(logger_sub, local_filename="opt_RI_basis_localLog")
246  CALL cp_add_default_logger(logger_sub)
247 
248  CALL generate_ri_init_basis(qs_env, mp2_env, nkind, max_rel_dev, basis_was_assoc)
249 
250  CALL read_ri_basis_set(qs_env, ri_basis_parameter, ri_basis_info, &
251  natom, nkind, kind_of, index_table_ri, dimen_ri, &
252  basis_s0)
253 
254  ndof = 0
255  max_l_am = 0
256  DO ikind = 1, nkind
257  DO iset = 1, ri_basis_parameter(ikind)%nset
258  ndof = ndof + 1
259  max_l_am = max(max_l_am, maxval(ri_basis_parameter(ikind)%lmax))
260  END DO
261  END DO
262 
263  ALLOCATE (exp_limits(2, nkind))
264  exp_limits(1, :) = huge(0)
265  exp_limits(2, :) = -huge(0)
266  DO ikind = 1, nkind
267  DO iset = 1, ri_basis_parameter(ikind)%nset
268  expon = ri_basis_parameter(ikind)%zet(1, iset)
269  IF (expon <= exp_limits(1, ikind)) exp_limits(1, ikind) = expon
270  IF (expon >= exp_limits(2, ikind)) exp_limits(2, ikind) = expon
271  END DO
272  IF (basis_was_assoc(ikind)) THEN
273  exp_limits(1, ikind) = exp_limits(1, ikind)*0.5_dp
274  exp_limits(2, ikind) = exp_limits(2, ikind)*1.5_dp
275  ELSE
276  exp_limits(1, ikind) = exp_limits(1, ikind)*0.8_dp*0.5_dp
277  exp_limits(2, ikind) = exp_limits(2, ikind)*1.2_dp*1.5_dp
278  END IF
279  END DO
280  DEALLOCATE (basis_was_assoc)
281 
282  ! check if the max angular momentum exceed the libint one
283  IF (max_l_am > build_eri_size) THEN
284  cpabort("The angular momentum needed exceeds the value assumed when configuring libint.")
285  END IF
286 
287  ! Allocate stuff
288  ALLOCATE (p(ndof))
289  p = 0.0_dp
290  ALLOCATE (xi(ndof))
291  xi = 0.0_dp
292  ALLOCATE (g(ndof))
293  g = 0.0_dp
294  ALLOCATE (dg(ndof))
295  dg = 0.0_dp
296  ALLOCATE (hdg(ndof))
297  hdg = 0.0_dp
298  ALLOCATE (pnew(ndof))
299  pnew = 0.0_dp
300  ALLOCATE (deriv(ndof))
301  deriv = 0.0_dp
302  ALLOCATE (hessin(ndof, ndof))
303  hessin = 0.0_dp
304  DO i = 1, ndof
305  hessin(i, i) = 1.0_dp
306  END DO
307 
308  ! initialize transformation function
309  ALLOCATE (lower_b(ndof))
310  lower_b = 0.0_dp
311  ALLOCATE (max_dev(ndof))
312  max_dev = 0.0_dp
313 
314  ! Initialize the transformation function
315  CALL init_transf(nkind, ri_basis_parameter, lower_b, max_dev, max_rel_dev)
316 
317  ! get the atomic kind set for writing the basis
318  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
319 
320  ! Calculate RI-MO-ERI's
321  CALL calc_energy_func(emp2, emp2_aa, emp2_bb, emp2_ab, emp2_ri, dri, di, &
322  integ_mp2, integ_mp2_aa, integ_mp2_bb, integ_mp2_ab, &
323  qs_env, natom, dimen, dimen_ri, homo, virtual, &
324  kind_of, index_table_ri, mp2_biel, mp2_env, auto, c, &
325  ri_basis_parameter, ri_basis_info, basis_s0, &
326  open_shell_case, homo_beta, virtual_beta, auto_beta, c_beta, para_env, unit_nr, &
327  .true.)
328 
329  ! ! Calculate function (DI) derivatives with respect to the RI basis exponent
330  CALL calc_energy_func_der(emp2, emp2_aa, emp2_bb, emp2_ab, di, &
331  integ_mp2, integ_mp2_aa, integ_mp2_bb, integ_mp2_ab, eps_step, &
332  qs_env, nkind, natom, dimen, dimen_ri, homo, virtual, &
333  kind_of, index_table_ri, mp2_biel, mp2_env, auto, c, &
334  ri_basis_parameter, ri_basis_info, basis_s0, &
335  open_shell_case, homo_beta, virtual_beta, auto_beta, c_beta, &
336  para_env, para_env_sub, number_groups, color_sub, unit_nr, &
337  p, lower_b, max_dev, &
338  deriv)
339 
340  g(:) = deriv
341  xi(:) = -g
342 
343  reset_edge = 1.5_dp
344  DO iiter = 1, max_num_iter
345  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,I5)') 'OPTIMIZATION STEP NUMBER', iiter
346 
347  ! perform step
348  pnew(:) = p + xi
349  CALL p2basis(nkind, ri_basis_parameter, lower_b, max_dev, pnew)
350 
351  ! check if we have to reset boundaries
352  reset_boundary = .false.
353  i = 0
354  DO ikind = 1, nkind
355  DO iset = 1, ri_basis_parameter(ikind)%nset
356  i = i + 1
357  expon = transf_val(lower_b(i), max_dev(i), pnew(i))
358  IF (abs(pnew(i)) > reset_edge .OR. expon < exp_limits(1, ikind) .OR. expon > exp_limits(2, ikind)) THEN
359  reset_boundary = .true.
360  EXIT
361  END IF
362  END DO
363  END DO
364  ! IF(nreset>1) reset_boundary=.TRUE.
365 
366  IF (reset_boundary) THEN
367  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'RESET BASIS: one of the exponent hits the boundary'
368  CALL reset_basis(nkind, ndof, ri_basis_parameter, reset_edge, &
369  pnew, lower_b, max_dev, max_rel_dev, exp_limits)
370  p(:) = pnew
371  xi = 0.0_dp
372  g = 0.0_dp
373  dg = 0.0_dp
374  hdg = 0.0_dp
375  pnew = 0.0_dp
376  hessin = 0.0_dp
377  DO i = 1, ndof
378  hessin(i, i) = 1.0_dp
379  END DO
380  deriv = 0.0_dp
381  ! Calculate RI-MO-ERI's
382  CALL calc_energy_func(emp2, emp2_aa, emp2_bb, emp2_ab, emp2_ri, dri, di, &
383  integ_mp2, integ_mp2_aa, integ_mp2_bb, integ_mp2_ab, &
384  qs_env, natom, dimen, dimen_ri, homo, virtual, &
385  kind_of, index_table_ri, mp2_biel, mp2_env, auto, c, &
386  ri_basis_parameter, ri_basis_info, basis_s0, &
387  open_shell_case, homo_beta, virtual_beta, auto_beta, c_beta, para_env, unit_nr, &
388  .false.)
389  ! ! Calculate function (DI) derivatives with respect to the RI basis exponent
390  CALL calc_energy_func_der(emp2, emp2_aa, emp2_bb, emp2_ab, di, &
391  integ_mp2, integ_mp2_aa, integ_mp2_bb, integ_mp2_ab, eps_step, &
392  qs_env, nkind, natom, dimen, dimen_ri, homo, virtual, &
393  kind_of, index_table_ri, mp2_biel, mp2_env, auto, c, &
394  ri_basis_parameter, ri_basis_info, basis_s0, &
395  open_shell_case, homo_beta, virtual_beta, auto_beta, c_beta, &
396  para_env, para_env_sub, number_groups, color_sub, unit_nr, &
397  p, lower_b, max_dev, &
398  deriv)
399 
400  g(:) = deriv
401  xi(:) = -g
402  pnew(:) = p + xi
403  CALL p2basis(nkind, ri_basis_parameter, lower_b, max_dev, pnew)
404  END IF
405 
406  ! calculate energy at the new point
407  CALL calc_energy_func(emp2, emp2_aa, emp2_bb, emp2_ab, emp2_ri_new, dri_new, di_new, &
408  integ_mp2, integ_mp2_aa, integ_mp2_bb, integ_mp2_ab, &
409  qs_env, natom, dimen, dimen_ri, homo, virtual, &
410  kind_of, index_table_ri, mp2_biel, mp2_env, auto, c, &
411  ri_basis_parameter, ri_basis_info, basis_s0, &
412  open_shell_case, homo_beta, virtual_beta, auto_beta, c_beta, para_env, unit_nr, &
413  .false.)
414 
415  ! update energy and direction
416  di = di_new
417  xi(:) = pnew - p
418  p(:) = pnew
419 
420  ! check for convergence
421  IF (unit_nr > 0) THEN
422  WRITE (unit_nr, *)
423  DO ikind = 1, nkind
424  CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=ri_aux_basis, &
425  basis_type="RI_AUX")
426  WRITE (unit_nr, '(T3,A,A)') atomic_kind_set(ikind)%element_symbol, ' RI_opt_basis'
427  WRITE (unit_nr, '(T3,I3)') ri_basis_parameter(ikind)%nset
428  DO iset = 1, ri_basis_parameter(ikind)%nset
429  WRITE (unit_nr, '(T3,10I4)') iset, &
430  ri_basis_parameter(ikind)%lmin(iset), &
431  ri_basis_parameter(ikind)%lmax(iset), &
432  ri_basis_parameter(ikind)%npgf(iset), &
433  (1, ishell=1, ri_basis_parameter(ikind)%nshell(iset))
434  DO ipgf = 1, ri_basis_parameter(ikind)%npgf(iset)
435  WRITE (unit_nr, '(T3,10F16.10)') ri_basis_parameter(ikind)%zet(ipgf, iset), &
436  (ri_aux_basis%gcc(ipgf, ishell, iset), &
437  ishell=1, ri_aux_basis%nshell(iset))
438  END DO
439  END DO
440  WRITE (unit_nr, *)
441  END DO
442  WRITE (unit_nr, *)
443  CALL m_flush(unit_nr)
444  END IF
445  IF (di/abs(emp2) <= eps_di_rel .AND. abs(dri_new) <= eps_dri) THEN
446  IF (unit_nr > 0) THEN
447  WRITE (unit_nr, '(T3,A,/)') 'OPTIMIZATION CONVERGED'
448  CALL m_flush(unit_nr)
449  END IF
450  EXIT
451  END IF
452 
453  ! calculate gradients
454  CALL calc_energy_func_der(emp2, emp2_aa, emp2_bb, emp2_ab, di, &
455  integ_mp2, integ_mp2_aa, integ_mp2_bb, integ_mp2_ab, eps_step, &
456  qs_env, nkind, natom, dimen, dimen_ri, homo, virtual, &
457  kind_of, index_table_ri, mp2_biel, mp2_env, auto, c, &
458  ri_basis_parameter, ri_basis_info, basis_s0, &
459  open_shell_case, homo_beta, virtual_beta, auto_beta, c_beta, &
460  para_env, para_env_sub, number_groups, color_sub, unit_nr, &
461  p, lower_b, max_dev, &
462  deriv)
463 
464  ! g is the vector containing the old gradient
465  dg(:) = deriv - g
466  g(:) = deriv
467  hdg(:) = matmul(hessin, dg)
468 
469  fac = sum(dg*xi)
470  fae = sum(dg*hdg)
471  sumdg = sum(dg*dg)
472  sumxi = sum(xi*xi)
473 
474  hess_up = .true.
475  IF (fac**2 > sumdg*sumxi*3.0e-8_dp) THEN
476  fac = 1.0_dp/fac
477  fad = 1.0_dp/fae
478  dg(:) = fac*xi - fad*hdg
479  DO i = 1, ndof
480  DO j = 1, ndof
481  hessin(i, j) = hessin(i, j) + fac*xi(i)*xi(j) &
482  - fad*hdg(i)*hdg(j) &
483  + fae*dg(i)*dg(j)
484  END DO
485  END DO
486  ELSE
487  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Skip Hessian Update'
488  hess_up = .false.
489  END IF
490 
491  ! new direction
492  xi(:) = -matmul(hessin, g)
493 
494  END DO
495 
496  IF (.NOT. (di/abs(emp2) <= eps_di_rel .AND. abs(dri_new) <= eps_dri)) THEN
497  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,I5,A)') 'OPTIMIZATION NOT CONVERGED IN', max_num_iter, ' STEPS.'
498  IF (unit_nr > 0) WRITE (unit_nr, *)
499  END IF
500 
501  DEALLOCATE (max_rel_dev)
502 
503  DEALLOCATE (p)
504  DEALLOCATE (xi)
505  DEALLOCATE (g)
506  DEALLOCATE (pnew)
507  DEALLOCATE (dg)
508  DEALLOCATE (hdg)
509  DEALLOCATE (deriv)
510  DEALLOCATE (hessin)
511  DEALLOCATE (lower_b)
512  DEALLOCATE (max_dev)
513  DEALLOCATE (exp_limits)
514 
515  IF (open_shell_case) THEN
516  DEALLOCATE (integ_mp2_aa)
517  DEALLOCATE (integ_mp2_bb)
518  DEALLOCATE (integ_mp2_ab)
519  ELSE
520  DEALLOCATE (integ_mp2)
521  END IF
522  DEALLOCATE (index_table_ri)
523 
524  ! Release RI basis set
525  CALL release_ri_basis_set(ri_basis_parameter, basis_s0)
526 
527  CALL mp_para_env_release(para_env_sub)
528  CALL cp_rm_default_logger()
529  CALL cp_logger_release(logger_sub)
530 
531  CALL timestop(handle)
532 
533  END SUBROUTINE optimize_ri_basis_main
534 
535 ! **************************************************************************************************
536 !> \brief ...
537 !> \param Emp2 ...
538 !> \param Emp2_AA ...
539 !> \param Emp2_BB ...
540 !> \param Emp2_AB ...
541 !> \param DI_ref ...
542 !> \param Integ_MP2 ...
543 !> \param Integ_MP2_AA ...
544 !> \param Integ_MP2_BB ...
545 !> \param Integ_MP2_AB ...
546 !> \param eps ...
547 !> \param qs_env ...
548 !> \param nkind ...
549 !> \param natom ...
550 !> \param dimen ...
551 !> \param dimen_RI ...
552 !> \param homo ...
553 !> \param virtual ...
554 !> \param kind_of ...
555 !> \param index_table_RI ...
556 !> \param mp2_biel ...
557 !> \param mp2_env ...
558 !> \param Auto ...
559 !> \param C ...
560 !> \param RI_basis_parameter ...
561 !> \param RI_basis_info ...
562 !> \param basis_S0 ...
563 !> \param open_shell_case ...
564 !> \param homo_beta ...
565 !> \param virtual_beta ...
566 !> \param Auto_beta ...
567 !> \param C_beta ...
568 !> \param para_env ...
569 !> \param para_env_sub ...
570 !> \param number_groups ...
571 !> \param color_sub ...
572 !> \param unit_nr ...
573 !> \param p ...
574 !> \param lower_B ...
575 !> \param max_dev ...
576 !> \param deriv ...
577 ! **************************************************************************************************
578  SUBROUTINE calc_energy_func_der(Emp2, Emp2_AA, Emp2_BB, Emp2_AB, DI_ref, &
579  Integ_MP2, Integ_MP2_AA, Integ_MP2_BB, Integ_MP2_AB, eps, &
580  qs_env, nkind, natom, dimen, dimen_RI, homo, virtual, &
581  kind_of, index_table_RI, mp2_biel, mp2_env, Auto, C, &
582  RI_basis_parameter, RI_basis_info, basis_S0, &
583  open_shell_case, homo_beta, virtual_beta, Auto_beta, C_beta, &
584  para_env, para_env_sub, number_groups, color_sub, unit_nr, &
585  p, lower_B, max_dev, &
586  deriv)
587  REAL(kind=dp), INTENT(IN) :: emp2
588  REAL(kind=dp), INTENT(INOUT) :: emp2_aa, emp2_bb, emp2_ab
589  REAL(kind=dp), INTENT(IN) :: di_ref
590  REAL(kind=dp), ALLOCATABLE, &
591  DIMENSION(:, :, :, :), INTENT(IN) :: integ_mp2, integ_mp2_aa, integ_mp2_bb, &
592  integ_mp2_ab
593  REAL(kind=dp), INTENT(IN) :: eps
594  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
595  INTEGER, INTENT(IN) :: nkind, natom, dimen, dimen_ri, homo, &
596  virtual
597  INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
598  INTEGER, ALLOCATABLE, DIMENSION(:, :), &
599  INTENT(INOUT) :: index_table_ri
600  TYPE(mp2_biel_type), INTENT(IN) :: mp2_biel
601  TYPE(mp2_type), INTENT(INOUT) :: mp2_env
602  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: auto
603  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: c
604  TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
605  POINTER :: ri_basis_parameter
606  TYPE(hfx_basis_info_type), INTENT(IN) :: ri_basis_info
607  TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
608  POINTER :: basis_s0
609  LOGICAL, INTENT(IN) :: open_shell_case
610  INTEGER, INTENT(IN) :: homo_beta, virtual_beta
611  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: auto_beta
612  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: c_beta
613  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env, para_env_sub
614  INTEGER, INTENT(IN) :: number_groups, color_sub, unit_nr
615  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: p, lower_b, max_dev
616  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: deriv
617 
618  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_energy_func_der'
619 
620  INTEGER :: handle, ideriv, ikind, iset, nseta
621  REAL(kind=dp) :: di, dri, emp2_ri, new_basis_val, &
622  orig_basis_val
623  REAL(kind=dp), VOLATILE :: step, temp
624 
625  CALL timeset(routinen, handle)
626 
627  step = eps
628 
629  ! cycle over the RI basis set exponent
630  deriv = 0.0_dp
631  ideriv = 0
632  DO ikind = 1, nkind
633  nseta = ri_basis_parameter(ikind)%nset
634  DO iset = 1, nseta
635  ! for now only uncontracted aux basis set
636  ideriv = ideriv + 1
637  IF (mod(ideriv, number_groups) /= color_sub) cycle
638 
639  ! calculate the numerical derivative
640  ! The eps is the relative change of the exponent for the
641  ! calculation of the numerical derivative
642 
643  ! in the new case eps is just the step length for calculating the numerical derivative
644 
645  cpassert(ri_basis_parameter(ikind)%npgf(iset) == 1)
646  orig_basis_val = ri_basis_parameter(ikind)%zet(1, iset)
647  temp = p(ideriv) + step
648  new_basis_val = transf_val(lower_b(ideriv), max_dev(ideriv), temp)
649  ri_basis_parameter(ikind)%zet(1, iset) = new_basis_val
650 
651  CALL calc_energy_func(emp2, emp2_aa, emp2_bb, emp2_ab, emp2_ri, dri, di, &
652  integ_mp2, integ_mp2_aa, integ_mp2_bb, integ_mp2_ab, &
653  qs_env, natom, dimen, dimen_ri, homo, virtual, &
654  kind_of, index_table_ri, mp2_biel, mp2_env, auto, c, &
655  ri_basis_parameter, ri_basis_info, basis_s0, &
656  open_shell_case, homo_beta, virtual_beta, auto_beta, c_beta, &
657  para_env_sub, unit_nr, .true.)
658 
659  ri_basis_parameter(ikind)%zet(1, iset) = orig_basis_val
660 
661  IF (para_env_sub%mepos == 0) THEN
662  temp = exp(di)
663  temp = temp/exp(di_ref)
664  deriv(ideriv) = log(temp)/step
665  END IF
666 
667  END DO
668  END DO
669 
670  CALL para_env%sum(deriv)
671 
672  CALL timestop(handle)
673 
674  END SUBROUTINE
675 
676 ! **************************************************************************************************
677 !> \brief ...
678 !> \param Emp2 ...
679 !> \param Emp2_AA ...
680 !> \param Emp2_BB ...
681 !> \param Emp2_AB ...
682 !> \param Emp2_RI ...
683 !> \param DRI ...
684 !> \param DI ...
685 !> \param Integ_MP2 ...
686 !> \param Integ_MP2_AA ...
687 !> \param Integ_MP2_BB ...
688 !> \param Integ_MP2_AB ...
689 !> \param qs_env ...
690 !> \param natom ...
691 !> \param dimen ...
692 !> \param dimen_RI ...
693 !> \param homo ...
694 !> \param virtual ...
695 !> \param kind_of ...
696 !> \param index_table_RI ...
697 !> \param mp2_biel ...
698 !> \param mp2_env ...
699 !> \param Auto ...
700 !> \param C ...
701 !> \param RI_basis_parameter ...
702 !> \param RI_basis_info ...
703 !> \param basis_S0 ...
704 !> \param open_shell_case ...
705 !> \param homo_beta ...
706 !> \param virtual_beta ...
707 !> \param Auto_beta ...
708 !> \param C_beta ...
709 !> \param para_env ...
710 !> \param unit_nr ...
711 !> \param no_write ...
712 ! **************************************************************************************************
713  SUBROUTINE calc_energy_func(Emp2, Emp2_AA, Emp2_BB, Emp2_AB, Emp2_RI, DRI, DI, &
714  Integ_MP2, Integ_MP2_AA, Integ_MP2_BB, Integ_MP2_AB, &
715  qs_env, natom, dimen, dimen_RI, homo, virtual, &
716  kind_of, index_table_RI, mp2_biel, mp2_env, Auto, C, &
717  RI_basis_parameter, RI_basis_info, basis_S0, &
718  open_shell_case, homo_beta, virtual_beta, Auto_beta, C_beta, para_env, unit_nr, &
719  no_write)
720  REAL(kind=dp), INTENT(IN) :: emp2, emp2_aa, emp2_bb, emp2_ab
721  REAL(kind=dp), INTENT(OUT) :: emp2_ri, dri, di
722  REAL(kind=dp), ALLOCATABLE, &
723  DIMENSION(:, :, :, :), INTENT(IN) :: integ_mp2, integ_mp2_aa, integ_mp2_bb, &
724  integ_mp2_ab
725  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
726  INTEGER, INTENT(IN) :: natom, dimen, dimen_ri, homo, virtual
727  INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
728  INTEGER, ALLOCATABLE, DIMENSION(:, :), &
729  INTENT(INOUT) :: index_table_ri
730  TYPE(mp2_biel_type), INTENT(IN) :: mp2_biel
731  TYPE(mp2_type), INTENT(INOUT) :: mp2_env
732  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: auto
733  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: c
734  TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
735  POINTER :: ri_basis_parameter
736  TYPE(hfx_basis_info_type), INTENT(IN) :: ri_basis_info
737  TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
738  POINTER :: basis_s0
739  LOGICAL, INTENT(IN) :: open_shell_case
740  INTEGER, INTENT(IN) :: homo_beta, virtual_beta
741  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: auto_beta
742  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: c_beta
743  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
744  INTEGER, INTENT(IN) :: unit_nr
745  LOGICAL, INTENT(IN) :: no_write
746 
747  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_energy_func'
748 
749  INTEGER :: handle
750  REAL(kind=dp) :: di_aa, di_ab, di_bb, dri_aa, dri_ab, &
751  dri_bb, emp2_ri_aa, emp2_ri_ab, &
752  emp2_ri_bb
753  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: lai, lai_beta
754 
755  CALL timeset(routinen, handle)
756 
757  CALL libint_ri_mp2(dimen, dimen_ri, homo, natom, mp2_biel, mp2_env, c, &
758  kind_of, ri_basis_parameter, ri_basis_info, basis_s0, index_table_ri, &
759  qs_env, para_env, lai)
760  IF (open_shell_case) THEN
761  CALL libint_ri_mp2(dimen, dimen_ri, homo_beta, natom, mp2_biel, mp2_env, c_beta, &
762  kind_of, ri_basis_parameter, ri_basis_info, basis_s0, index_table_ri, &
763  qs_env, para_env, lai_beta)
764  END IF
765 
766  ! Contract integrals into energy
767  IF (open_shell_case) THEN
768  ! alpha-alpha
769  CALL contract_integrals(di_aa, emp2_ri_aa, dri_aa, emp2_aa, homo, homo, virtual, virtual, &
770  1.0_dp, 0.5_dp, .true., &
771  auto, auto, integ_mp2_aa, &
772  lai, lai, para_env)
773 
774  ! beta-beta
775  CALL contract_integrals(di_bb, emp2_ri_bb, dri_bb, emp2_bb, homo_beta, homo_beta, virtual_beta, virtual_beta, &
776  1.0_dp, 0.5_dp, .true., &
777  auto_beta, auto_beta, integ_mp2_bb, &
778  lai_beta, lai_beta, para_env)
779 
780  ! alpha-beta
781  CALL contract_integrals(di_ab, emp2_ri_ab, dri_ab, emp2_ab*2.0_dp, homo, homo_beta, virtual, virtual_beta, &
782  1.0_dp, 1.0_dp, .false., &
783  auto, auto_beta, integ_mp2_ab, &
784  lai, lai_beta, para_env)
785 
786  emp2_ri = emp2_ri_aa + emp2_ri_bb + emp2_ri_ab
787  dri = dri_aa + dri_bb + dri_ab
788  di = di_aa + di_bb + di_ab
789  ELSE
790  CALL contract_integrals(di, emp2_ri, dri, emp2, homo, homo, virtual, virtual, &
791  2.0_dp, 1.0_dp, .true., &
792  auto, auto, integ_mp2, &
793  lai, lai, para_env)
794  END IF
795 
796  IF (.NOT. no_write) THEN
797  IF (unit_nr > 0) THEN
798  WRITE (unit_nr, '(/,(T3,A,T56,F25.14))') &
799  'Emp2 =', emp2, &
800  'Emp2-RI =', emp2_ri
801  WRITE (unit_nr, '(T3,A,T56,ES25.10)') &
802  'DRI =', dri, &
803  'DI =', di, &
804  'DI/|Emp2| =', di/abs(emp2)
805  CALL m_flush(unit_nr)
806  END IF
807  END IF
808 
809  DEALLOCATE (lai)
810  IF (open_shell_case) DEALLOCATE (lai_beta)
811 
812  CALL timestop(handle)
813 
814  END SUBROUTINE
815 
816 ! **************************************************************************************************
817 !> \brief ...
818 !> \param nkind ...
819 !> \param RI_basis_parameter ...
820 !> \param lower_B ...
821 !> \param max_dev ...
822 !> \param max_rel_dev ...
823 ! **************************************************************************************************
824  PURE SUBROUTINE init_transf(nkind, RI_basis_parameter, lower_B, max_dev, max_rel_dev)
825  INTEGER, INTENT(IN) :: nkind
826  TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
827  POINTER :: ri_basis_parameter
828  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: lower_b, max_dev
829  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: max_rel_dev
830 
831  INTEGER :: ikind, ipos, iset
832 
833  ipos = 0
834  DO ikind = 1, nkind
835  DO iset = 1, ri_basis_parameter(ikind)%nset
836  ipos = ipos + 1
837  lower_b(ipos) = ri_basis_parameter(ikind)%zet(1, iset)*(1.0_dp - max_rel_dev(ipos))
838  max_dev(ipos) = ri_basis_parameter(ikind)%zet(1, iset)*2.0_dp*max_rel_dev(ipos)
839  END DO
840  END DO
841 
842  END SUBROUTINE
843 
844 ! **************************************************************************************************
845 !> \brief ...
846 !> \param nkind ...
847 !> \param RI_basis_parameter ...
848 !> \param Lower_B ...
849 !> \param max_dev ...
850 !> \param p ...
851 ! **************************************************************************************************
852  SUBROUTINE p2basis(nkind, RI_basis_parameter, Lower_B, max_dev, p)
853  INTEGER, INTENT(IN) :: nkind
854  TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
855  POINTER :: ri_basis_parameter
856  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
857  INTENT(IN) :: lower_b, max_dev, p
858 
859  INTEGER :: ikind, ipos, iset
860  REAL(kind=dp) :: valout
861 
862  ipos = 0
863  DO ikind = 1, nkind
864  DO iset = 1, ri_basis_parameter(ikind)%nset
865  ipos = ipos + 1
866  valout = transf_val(lower_b(ipos), max_dev(ipos), p(ipos))
867  ri_basis_parameter(ikind)%zet(1, iset) = valout
868  END DO
869  END DO
870 
871  END SUBROUTINE
872 
873 ! **************************************************************************************************
874 !> \brief ...
875 !> \param nkind ...
876 !> \param ndof ...
877 !> \param RI_basis_parameter ...
878 !> \param reset_edge ...
879 !> \param pnew ...
880 !> \param lower_B ...
881 !> \param max_dev ...
882 !> \param max_rel_dev ...
883 !> \param exp_limits ...
884 ! **************************************************************************************************
885  SUBROUTINE reset_basis(nkind, ndof, RI_basis_parameter, reset_edge, &
886  pnew, lower_B, max_dev, max_rel_dev, exp_limits)
887  INTEGER, INTENT(IN) :: nkind, ndof
888  TYPE(hfx_basis_type), DIMENSION(:), INTENT(IN), &
889  POINTER :: ri_basis_parameter
890  REAL(kind=dp), INTENT(IN) :: reset_edge
891  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: pnew, lower_b, max_dev, max_rel_dev
892  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: exp_limits
893 
894  INTEGER :: am_max, iexpo, ikind, ipos, ipos_p, &
895  iset, la
896  INTEGER, ALLOCATABLE, DIMENSION(:) :: nf_per_l
897  INTEGER, DIMENSION(ndof) :: change_expo
898  LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_to_be_changed
899  REAL(kind=dp) :: expo, geom_fact, pmax
900  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: max_min_exp_per_l
901  REAL(kind=dp), DIMENSION(ndof) :: new_expo, old_expo, old_lower_b, &
902  old_max_dev, old_max_rel_dev, old_pnew
903 
904 ! make a copy of the original parameters
905 
906  old_pnew = pnew
907  old_lower_b = lower_b
908  old_max_dev = max_dev
909  old_max_rel_dev = max_rel_dev
910  old_expo = 0.0_dp
911  ipos = 0
912  DO ikind = 1, nkind
913  DO iset = 1, ri_basis_parameter(ikind)%nset
914  ipos = ipos + 1
915  old_expo(ipos) = ri_basis_parameter(ikind)%zet(1, iset)
916  END DO
917  END DO
918 
919  pnew = 0.0_dp
920  lower_b = 0.0_dp
921  max_dev = 0.0_dp
922  max_rel_dev = 0.0_dp
923 
924  change_expo = 0
925 
926  new_expo = 0.0_dp
927  ipos = 0
928  ipos_p = 0
929  DO ikind = 1, nkind
930  am_max = maxval(ri_basis_parameter(ikind)%lmax(:))
931  ALLOCATE (nf_per_l(0:am_max))
932  nf_per_l = 0
933  ALLOCATE (max_min_exp_per_l(2, 0:am_max))
934  max_min_exp_per_l(1, :) = huge(0)
935  max_min_exp_per_l(2, :) = -huge(0)
936 
937  DO iset = 1, ri_basis_parameter(ikind)%nset
938  la = ri_basis_parameter(ikind)%lmax(iset)
939  expo = ri_basis_parameter(ikind)%zet(1, iset)
940  nf_per_l(la) = nf_per_l(la) + 1
941  IF (expo <= max_min_exp_per_l(1, la)) max_min_exp_per_l(1, la) = expo
942  IF (expo >= max_min_exp_per_l(2, la)) max_min_exp_per_l(2, la) = expo
943  END DO
944 
945  max_min_exp_per_l(1, la) = max(max_min_exp_per_l(1, la), exp_limits(1, ikind))
946  max_min_exp_per_l(2, la) = min(max_min_exp_per_l(2, la), exp_limits(2, ikind))
947 
948  ! always s exponents as maximum and minimu
949  ! expo=MAXVAL(max_min_exp_per_l(2,:))
950  ! max_min_exp_per_l(2,0)=expo
951  ! expo=MINVAL(max_min_exp_per_l(1,:))
952  ! max_min_exp_per_l(1,0)=expo
953 
954  ALLOCATE (has_to_be_changed(0:am_max))
955  has_to_be_changed = .false.
956  DO la = 0, am_max
957  pmax = -huge(0)
958  DO iexpo = 1, nf_per_l(la)
959  ipos_p = ipos_p + 1
960  IF (abs(old_pnew(ipos_p)) >= pmax) pmax = abs(old_pnew(ipos_p))
961  ! check if any of the exponents go out of range
962  expo = transf_val(old_lower_b(ipos_p), old_max_dev(ipos_p), old_pnew(ipos_p))
963  IF (expo < exp_limits(1, ikind) .OR. expo > exp_limits(2, ikind)) has_to_be_changed(la) = .true.
964  END DO
965  IF (pmax > reset_edge) has_to_be_changed(la) = .true.
966  END DO
967 
968  DO la = 0, am_max
969  IF (nf_per_l(la) == 1) THEN
970  ipos = ipos + 1
971  new_expo(ipos) = max_min_exp_per_l(1, la)
972  IF (new_expo(ipos) >= exp_limits(1, ikind) .AND. new_expo(ipos) <= exp_limits(2, ikind)) THEN
973  max_rel_dev(ipos) = (new_expo(ipos) - old_lower_b(ipos))/new_expo(ipos)
974  IF (max_rel_dev(ipos) <= 0.1_dp) max_rel_dev(ipos) = 0.8_dp
975  ELSE
976  new_expo(ipos) = (exp_limits(2, ikind) + exp_limits(1, ikind))/2.0_dp
977  max_rel_dev(ipos) = (new_expo(ipos) - exp_limits(1, ikind))/new_expo(ipos)
978  END IF
979  IF (has_to_be_changed(la)) change_expo(ipos) = 1
980  ELSE
981  IF (max_min_exp_per_l(2, la)/max_min_exp_per_l(1, la) < 2.0_dp) THEN
982  max_min_exp_per_l(1, la) = max_min_exp_per_l(1, la)*0.5
983  max_min_exp_per_l(2, la) = max_min_exp_per_l(2, la)*1.5
984  END IF
985  geom_fact = (max_min_exp_per_l(2, la)/max_min_exp_per_l(1, la))**(1.0_dp/real(nf_per_l(la) - 1, dp))
986  DO iexpo = 1, nf_per_l(la)
987  ipos = ipos + 1
988  new_expo(ipos) = max_min_exp_per_l(1, la)*(geom_fact**(iexpo - 1))
989  max_rel_dev(ipos) = (geom_fact - 1.0_dp)/(geom_fact + 1.0_dp)*0.9_dp
990  IF (has_to_be_changed(la)) change_expo(ipos) = 1
991  END DO
992  END IF
993  END DO
994 
995  DEALLOCATE (has_to_be_changed)
996 
997  DEALLOCATE (nf_per_l)
998  DEALLOCATE (max_min_exp_per_l)
999  END DO
1000 
1001  ipos = 0
1002  DO ikind = 1, nkind
1003  DO iset = 1, ri_basis_parameter(ikind)%nset
1004  ipos = ipos + 1
1005  ri_basis_parameter(ikind)%zet(1, iset) = new_expo(ipos)
1006  END DO
1007  END DO
1008 
1009  CALL init_transf(nkind, ri_basis_parameter, lower_b, max_dev, max_rel_dev)
1010 
1011  ipos = 0
1012  DO ikind = 1, nkind
1013  DO iset = 1, ri_basis_parameter(ikind)%nset
1014  ipos = ipos + 1
1015  IF (change_expo(ipos) == 0) THEN
1016  ! restore original
1017  pnew(ipos) = old_pnew(ipos)
1018  lower_b(ipos) = old_lower_b(ipos)
1019  max_dev(ipos) = old_max_dev(ipos)
1020  max_rel_dev(ipos) = old_max_rel_dev(ipos)
1021  ri_basis_parameter(ikind)%zet(1, iset) = old_expo(ipos)
1022  END IF
1023  END DO
1024  END DO
1025 
1026  END SUBROUTINE
1027 
1028 ! **************************************************************************************************
1029 !> \brief ...
1030 !> \param DI ...
1031 !> \param Emp2_RI ...
1032 !> \param DRI ...
1033 !> \param Emp2 ...
1034 !> \param homo ...
1035 !> \param homo_beta ...
1036 !> \param virtual ...
1037 !> \param virtual_beta ...
1038 !> \param fact ...
1039 !> \param fact2 ...
1040 !> \param calc_ex ...
1041 !> \param MOenerg ...
1042 !> \param MOenerg_beta ...
1043 !> \param abij ...
1044 !> \param Lai ...
1045 !> \param Lai_beta ...
1046 !> \param para_env ...
1047 ! **************************************************************************************************
1048  SUBROUTINE contract_integrals(DI, Emp2_RI, DRI, Emp2, homo, homo_beta, virtual, virtual_beta, &
1049  fact, fact2, calc_ex, &
1050  MOenerg, MOenerg_beta, abij, &
1051  Lai, Lai_beta, para_env)
1052  REAL(kind=dp), INTENT(OUT) :: di, emp2_ri, dri
1053  REAL(kind=dp), INTENT(IN) :: emp2
1054  INTEGER, INTENT(IN) :: homo, homo_beta, virtual, virtual_beta
1055  REAL(kind=dp), INTENT(IN) :: fact, fact2
1056  LOGICAL, INTENT(IN) :: calc_ex
1057  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: moenerg, moenerg_beta
1058  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: abij
1059  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: lai, lai_beta
1060  TYPE(mp_para_env_type), INTENT(IN) :: para_env
1061 
1062  INTEGER :: a, b, i, ij_counter, j
1063  REAL(kind=dp) :: t_iajb, t_iajb_ri
1064  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_ab
1065 
1066  ALLOCATE (mat_ab(virtual, virtual_beta))
1067 
1068  di = 0.0_dp
1069  emp2_ri = 0.0_dp
1070  ij_counter = 0
1071  DO j = 1, homo_beta
1072  DO i = 1, homo
1073  ij_counter = ij_counter + 1
1074  IF (mod(ij_counter, para_env%num_pe) /= para_env%mepos) cycle
1075  mat_ab = 0.0_dp
1076  mat_ab(:, :) = matmul(transpose(lai(:, :, i)), lai_beta(:, :, j))
1077  DO b = 1, virtual_beta
1078  DO a = 1, virtual
1079  IF (calc_ex) THEN
1080  t_iajb = fact*abij(a, b, i, j) - abij(b, a, i, j)
1081  t_iajb_ri = fact*mat_ab(a, b) - mat_ab(b, a)
1082  ELSE
1083  t_iajb = fact*abij(a, b, i, j)
1084  t_iajb_ri = fact*mat_ab(a, b)
1085  END IF
1086  t_iajb = t_iajb/(moenerg(a + homo) + moenerg_beta(b + homo_beta) - moenerg(i) - moenerg_beta(j))
1087  t_iajb_ri = t_iajb_ri/(moenerg(a + homo) + moenerg_beta(b + homo_beta) - moenerg(i) - moenerg_beta(j))
1088 
1089  emp2_ri = emp2_ri - t_iajb_ri*mat_ab(a, b)*fact2
1090 
1091  di = di - t_iajb*mat_ab(a, b)*fact2
1092 
1093  END DO
1094  END DO
1095  END DO
1096  END DO
1097  CALL para_env%sum(di)
1098  CALL para_env%sum(emp2_ri)
1099 
1100  dri = emp2 - emp2_ri
1101  di = 2.0d+00*di - emp2 - emp2_ri
1102 
1103  DEALLOCATE (mat_ab)
1104 
1105  END SUBROUTINE
1106 
1107 ! **************************************************************************************************
1108 !> \brief ...
1109 !> \param homo ...
1110 !> \param homo_beta ...
1111 !> \param para_env ...
1112 !> \param elements_ij_proc ...
1113 !> \param ij_list_proc ...
1114 ! **************************************************************************************************
1115  PURE SUBROUTINE calc_elem_ij_proc(homo, homo_beta, para_env, elements_ij_proc, ij_list_proc)
1116  INTEGER, INTENT(IN) :: homo, homo_beta
1117  TYPE(mp_para_env_type), INTENT(IN) :: para_env
1118  INTEGER, INTENT(OUT) :: elements_ij_proc
1119  INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: ij_list_proc
1120 
1121  INTEGER :: i, ij_counter, j
1122 
1123  elements_ij_proc = 0
1124  ij_counter = -1
1125  DO i = 1, homo
1126  DO j = 1, homo_beta
1127  ij_counter = ij_counter + 1
1128  IF (mod(ij_counter, para_env%num_pe) == para_env%mepos) elements_ij_proc = elements_ij_proc + 1
1129  END DO
1130  END DO
1131 
1132  ALLOCATE (ij_list_proc(elements_ij_proc, 2))
1133  ij_list_proc = 0
1134  ij_counter = -1
1135  elements_ij_proc = 0
1136  DO i = 1, homo
1137  DO j = 1, homo_beta
1138  ij_counter = ij_counter + 1
1139  IF (mod(ij_counter, para_env%num_pe) == para_env%mepos) THEN
1140  elements_ij_proc = elements_ij_proc + 1
1141  ij_list_proc(elements_ij_proc, 1) = i
1142  ij_list_proc(elements_ij_proc, 2) = j
1143  END IF
1144  END DO
1145  END DO
1146 
1147  END SUBROUTINE calc_elem_ij_proc
1148 
1149 ! **************************************************************************************************
1150 !> \brief ...
1151 !> \param lower_B ...
1152 !> \param max_dev ...
1153 !> \param valin ...
1154 !> \return ...
1155 ! **************************************************************************************************
1156  ELEMENTAL FUNCTION transf_val(lower_B, max_dev, valin) RESULT(valout)
1157  REAL(kind=dp), INTENT(IN) :: lower_b, max_dev, valin
1158  REAL(kind=dp) :: valout
1159 
1160  REAL(kind=dp), PARAMETER :: alpha = 2.633915794_dp
1161 
1162  valout = 0.0_dp
1163  valout = lower_b + max_dev/(1.0_dp + exp(-alpha*valin))
1164 
1165  END FUNCTION
1166 
1167 ! **************************************************************************************************
1168 !> \brief ...
1169 !> \param qs_env ...
1170 !> \param mp2_env ...
1171 !> \param nkind ...
1172 !> \param max_rel_dev_output ...
1173 !> \param basis_was_assoc ...
1174 ! **************************************************************************************************
1175  SUBROUTINE generate_ri_init_basis(qs_env, mp2_env, nkind, max_rel_dev_output, basis_was_assoc)
1176  TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1177  TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1178  INTEGER, INTENT(OUT) :: nkind
1179  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1180  INTENT(INOUT) :: max_rel_dev_output
1181  LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: basis_was_assoc
1182 
1183  CHARACTER(LEN=*), PARAMETER :: routinen = 'generate_RI_init_basis'
1184 
1185  INTEGER :: basis_quality, handle, iexpo, iii, ikind, ipgf, iset, ishell, jexpo, jjj, la, &
1186  max_am, nexpo_shell, nseta, nsgfa, nsgfa_ri, prog_func, prog_l, ri_max_am, ri_nset, &
1187  ri_prev_size
1188  INTEGER, ALLOCATABLE, DIMENSION(:) :: l_expo, num_sgf_per_l, ordered_pos, &
1189  ri_l_expo, ri_num_sgf_per_l, &
1190  tot_num_exp_per_l
1191  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: l_tab
1192  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell
1193  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl
1194  LOGICAL :: external_num_of_func
1195  REAL(dp) :: geom_fact
1196  REAL(dp), ALLOCATABLE, DIMENSION(:) :: exponents, ri_exponents
1197  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: exp_tab, max_min_exp_l
1198  REAL(dp), DIMENSION(:, :), POINTER :: sphi_a, zet
1199  REAL(dp), DIMENSION(:, :, :), POINTER :: gcc
1200  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: max_rel_dev, max_rel_dev_prev
1201  TYPE(gto_basis_set_type), POINTER :: orb_basis_a, tmp_basis
1202  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1203  TYPE(qs_kind_type), POINTER :: atom_kind
1204 
1205  CALL timeset(routinen, handle)
1206 
1207  basis_quality = mp2_env%ri_opt_param%basis_quality
1208  external_num_of_func = .false.
1209  IF (ALLOCATED(mp2_env%ri_opt_param%RI_nset_per_l)) external_num_of_func = .true.
1210 
1211  NULLIFY (qs_kind_set)
1212  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
1213  nkind = SIZE(qs_kind_set, 1)
1214 
1215  ALLOCATE (basis_was_assoc(nkind))
1216  basis_was_assoc = .false.
1217 
1218  IF (external_num_of_func .AND. nkind > 1) THEN
1219  CALL cp_warn(__location__, &
1220  "There are more than one kind of atom. The same pattern of functions, "// &
1221  "as specified by NUM_FUNC, will be used for all kinds.")
1222  END IF
1223 
1224  DO ikind = 1, nkind
1225  NULLIFY (atom_kind)
1226  atom_kind => qs_kind_set(ikind)
1227 
1228  CALL get_qs_kind(qs_kind=atom_kind, basis_set=orb_basis_a, basis_type="RI_AUX")
1229  ! save info if the basis was or not associated
1230  basis_was_assoc(ikind) = ASSOCIATED(orb_basis_a)
1231 
1232  ! skip the generation of the initial guess if the RI basis is
1233  ! provided in input
1234  IF (basis_was_assoc(ikind)) THEN
1235  nseta = orb_basis_a%nset
1236  la_max => orb_basis_a%lmax
1237  la_min => orb_basis_a%lmin
1238  npgfa => orb_basis_a%npgf
1239  nshell => orb_basis_a%nshell
1240  zet => orb_basis_a%zet
1241  nl => orb_basis_a%l
1242 
1243  max_am = maxval(la_max)
1244 
1245  ri_max_am = max_am
1246  ALLOCATE (ri_num_sgf_per_l(0:ri_max_am))
1247  ri_num_sgf_per_l = 0
1248  ri_nset = 0
1249  DO iset = 1, nseta
1250  DO la = la_min(iset), la_max(iset)
1251  ri_nset = ri_nset + 1
1252  ri_num_sgf_per_l(la) = ri_num_sgf_per_l(la) + 1
1253  IF (npgfa(iset) > 1) THEN
1254  CALL cp_warn(__location__, &
1255  "The RI basis set optimizer can not handle contracted Gaussian. "// &
1256  "Calculation continue with only uncontracted functions.")
1257  END IF
1258  END DO
1259  END DO
1260 
1261  ALLOCATE (exp_tab(maxval(ri_num_sgf_per_l), 0:ri_max_am))
1262  exp_tab = 0.0_dp
1263  DO iii = 0, ri_max_am
1264  iexpo = 0
1265  DO iset = 1, nseta
1266  DO la = la_min(iset), la_max(iset)
1267  IF (la /= iii) cycle
1268  iexpo = iexpo + 1
1269  exp_tab(iexpo, iii) = zet(1, iset)
1270  END DO
1271  END DO
1272  END DO
1273 
1274  ! sort exponents
1275  DO iii = 0, ri_max_am
1276  ALLOCATE (ordered_pos(ri_num_sgf_per_l(iii)))
1277  ordered_pos = 0
1278  CALL sort(exp_tab(1:ri_num_sgf_per_l(iii), iii), ri_num_sgf_per_l(iii), ordered_pos)
1279  DEALLOCATE (ordered_pos)
1280  END DO
1281 
1282  ALLOCATE (ri_l_expo(ri_nset))
1283  ri_l_expo = 0
1284  ALLOCATE (ri_exponents(ri_nset))
1285  ri_exponents = 0.0_dp
1286 
1287  iset = 0
1288  DO iii = 0, ri_max_am
1289  DO iexpo = 1, ri_num_sgf_per_l(iii)
1290  iset = iset + 1
1291  ri_l_expo(iset) = iii
1292  ri_exponents(iset) = exp_tab(iexpo, iii)
1293  END DO
1294  END DO
1295  DEALLOCATE (exp_tab)
1296 
1297  ALLOCATE (max_rel_dev(ri_nset))
1298  max_rel_dev = 1.0_dp
1299  iset = 0
1300  DO iii = 0, ri_max_am
1301  IF (ri_num_sgf_per_l(iii) == 1) THEN
1302  iset = iset + 1
1303  max_rel_dev(iset) = 0.35_dp
1304  ELSE
1305  iset = iset + 1
1306  max_rel_dev(iset) = (ri_exponents(iset + 1) + ri_exponents(iset))/2.0_dp
1307  max_rel_dev(iset) = max_rel_dev(iset)/ri_exponents(iset) - 1.0_dp
1308  DO iexpo = 2, ri_num_sgf_per_l(iii)
1309  iset = iset + 1
1310  max_rel_dev(iset) = (ri_exponents(iset) + ri_exponents(iset - 1))/2.0_dp
1311  max_rel_dev(iset) = 1.0_dp - max_rel_dev(iset)/ri_exponents(iset)
1312  END DO
1313  END IF
1314  END DO
1315  max_rel_dev(:) = max_rel_dev(:)*0.9_dp
1316  DEALLOCATE (ri_num_sgf_per_l)
1317 
1318  ! deallocate the old basis before moving on
1319  CALL remove_basis_from_container(qs_kind_set(ikind)%basis_sets, basis_type="RI_AUX")
1320 
1321  ELSE
1322 
1323  CALL get_qs_kind(qs_kind=atom_kind, basis_set=orb_basis_a)
1324 
1325  sphi_a => orb_basis_a%sphi
1326  nseta = orb_basis_a%nset
1327  la_max => orb_basis_a%lmax
1328  la_min => orb_basis_a%lmin
1329  npgfa => orb_basis_a%npgf
1330  first_sgfa => orb_basis_a%first_sgf
1331  nshell => orb_basis_a%nshell
1332  zet => orb_basis_a%zet
1333  gcc => orb_basis_a%gcc
1334  nl => orb_basis_a%l
1335  nsgfa = orb_basis_a%nsgf
1336 
1337  max_am = maxval(la_max)
1338 
1339  nexpo_shell = 0
1340  DO iset = 1, nseta
1341  DO ishell = 1, nshell(iset)
1342  DO la = la_min(iset), la_max(iset)
1343  IF (ishell > 1) THEN
1344  IF (nl(ishell, iset) == nl(ishell - 1, iset)) cycle
1345  END IF
1346  IF (la /= nl(ishell, iset)) cycle
1347  nexpo_shell = nexpo_shell + npgfa(iset)
1348  END DO
1349  END DO
1350  END DO
1351 
1352  ALLOCATE (exponents(nexpo_shell))
1353  exponents = 0.0_dp
1354  ALLOCATE (l_expo(nexpo_shell))
1355  l_expo = 0
1356  ALLOCATE (num_sgf_per_l(0:max_am))
1357  num_sgf_per_l = 0
1358  iexpo = 0
1359  DO iset = 1, nseta
1360  DO ishell = 1, nshell(iset)
1361  DO la = la_min(iset), la_max(iset)
1362  IF (ishell > 1) THEN
1363  IF (nl(ishell, iset) == nl(ishell - 1, iset)) cycle
1364  END IF
1365  IF (la /= nl(ishell, iset)) cycle
1366  DO ipgf = 1, npgfa(iset)
1367  iexpo = iexpo + 1
1368  exponents(iexpo) = zet(ipgf, iset)
1369  l_expo(iexpo) = la
1370  END DO
1371  num_sgf_per_l(la) = num_sgf_per_l(la) + 1
1372  END DO
1373  END DO
1374  END DO
1375 
1376  ALLOCATE (exp_tab(nexpo_shell, nexpo_shell))
1377  exp_tab = 0.0_dp
1378  ALLOCATE (l_tab(nexpo_shell, nexpo_shell))
1379  l_tab = 0
1380  ALLOCATE (tot_num_exp_per_l(0:max_am*2))
1381  tot_num_exp_per_l = 0
1382  DO iexpo = 1, nexpo_shell
1383  DO jexpo = iexpo, nexpo_shell
1384  exp_tab(jexpo, iexpo) = exponents(jexpo) + exponents(iexpo)
1385  exp_tab(iexpo, jexpo) = exp_tab(jexpo, iexpo)
1386  l_tab(jexpo, iexpo) = l_expo(jexpo) + l_expo(iexpo)
1387  l_tab(iexpo, jexpo) = l_tab(jexpo, iexpo)
1388  tot_num_exp_per_l(l_tab(jexpo, iexpo)) = tot_num_exp_per_l(l_tab(jexpo, iexpo)) + 1
1389  END DO
1390  END DO
1391  DEALLOCATE (l_expo)
1392  DEALLOCATE (exponents)
1393 
1394  ALLOCATE (max_min_exp_l(2, 0:max_am*2))
1395  max_min_exp_l(1, :) = huge(0)
1396  max_min_exp_l(2, :) = -huge(0)
1397 
1398  DO la = 0, max_am*2
1399  DO iexpo = 1, nexpo_shell
1400  DO jexpo = iexpo, nexpo_shell
1401  IF (la /= l_tab(jexpo, iexpo)) cycle
1402  IF (exp_tab(jexpo, iexpo) <= max_min_exp_l(1, la)) max_min_exp_l(1, la) = exp_tab(jexpo, iexpo)
1403  IF (exp_tab(jexpo, iexpo) >= max_min_exp_l(2, la)) max_min_exp_l(2, la) = exp_tab(jexpo, iexpo)
1404  END DO
1405  END DO
1406  END DO
1407  DEALLOCATE (exp_tab)
1408  DEALLOCATE (l_tab)
1409 
1410  ! ! scale the limits of the exponents to avoid the optimizer to go out-of-range
1411  ! ! (0.2 just empirical parameter)
1412  max_min_exp_l(1, :) = max_min_exp_l(1, :)/0.80_dp
1413  max_min_exp_l(2, :) = max_min_exp_l(2, :)/1.20_dp
1414 
1415  ALLOCATE (ri_num_sgf_per_l(0:max_am*2))
1416  ri_num_sgf_per_l = 0
1417 
1418  SELECT CASE (basis_quality)
1419  CASE (0)
1420  ! normal
1421  prog_func = 0
1422  prog_l = 2
1423  CASE (1)
1424  ! large
1425  prog_func = 1
1426  prog_l = 2
1427  CASE (2)
1428  prog_func = 2
1429  prog_l = 2
1430  CASE DEFAULT
1431  prog_func = 0
1432  prog_l = 2
1433  END SELECT
1434 
1435  IF (external_num_of_func) THEN
1436  ! cp2k can not exceed angular momentum 7
1437  ri_max_am = min(SIZE(mp2_env%ri_opt_param%RI_nset_per_l) - 1, 7)
1438  IF (ri_max_am > max_am*2) THEN
1439  DEALLOCATE (ri_num_sgf_per_l)
1440  ALLOCATE (ri_num_sgf_per_l(0:ri_max_am))
1441  ri_num_sgf_per_l = 0
1442  END IF
1443  DO la = 0, ri_max_am
1444  ri_num_sgf_per_l(la) = mp2_env%ri_opt_param%RI_nset_per_l(la)
1445  END DO
1446  ELSE
1447  ri_num_sgf_per_l(0) = num_sgf_per_l(0)*2 + prog_func
1448  DO la = 1, max_am
1449  ri_num_sgf_per_l(la) = ri_num_sgf_per_l(la - 1) - 1
1450  IF (ri_num_sgf_per_l(la) == 0) THEN
1451  ri_num_sgf_per_l(la) = 1
1452  EXIT
1453  END IF
1454  END DO
1455  DO la = max_am + 1, max_am*2
1456  ri_num_sgf_per_l(la) = ri_num_sgf_per_l(la - 1) - prog_l
1457  IF (ri_num_sgf_per_l(la) == 0) THEN
1458  ri_num_sgf_per_l(la) = 1
1459  END IF
1460  IF (ri_num_sgf_per_l(la) == 1) EXIT
1461  END DO
1462  ri_max_am = min(max_am*2, 7)
1463  DO la = 0, min(max_am*2, 7)
1464  IF (ri_num_sgf_per_l(la) == 0) THEN
1465  ri_max_am = la - 1
1466  EXIT
1467  END IF
1468  END DO
1469 
1470  iii = 0
1471  jjj = 0
1472  nsgfa_ri = 0
1473  DO la = 1, max_am*2
1474  IF (tot_num_exp_per_l(la) >= iii) THEN
1475  iii = tot_num_exp_per_l(la)
1476  jjj = la
1477  END IF
1478  nsgfa_ri = nsgfa_ri + ri_num_sgf_per_l(la)*(la*2 + 1)
1479  END DO
1480  DEALLOCATE (tot_num_exp_per_l)
1481  IF (real(nsgfa_ri, kind=dp)/real(nsgfa, kind=dp) <= 2.5_dp) THEN
1482  ri_num_sgf_per_l(jjj) = ri_num_sgf_per_l(jjj) + 1
1483  END IF
1484  END IF
1485 
1486  ri_nset = sum(ri_num_sgf_per_l)
1487 
1488  ALLOCATE (ri_exponents(ri_nset))
1489  ri_exponents = 0.0_dp
1490 
1491  ALLOCATE (ri_l_expo(ri_nset))
1492  ri_l_expo = 0
1493 
1494  ALLOCATE (max_rel_dev(ri_nset))
1495  max_rel_dev = 1.0_dp
1496 
1497  iset = 0
1498  DO la = 0, ri_max_am
1499  IF (ri_num_sgf_per_l(la) == 1) THEN
1500  iset = iset + 1
1501  ri_exponents(iset) = (max_min_exp_l(2, la) + max_min_exp_l(1, la))/2.0_dp
1502  ri_l_expo(iset) = la
1503  geom_fact = max_min_exp_l(2, la)/max_min_exp_l(1, la)
1504 
1505  max_rel_dev(iset) = (geom_fact - 1.0_dp)/(geom_fact + 1.0_dp)*0.9_dp
1506  ELSE
1507  geom_fact = (max_min_exp_l(2, la)/max_min_exp_l(1, la))**(1.0_dp/real(ri_num_sgf_per_l(la) - 1, dp))
1508  DO iexpo = 1, ri_num_sgf_per_l(la)
1509  iset = iset + 1
1510  ri_exponents(iset) = max_min_exp_l(1, la)*(geom_fact**(iexpo - 1))
1511  ri_l_expo(iset) = la
1512  max_rel_dev(iset) = (geom_fact - 1.0_dp)/(geom_fact + 1.0_dp)*0.9_dp
1513  END DO
1514  END IF
1515  END DO
1516  DEALLOCATE (num_sgf_per_l)
1517  DEALLOCATE (max_min_exp_l)
1518  DEALLOCATE (ri_num_sgf_per_l)
1519 
1520  END IF ! RI basis not associated
1521 
1522  ! create the new basis
1523  NULLIFY (tmp_basis)
1524  CALL create_ri_basis(tmp_basis, ri_nset, ri_l_expo, ri_exponents)
1525  CALL add_basis_set_to_container(qs_kind_set(ikind)%basis_sets, tmp_basis, "RI_AUX")
1526 
1527  DEALLOCATE (ri_exponents)
1528  DEALLOCATE (ri_l_expo)
1529 
1530  IF (.NOT. ALLOCATED(max_rel_dev_output)) THEN
1531  ALLOCATE (max_rel_dev_output(ri_nset))
1532  max_rel_dev_output(:) = max_rel_dev
1533  ELSE
1534  ! make a copy
1535  ri_prev_size = SIZE(max_rel_dev_output)
1536  ALLOCATE (max_rel_dev_prev(ri_prev_size))
1537  max_rel_dev_prev(:) = max_rel_dev_output
1538  DEALLOCATE (max_rel_dev_output)
1539  ! reallocate and copy
1540  ALLOCATE (max_rel_dev_output(ri_prev_size + ri_nset))
1541  max_rel_dev_output(1:ri_prev_size) = max_rel_dev_prev
1542  max_rel_dev_output(ri_prev_size + 1:ri_prev_size + ri_nset) = max_rel_dev
1543  DEALLOCATE (max_rel_dev_prev)
1544  END IF
1545  DEALLOCATE (max_rel_dev)
1546 
1547  END DO ! ikind
1548 
1549  IF (external_num_of_func) DEALLOCATE (mp2_env%ri_opt_param%RI_nset_per_l)
1550 
1551  CALL timestop(handle)
1552 
1553  END SUBROUTINE generate_ri_init_basis
1554 
1555 ! **************************************************************************************************
1556 !> \brief ...
1557 !> \param gto_basis_set ...
1558 !> \param RI_nset ...
1559 !> \param RI_l_expo ...
1560 !> \param RI_exponents ...
1561 ! **************************************************************************************************
1562  SUBROUTINE create_ri_basis(gto_basis_set, RI_nset, RI_l_expo, RI_exponents)
1563  TYPE(gto_basis_set_type), POINTER :: gto_basis_set
1564  INTEGER, INTENT(IN) :: ri_nset
1565  INTEGER, DIMENSION(:), INTENT(IN) :: ri_l_expo
1566  REAL(dp), DIMENSION(:), INTENT(IN) :: ri_exponents
1567 
1568  CHARACTER(LEN=*), PARAMETER :: routinen = 'create_ri_basis'
1569 
1570  INTEGER :: handle, i, ico, ipgf, iset, ishell, &
1571  lshell, m, maxco, maxl, maxpgf, &
1572  maxshell, ncgf, nmin, nset, nsgf
1573  INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
1574  INTEGER, DIMENSION(:, :), POINTER :: l, n
1575  REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
1576  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
1577 
1578  CALL timeset(routinen, handle)
1579  NULLIFY (lmax, lmin, npgf, nshell, l, n, zet, gcc)
1580 
1581  ! allocate the basis
1582  CALL allocate_gto_basis_set(gto_basis_set)
1583 
1584  ! brute force
1585  nset = 0
1586  maxshell = 0
1587  maxpgf = 0
1588  maxco = 0
1589  ncgf = 0
1590  nsgf = 0
1591  gto_basis_set%nset = nset
1592  gto_basis_set%ncgf = ncgf
1593  gto_basis_set%nsgf = nsgf
1594  CALL reallocate(gto_basis_set%lmax, 1, nset)
1595  CALL reallocate(gto_basis_set%lmin, 1, nset)
1596  CALL reallocate(gto_basis_set%npgf, 1, nset)
1597  CALL reallocate(gto_basis_set%nshell, 1, nset)
1598  CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
1599  CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
1600  CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
1601  CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
1602  CALL reallocate(gto_basis_set%set_radius, 1, nset)
1603  CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
1604  CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
1605  CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
1606  CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
1607  CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
1608  CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
1609  CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
1610  CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
1611  CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
1612  CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
1613  CALL reallocate(gto_basis_set%lx, 1, ncgf)
1614  CALL reallocate(gto_basis_set%ly, 1, ncgf)
1615  CALL reallocate(gto_basis_set%lz, 1, ncgf)
1616  CALL reallocate(gto_basis_set%m, 1, nsgf)
1617  CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
1618 
1619  nset = ri_nset
1620 
1621  ! locals
1622  CALL reallocate(npgf, 1, nset)
1623  CALL reallocate(nshell, 1, nset)
1624  CALL reallocate(lmax, 1, nset)
1625  CALL reallocate(lmin, 1, nset)
1626  CALL reallocate(n, 1, 1, 1, nset)
1627 
1628  maxl = 0
1629  maxpgf = 0
1630  maxshell = 0
1631  DO iset = 1, nset
1632  n(1, iset) = iset
1633  lmin(iset) = ri_l_expo(iset)
1634  lmax(iset) = ri_l_expo(iset)
1635  npgf(iset) = 1
1636 
1637  maxl = max(maxl, lmax(iset))
1638 
1639  IF (npgf(iset) > maxpgf) THEN
1640  maxpgf = npgf(iset)
1641  CALL reallocate(zet, 1, maxpgf, 1, nset)
1642  CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
1643  END IF
1644  nshell(iset) = 0
1645  DO lshell = lmin(iset), lmax(iset)
1646  nmin = n(1, iset) + lshell - lmin(iset)
1647  ishell = 1
1648  nshell(iset) = nshell(iset) + ishell
1649  IF (nshell(iset) > maxshell) THEN
1650  maxshell = nshell(iset)
1651  CALL reallocate(n, 1, maxshell, 1, nset)
1652  CALL reallocate(l, 1, maxshell, 1, nset)
1653  CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
1654  END IF
1655  DO i = 1, ishell
1656  n(nshell(iset) - ishell + i, iset) = nmin + i - 1
1657  l(nshell(iset) - ishell + i, iset) = lshell
1658  END DO
1659  END DO
1660 
1661  DO ipgf = 1, npgf(iset)
1662  zet(ipgf, iset) = ri_exponents(iset)
1663  DO ishell = 1, nshell(iset)
1664  gcc(ipgf, ishell, iset) = 1.0_dp
1665  END DO
1666  END DO
1667  END DO
1668 
1669  CALL init_orbital_pointers(maxl)
1670 
1671  gto_basis_set%nset = nset
1672  CALL reallocate(gto_basis_set%lmax, 1, nset)
1673  CALL reallocate(gto_basis_set%lmin, 1, nset)
1674  CALL reallocate(gto_basis_set%npgf, 1, nset)
1675  CALL reallocate(gto_basis_set%nshell, 1, nset)
1676  CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
1677  CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
1678  CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
1679  CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
1680 
1681  ! Copy the basis set information into the data structure
1682 
1683  DO iset = 1, nset
1684  gto_basis_set%lmax(iset) = lmax(iset)
1685  gto_basis_set%lmin(iset) = lmin(iset)
1686  gto_basis_set%npgf(iset) = npgf(iset)
1687  gto_basis_set%nshell(iset) = nshell(iset)
1688  DO ishell = 1, nshell(iset)
1689  gto_basis_set%n(ishell, iset) = n(ishell, iset)
1690  gto_basis_set%l(ishell, iset) = l(ishell, iset)
1691  DO ipgf = 1, npgf(iset)
1692  gto_basis_set%gcc(ipgf, ishell, iset) = gcc(ipgf, ishell, iset)
1693  END DO
1694  END DO
1695  DO ipgf = 1, npgf(iset)
1696  gto_basis_set%zet(ipgf, iset) = zet(ipgf, iset)
1697  END DO
1698  END DO
1699 
1700  ! Initialise the depending atomic kind information
1701 
1702  CALL reallocate(gto_basis_set%set_radius, 1, nset)
1703  CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
1704  CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
1705  CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
1706  CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
1707  CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
1708  CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
1709  CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
1710 
1711  maxco = 0
1712  ncgf = 0
1713  nsgf = 0
1714 
1715  DO iset = 1, nset
1716  gto_basis_set%ncgf_set(iset) = 0
1717  gto_basis_set%nsgf_set(iset) = 0
1718  DO ishell = 1, nshell(iset)
1719  lshell = gto_basis_set%l(ishell, iset)
1720  gto_basis_set%first_cgf(ishell, iset) = ncgf + 1
1721  ncgf = ncgf + nco(lshell)
1722  gto_basis_set%last_cgf(ishell, iset) = ncgf
1723  gto_basis_set%ncgf_set(iset) = &
1724  gto_basis_set%ncgf_set(iset) + nco(lshell)
1725  gto_basis_set%first_sgf(ishell, iset) = nsgf + 1
1726  nsgf = nsgf + nso(lshell)
1727  gto_basis_set%last_sgf(ishell, iset) = nsgf
1728  gto_basis_set%nsgf_set(iset) = &
1729  gto_basis_set%nsgf_set(iset) + nso(lshell)
1730  END DO
1731  maxco = max(maxco, npgf(iset)*ncoset(lmax(iset)))
1732  END DO
1733 
1734  gto_basis_set%ncgf = ncgf
1735  gto_basis_set%nsgf = nsgf
1736 
1737  CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
1738  CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
1739  CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
1740  CALL reallocate(gto_basis_set%lx, 1, ncgf)
1741  CALL reallocate(gto_basis_set%ly, 1, ncgf)
1742  CALL reallocate(gto_basis_set%lz, 1, ncgf)
1743  CALL reallocate(gto_basis_set%m, 1, nsgf)
1744  CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
1745  ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
1746 
1747  ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
1748 
1749  ncgf = 0
1750  nsgf = 0
1751 
1752  DO iset = 1, nset
1753  DO ishell = 1, nshell(iset)
1754  lshell = gto_basis_set%l(ishell, iset)
1755  DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
1756  ncgf = ncgf + 1
1757  gto_basis_set%lx(ncgf) = indco(1, ico)
1758  gto_basis_set%ly(ncgf) = indco(2, ico)
1759  gto_basis_set%lz(ncgf) = indco(3, ico)
1760  gto_basis_set%cgf_symbol(ncgf) = &
1761  cgf_symbol(n(ishell, iset), (/gto_basis_set%lx(ncgf), &
1762  gto_basis_set%ly(ncgf), &
1763  gto_basis_set%lz(ncgf)/))
1764  END DO
1765  DO m = -lshell, lshell
1766  nsgf = nsgf + 1
1767  gto_basis_set%m(nsgf) = m
1768  gto_basis_set%sgf_symbol(nsgf) = &
1769  sgf_symbol(n(ishell, iset), lshell, m)
1770  END DO
1771  END DO
1772  END DO
1773 
1774  DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet)
1775 
1776  CALL timestop(handle)
1777 
1778  END SUBROUTINE
1779 
1780 END MODULE mp2_optimize_ri_basis
1781 
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
basis_set_container_types
Definition: basis_set_container_types.F:13
basis_set_container_types::remove_basis_from_container
subroutine, public remove_basis_from_container(container, inum, basis_type)
...
Definition: basis_set_container_types.F:178
basis_set_container_types::add_basis_set_to_container
subroutine, public add_basis_set_to_container(container, basis_set, basis_set_type)
...
Definition: basis_set_container_types.F:149
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::allocate_gto_basis_set
subroutine, public allocate_gto_basis_set(gto_basis_set)
...
Definition: basis_set_types.F:150
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition: cp_log_handling.F:567
cp_log_handling::cp_logger_set
subroutine, public cp_logger_set(logger, local_filename, global_filename)
sets various attributes of the given logger
Definition: cp_log_handling.F:717
cp_log_handling::cp_rm_default_logger
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
Definition: cp_log_handling.F:213
cp_log_handling::cp_logger_release
subroutine, public cp_logger_release(logger)
releases this logger
Definition: cp_log_handling.F:423
cp_log_handling::cp_logger_create
subroutine, public cp_logger_create(logger, para_env, print_level, default_global_unit_nr, default_local_unit_nr, global_filename, local_filename, close_global_unit_on_dealloc, iter_info, close_local_unit_on_dealloc, suffix, template_logger)
initializes a logger
Definition: cp_log_handling.F:285
cp_log_handling::cp_add_default_logger
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
Definition: cp_log_handling.F:188
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition: libint_wrapper.F:15
libint_wrapper::build_eri_size
integer, parameter, public build_eri_size
Definition: libint_wrapper.F:66
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::default_output_unit
integer, parameter, public default_output_unit
Definition: machine.F:45
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition: message_passing.F:2027
mp2_direct_method
Routines to calculate MP2 energy.
Definition: mp2_direct_method.F:14
mp2_direct_method::mp2_canonical_direct_single_batch
subroutine, public mp2_canonical_direct_single_batch(Emp2, Emp2_Cou, Emp2_ex, mp2_env, qs_env, para_env, mp2_biel, dimen, C, Auto, i_batch_start, Ni_occupied, occupied, elements_ij_proc, ij_list_proc, Nj_occupied, j_batch_start, occupied_beta, C_beta, Auto_beta, Integ_MP2)
...
Definition: mp2_direct_method.F:386
mp2_optimize_ri_basis
Routines to optimize the RI-MP2 basis. Only exponents of non-contracted auxiliary basis basis are opt...
Definition: mp2_optimize_ri_basis.F:17
mp2_optimize_ri_basis::optimize_ri_basis_main
subroutine, public optimize_ri_basis_main(Emp2, Emp2_Cou, Emp2_ex, Emp2_S, Emp2_T, dimen, natom, homo, mp2_biel, mp2_env, C, Auto, kind_of, qs_env, para_env, unit_nr, homo_beta, C_beta, Auto_beta)
optimize RI-MP2 basis set
Definition: mp2_optimize_ri_basis.F:97
mp2_ri_libint
Routines to calculate the 3 and 2 center ERI's needed in the RI approximation using libint.
Definition: mp2_ri_libint.F:15
mp2_ri_libint::read_ri_basis_set
subroutine, public read_ri_basis_set(qs_env, RI_basis_parameter, RI_basis_info, natom, nkind, kind_of, RI_index_table, RI_dimen, basis_S0)
Read the auxiliary basis set for RI approxiamtion.
Definition: mp2_ri_libint.F:159
mp2_ri_libint::release_ri_basis_set
subroutine, public release_ri_basis_set(RI_basis_parameter, basis_S0)
Release the auxiliary basis set for RI approxiamtion (to be used only in the case of basis optimizati...
Definition: mp2_ri_libint.F:354
mp2_ri_libint::libint_ri_mp2
subroutine, public libint_ri_mp2(dimen, RI_dimen, occupied, natom, mp2_biel, mp2_env, C, kind_of, RI_basis_parameter, RI_basis_info, basis_S0, RI_index_table, qs_env, para_env, Lai)
...
Definition: mp2_ri_libint.F:104
mp2_types
Types needed for MP2 calculations.
Definition: mp2_types.F:14
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition: orbital_pointers.F:253
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition: orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition: orbital_pointers.F:43
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition: orbital_pointers.F:42
orbital_symbols
orbital_symbols
Definition: orbital_symbols.F:14
orbital_symbols::cgf_symbol
character(len=12) function, public cgf_symbol(n, lxyz)
Build a Cartesian orbital symbol (orbital labels for printing).
Definition: orbital_symbols.F:44
orbital_symbols::sgf_symbol
character(len=6) function, public sgf_symbol(n, l, m)
Build a spherical orbital symbol (orbital labels for printing).
Definition: orbital_symbols.F:107
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
util
All kind of helpful little routines.
Definition: util.F:14