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 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 6 : SUBROUTINE optimize_ri_basis_main(Emp2, Emp2_Cou, Emp2_ex, Emp2_S, Emp2_T, dimen, natom, homo, &
94 6 : mp2_biel, mp2_env, C, Auto, kind_of, &
95 : qs_env, para_env, &
96 6 : 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 6 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ij_list_proc, index_table_RI
119 : LOGICAL :: hess_up, open_shell_case, reset_boundary
120 6 : 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 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: deriv, dg, g, hdg, lower_B, max_dev, &
125 6 : max_rel_dev, p, pnew, xi
126 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: exp_limits, hessin
127 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: Integ_MP2, Integ_MP2_AA, Integ_MP2_AB, &
128 6 : Integ_MP2_BB
129 6 : 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 6 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_S0, RI_basis_parameter
134 : TYPE(mp_para_env_type), POINTER :: para_env_sub
135 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
136 :
137 6 : CALL timeset(routineN, handle)
138 6 : logger => cp_get_default_logger()
139 :
140 6 : open_shell_case = .FALSE.
141 6 : IF (PRESENT(homo_beta) .AND. PRESENT(C_beta) .AND. PRESENT(Auto_beta)) open_shell_case = .TRUE.
142 :
143 6 : virtual = dimen - homo
144 :
145 6 : eps_DRI = mp2_env%ri_opt_param%DRI
146 6 : eps_DI_rel = mp2_env%ri_opt_param%DI_rel
147 6 : eps_step = mp2_env%ri_opt_param%eps_step
148 6 : max_num_iter = mp2_env%ri_opt_param%max_num_iter
149 :
150 : ! calculate the ERI's over molecular integrals
151 6 : Emp2 = 0.0_dp
152 6 : Emp2_Cou = 0.0_dp
153 6 : Emp2_ex = 0.0_dp
154 6 : Emp2_S = 0.0_dp
155 6 : Emp2_T = 0.0_dp
156 6 : IF (open_shell_case) THEN
157 : ! open shell case
158 6 : virtual_beta = dimen - homo_beta
159 :
160 : ! alpha-aplha case
161 6 : Emp2_AA = 0.0_dp
162 6 : Emp2_AA_Cou = 0.0_dp
163 6 : Emp2_AA_ex = 0.0_dp
164 6 : 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 6 : Integ_MP2=Integ_MP2_AA)
169 6 : CALL para_env%sum(Emp2_AA_Cou)
170 6 : CALL para_env%sum(Emp2_AA_Ex)
171 6 : CALL para_env%sum(Emp2_AA)
172 6 : DEALLOCATE (ij_list_proc)
173 :
174 : ! beta-beta case
175 6 : Emp2_BB = 0.0_dp
176 6 : Emp2_BB_Cou = 0.0_dp
177 6 : Emp2_BB_ex = 0.0_dp
178 6 : 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 6 : Integ_MP2=Integ_MP2_BB)
183 6 : CALL para_env%sum(Emp2_BB_Cou)
184 6 : CALL para_env%sum(Emp2_BB_Ex)
185 6 : CALL para_env%sum(Emp2_BB)
186 6 : DEALLOCATE (ij_list_proc)
187 :
188 : ! aplha-beta case
189 6 : Emp2_AB = 0.0_dp
190 6 : Emp2_AB_Cou = 0.0_dp
191 6 : Emp2_AB_ex = 0.0_dp
192 6 : 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 6 : homo_beta, C_beta, Auto_beta, Integ_MP2=Integ_MP2_AB)
197 6 : CALL para_env%sum(Emp2_AB_Cou)
198 6 : CALL para_env%sum(Emp2_AB_Ex)
199 6 : CALL para_env%sum(Emp2_AB)
200 6 : DEALLOCATE (ij_list_proc)
201 :
202 6 : Emp2 = Emp2_AA + Emp2_BB + Emp2_AB*2.0_dp !+Emp2_BA
203 6 : Emp2_Cou = Emp2_AA_Cou + Emp2_BB_Cou + Emp2_AB_Cou*2.0_dp !+Emp2_BA
204 6 : Emp2_ex = Emp2_AA_ex + Emp2_BB_ex + Emp2_AB_ex*2.0_dp !+Emp2_BA
205 :
206 6 : Emp2_S = Emp2_AB*2.0_dp
207 6 : Emp2_T = Emp2_AA + Emp2_BB
208 :
209 : ! Replicate the MO-ERI's over all processes
210 6 : CALL para_env%sum(Integ_MP2_AA)
211 6 : CALL para_env%sum(Integ_MP2_BB)
212 6 : CALL para_env%sum(Integ_MP2_AB)
213 :
214 : ELSE
215 : ! close shell case
216 0 : 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 0 : Integ_MP2=Integ_MP2)
221 0 : CALL para_env%sum(Emp2_Cou)
222 0 : CALL para_env%sum(Emp2_Ex)
223 0 : CALL para_env%sum(Emp2)
224 0 : DEALLOCATE (ij_list_proc)
225 :
226 : ! Replicate the MO-ERI's over all processes
227 0 : CALL para_env%sum(Integ_MP2)
228 :
229 : END IF
230 :
231 : ! create the para_env_sub
232 6 : number_groups = para_env%num_pe/mp2_env%mp2_num_proc
233 6 : color_sub = para_env%mepos/mp2_env%mp2_num_proc
234 6 : ALLOCATE (para_env_sub)
235 6 : CALL para_env_sub%from_split(para_env, color_sub)
236 :
237 6 : IF (para_env%is_source()) THEN
238 3 : local_unit_nr = cp_logger_get_default_unit_nr(logger, local=.FALSE.)
239 : ELSE
240 3 : local_unit_nr = default_output_unit
241 : END IF
242 6 : NULLIFY (logger_sub)
243 : CALL cp_logger_create(logger_sub, para_env=para_env_sub, &
244 6 : default_global_unit_nr=local_unit_nr, close_global_unit_on_dealloc=.FALSE.)
245 6 : CALL cp_logger_set(logger_sub, local_filename="opt_RI_basis_localLog")
246 6 : CALL cp_add_default_logger(logger_sub)
247 :
248 6 : 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 6 : basis_S0)
253 :
254 6 : ndof = 0
255 6 : max_l_am = 0
256 12 : DO ikind = 1, nkind
257 98 : DO iset = 1, RI_basis_parameter(ikind)%nset
258 86 : ndof = ndof + 1
259 1330 : max_l_am = MAX(max_l_am, MAXVAL(RI_basis_parameter(ikind)%lmax))
260 : END DO
261 : END DO
262 :
263 18 : ALLOCATE (exp_limits(2, nkind))
264 12 : exp_limits(1, :) = HUGE(0)
265 12 : exp_limits(2, :) = -HUGE(0)
266 12 : DO ikind = 1, nkind
267 92 : DO iset = 1, RI_basis_parameter(ikind)%nset
268 86 : expon = RI_basis_parameter(ikind)%zet(1, iset)
269 86 : IF (expon <= exp_limits(1, ikind)) exp_limits(1, ikind) = expon
270 92 : IF (expon >= exp_limits(2, ikind)) exp_limits(2, ikind) = expon
271 : END DO
272 12 : IF (basis_was_assoc(ikind)) THEN
273 2 : exp_limits(1, ikind) = exp_limits(1, ikind)*0.5_dp
274 2 : exp_limits(2, ikind) = exp_limits(2, ikind)*1.5_dp
275 : ELSE
276 4 : exp_limits(1, ikind) = exp_limits(1, ikind)*0.8_dp*0.5_dp
277 4 : exp_limits(2, ikind) = exp_limits(2, ikind)*1.2_dp*1.5_dp
278 : END IF
279 : END DO
280 6 : DEALLOCATE (basis_was_assoc)
281 :
282 : ! check if the max angular momentum exceed the libint one
283 6 : IF (max_l_am > build_eri_size) THEN
284 0 : CPABORT("The angular momentum needed exceeds the value assumed when configuring libint.")
285 : END IF
286 :
287 : ! Allocate stuff
288 18 : ALLOCATE (p(ndof))
289 92 : p = 0.0_dp
290 18 : ALLOCATE (xi(ndof))
291 92 : xi = 0.0_dp
292 18 : ALLOCATE (g(ndof))
293 92 : g = 0.0_dp
294 18 : ALLOCATE (dg(ndof))
295 92 : dg = 0.0_dp
296 18 : ALLOCATE (hdg(ndof))
297 92 : hdg = 0.0_dp
298 18 : ALLOCATE (pnew(ndof))
299 92 : pnew = 0.0_dp
300 18 : ALLOCATE (deriv(ndof))
301 92 : deriv = 0.0_dp
302 24 : ALLOCATE (hessin(ndof, ndof))
303 1330 : hessin = 0.0_dp
304 92 : DO i = 1, ndof
305 92 : hessin(i, i) = 1.0_dp
306 : END DO
307 :
308 : ! initialize transformation function
309 18 : ALLOCATE (lower_B(ndof))
310 92 : lower_B = 0.0_dp
311 18 : ALLOCATE (max_dev(ndof))
312 92 : max_dev = 0.0_dp
313 :
314 : ! Initialize the transformation function
315 6 : 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 6 : 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 6 : .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 6 : deriv)
339 :
340 92 : g(:) = deriv
341 92 : xi(:) = -g
342 :
343 6 : reset_edge = 1.5_dp
344 180 : DO iiter = 1, max_num_iter
345 180 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,I5)') 'OPTIMIZATION STEP NUMBER', iiter
346 :
347 : ! perform step
348 2652 : pnew(:) = p + xi
349 180 : CALL p2basis(nkind, RI_basis_parameter, lower_B, max_dev, pnew)
350 :
351 : ! check if we have to reset boundaries
352 180 : reset_boundary = .FALSE.
353 180 : i = 0
354 360 : DO ikind = 1, nkind
355 2828 : DO iset = 1, RI_basis_parameter(ikind)%nset
356 2472 : i = i + 1
357 2472 : expon = transf_val(lower_B(i), max_dev(i), pnew(i))
358 2648 : 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 180 : IF (reset_boundary) THEN
367 4 : 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 4 : pnew, lower_B, max_dev, max_rel_dev, exp_limits)
370 60 : p(:) = pnew
371 60 : xi = 0.0_dp
372 60 : g = 0.0_dp
373 60 : dg = 0.0_dp
374 60 : hdg = 0.0_dp
375 60 : pnew = 0.0_dp
376 848 : hessin = 0.0_dp
377 60 : DO i = 1, ndof
378 60 : hessin(i, i) = 1.0_dp
379 : END DO
380 60 : 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 4 : .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 4 : deriv)
399 :
400 60 : g(:) = deriv
401 60 : xi(:) = -g
402 60 : pnew(:) = p + xi
403 4 : 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 180 : .FALSE.)
414 :
415 : ! update energy and direction
416 180 : DI = DI_new
417 2652 : xi(:) = pnew - p
418 2652 : p(:) = pnew
419 :
420 : ! check for convergence
421 180 : IF (unit_nr > 0) THEN
422 90 : WRITE (unit_nr, *)
423 180 : DO ikind = 1, nkind
424 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=ri_aux_basis, &
425 90 : basis_type="RI_AUX")
426 90 : WRITE (unit_nr, '(T3,A,A)') atomic_kind_set(ikind)%element_symbol, ' RI_opt_basis'
427 90 : WRITE (unit_nr, '(T3,I3)') RI_basis_parameter(ikind)%nset
428 1326 : DO iset = 1, RI_basis_parameter(ikind)%nset
429 1236 : WRITE (unit_nr, '(T3,10I4)') iset, &
430 1236 : RI_basis_parameter(ikind)%lmin(iset), &
431 1236 : RI_basis_parameter(ikind)%lmax(iset), &
432 1236 : RI_basis_parameter(ikind)%npgf(iset), &
433 3708 : (1, ishell=1, RI_basis_parameter(ikind)%nshell(iset))
434 2562 : DO ipgf = 1, RI_basis_parameter(ikind)%npgf(iset)
435 1236 : WRITE (unit_nr, '(T3,10F16.10)') RI_basis_parameter(ikind)%zet(ipgf, iset), &
436 2472 : (ri_aux_basis%gcc(ipgf, ishell, iset), &
437 6180 : ishell=1, ri_aux_basis%nshell(iset))
438 : END DO
439 : END DO
440 180 : WRITE (unit_nr, *)
441 : END DO
442 90 : WRITE (unit_nr, *)
443 90 : CALL m_flush(unit_nr)
444 : END IF
445 180 : IF (DI/ABS(Emp2) <= eps_DI_rel .AND. ABS(DRI_new) <= eps_DRI) THEN
446 6 : IF (unit_nr > 0) THEN
447 3 : WRITE (unit_nr, '(T3,A,/)') 'OPTIMIZATION CONVERGED'
448 3 : 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 174 : deriv)
463 :
464 : ! g is the vector containing the old gradient
465 2560 : dg(:) = deriv - g
466 2560 : g(:) = deriv
467 37824 : hdg(:) = MATMUL(hessin, dg)
468 :
469 2560 : fac = SUM(dg*xi)
470 2560 : fae = SUM(dg*hdg)
471 2560 : sumdg = SUM(dg*dg)
472 2560 : sumxi = SUM(xi*xi)
473 :
474 174 : hess_up = .TRUE.
475 174 : IF (fac**2 > sumdg*sumxi*3.0E-8_dp) THEN
476 174 : fac = 1.0_dp/fac
477 174 : fad = 1.0_dp/fae
478 2560 : dg(:) = fac*xi - fad*hdg
479 2560 : DO i = 1, ndof
480 35438 : DO j = 1, ndof
481 : hessin(i, j) = hessin(i, j) + fac*xi(i)*xi(j) &
482 : - fad*hdg(i)*hdg(j) &
483 35264 : + fae*dg(i)*dg(j)
484 : END DO
485 : END DO
486 : ELSE
487 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Skip Hessian Update'
488 174 : hess_up = .FALSE.
489 : END IF
490 :
491 : ! new direction
492 40558 : xi(:) = -MATMUL(hessin, g)
493 :
494 : END DO
495 :
496 6 : IF (.NOT. (DI/ABS(Emp2) <= eps_DI_rel .AND. ABS(DRI_new) <= eps_DRI)) THEN
497 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,I5,A)') 'OPTIMIZATION NOT CONVERGED IN', max_num_iter, ' STEPS.'
498 0 : IF (unit_nr > 0) WRITE (unit_nr, *)
499 : END IF
500 :
501 6 : DEALLOCATE (max_rel_dev)
502 :
503 6 : DEALLOCATE (p)
504 6 : DEALLOCATE (xi)
505 6 : DEALLOCATE (g)
506 6 : DEALLOCATE (pnew)
507 6 : DEALLOCATE (dg)
508 6 : DEALLOCATE (hdg)
509 6 : DEALLOCATE (deriv)
510 6 : DEALLOCATE (Hessin)
511 6 : DEALLOCATE (lower_B)
512 6 : DEALLOCATE (max_dev)
513 6 : DEALLOCATE (exp_limits)
514 :
515 6 : IF (open_shell_case) THEN
516 6 : DEALLOCATE (Integ_MP2_AA)
517 6 : DEALLOCATE (Integ_MP2_BB)
518 6 : DEALLOCATE (Integ_MP2_AB)
519 : ELSE
520 0 : DEALLOCATE (Integ_MP2)
521 : END IF
522 6 : DEALLOCATE (index_table_RI)
523 :
524 : ! Release RI basis set
525 6 : CALL release_RI_basis_set(RI_basis_parameter, basis_S0)
526 :
527 6 : CALL mp_para_env_release(para_env_sub)
528 6 : CALL cp_rm_default_logger()
529 6 : CALL cp_logger_release(logger_sub)
530 :
531 6 : CALL timestop(handle)
532 :
533 18 : 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 368 : 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 184 : kind_of, index_table_RI, mp2_biel, mp2_env, Auto, C, &
582 : RI_basis_parameter, RI_basis_info, basis_S0, &
583 184 : open_shell_case, homo_beta, virtual_beta, Auto_beta, C_beta, &
584 : para_env, para_env_sub, number_groups, color_sub, unit_nr, &
585 184 : p, lower_B, max_dev, &
586 184 : 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 184 : CALL timeset(routineN, handle)
626 :
627 184 : step = eps
628 :
629 : ! cycle over the RI basis set exponent
630 2712 : deriv = 0.0_dp
631 184 : ideriv = 0
632 368 : DO ikind = 1, nkind
633 184 : nseta = RI_basis_parameter(ikind)%nset
634 2896 : DO iset = 1, nseta
635 : ! for now only uncontracted aux basis set
636 2528 : ideriv = ideriv + 1
637 2528 : 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 1264 : CPASSERT(RI_basis_parameter(ikind)%npgf(iset) == 1)
646 1264 : orig_basis_val = RI_basis_parameter(ikind)%zet(1, iset)
647 1264 : temp = p(ideriv) + step
648 1264 : new_basis_val = transf_val(lower_B(ideriv), max_dev(ideriv), temp)
649 1264 : 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 1264 : para_env_sub, unit_nr, .TRUE.)
658 :
659 1264 : RI_basis_parameter(ikind)%zet(1, iset) = orig_basis_val
660 :
661 1448 : IF (para_env_sub%mepos == 0) THEN
662 1264 : temp = EXP(DI)
663 1264 : temp = temp/EXP(DI_ref)
664 1264 : deriv(ideriv) = LOG(temp)/step
665 : END IF
666 :
667 : END DO
668 : END DO
669 :
670 5240 : CALL para_env%sum(deriv)
671 :
672 184 : CALL timestop(handle)
673 :
674 184 : 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 1454 : 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 1454 : kind_of, index_table_RI, mp2_biel, mp2_env, Auto, C, &
717 : RI_basis_parameter, RI_basis_info, basis_S0, &
718 1454 : 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 1454 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Lai, Lai_beta
754 :
755 1454 : 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 1454 : qs_env, para_env, Lai)
760 1454 : 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 1454 : qs_env, para_env, Lai_beta)
764 : END IF
765 :
766 : ! Contract integrals into energy
767 1454 : 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 1454 : 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 1454 : 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 1454 : Lai, Lai_beta, para_env)
785 :
786 1454 : Emp2_RI = Emp2_RI_AA + Emp2_RI_BB + Emp2_RI_AB
787 1454 : DRI = DRI_AA + DRI_BB + DRI_AB
788 1454 : 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 0 : Lai, Lai, para_env)
794 : END IF
795 :
796 1454 : IF (.NOT. no_write) THEN
797 184 : IF (unit_nr > 0) THEN
798 : WRITE (unit_nr, '(/,(T3,A,T56,F25.14))') &
799 92 : 'Emp2 =', Emp2, &
800 184 : 'Emp2-RI =', Emp2_RI
801 : WRITE (unit_nr, '(T3,A,T56,ES25.10)') &
802 92 : 'DRI =', DRI, &
803 92 : 'DI =', DI, &
804 184 : 'DI/|Emp2| =', DI/ABS(Emp2)
805 92 : CALL m_flush(unit_nr)
806 : END IF
807 : END IF
808 :
809 1454 : DEALLOCATE (Lai)
810 1454 : IF (open_shell_case) DEALLOCATE (Lai_beta)
811 :
812 1454 : CALL timestop(handle)
813 :
814 1454 : 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 10 : 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 10 : ipos = 0
834 20 : DO ikind = 1, nkind
835 162 : DO iset = 1, RI_basis_parameter(ikind)%nset
836 142 : ipos = ipos + 1
837 142 : lower_B(ipos) = RI_basis_parameter(ikind)%zet(1, iset)*(1.0_dp - max_rel_dev(ipos))
838 152 : 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 10 : 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 184 : 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 184 : ipos = 0
863 368 : DO ikind = 1, nkind
864 2896 : DO iset = 1, RI_basis_parameter(ikind)%nset
865 2528 : ipos = ipos + 1
866 2528 : valout = transf_val(lower_B(ipos), max_dev(ipos), p(ipos))
867 2712 : RI_basis_parameter(ikind)%zet(1, iset) = valout
868 : END DO
869 : END DO
870 :
871 184 : 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 4 : SUBROUTINE reset_basis(nkind, ndof, RI_basis_parameter, reset_edge, &
886 4 : 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 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nf_per_l
897 8 : INTEGER, DIMENSION(ndof) :: change_expo
898 4 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_to_be_changed
899 : REAL(KIND=dp) :: expo, geom_fact, pmax
900 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: max_min_exp_per_l
901 12 : REAL(KIND=dp), DIMENSION(ndof) :: new_expo, old_expo, old_lower_B, &
902 4 : old_max_dev, old_max_rel_dev, old_pnew
903 :
904 : ! make a copy of the original parameters
905 :
906 60 : old_pnew = pnew
907 60 : old_lower_B = lower_B
908 60 : old_max_dev = max_dev
909 60 : old_max_rel_dev = max_rel_dev
910 60 : old_expo = 0.0_dp
911 4 : ipos = 0
912 8 : DO ikind = 1, nkind
913 64 : DO iset = 1, RI_basis_parameter(ikind)%nset
914 56 : ipos = ipos + 1
915 60 : old_expo(ipos) = RI_basis_parameter(ikind)%zet(1, iset)
916 : END DO
917 : END DO
918 :
919 60 : pnew = 0.0_dp
920 60 : lower_B = 0.0_dp
921 60 : max_dev = 0.0_dp
922 60 : max_rel_dev = 0.0_dp
923 :
924 60 : change_expo = 0
925 :
926 60 : new_expo = 0.0_dp
927 : ipos = 0
928 : ipos_p = 0
929 8 : DO ikind = 1, nkind
930 60 : am_max = MAXVAL(RI_basis_parameter(ikind)%lmax(:))
931 12 : ALLOCATE (nf_per_l(0:am_max))
932 20 : nf_per_l = 0
933 12 : ALLOCATE (max_min_exp_per_l(2, 0:am_max))
934 20 : max_min_exp_per_l(1, :) = HUGE(0)
935 20 : max_min_exp_per_l(2, :) = -HUGE(0)
936 :
937 60 : DO iset = 1, RI_basis_parameter(ikind)%nset
938 56 : la = RI_basis_parameter(ikind)%lmax(iset)
939 56 : expo = RI_basis_parameter(ikind)%zet(1, iset)
940 56 : nf_per_l(la) = nf_per_l(la) + 1
941 56 : IF (expo <= max_min_exp_per_l(1, la)) max_min_exp_per_l(1, la) = expo
942 60 : IF (expo >= max_min_exp_per_l(2, la)) max_min_exp_per_l(2, la) = expo
943 : END DO
944 :
945 4 : max_min_exp_per_l(1, la) = MAX(max_min_exp_per_l(1, la), exp_limits(1, ikind))
946 4 : 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 12 : ALLOCATE (has_to_be_changed(0:am_max))
955 20 : has_to_be_changed = .FALSE.
956 20 : DO la = 0, am_max
957 16 : pmax = -HUGE(0)
958 72 : DO iexpo = 1, nf_per_l(la)
959 56 : ipos_p = ipos_p + 1
960 56 : 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 56 : expo = transf_val(old_lower_B(ipos_p), old_max_dev(ipos_p), old_pnew(ipos_p))
963 72 : IF (expo < exp_limits(1, ikind) .OR. expo > exp_limits(2, ikind)) has_to_be_changed(la) = .TRUE.
964 : END DO
965 20 : IF (pmax > reset_edge) has_to_be_changed(la) = .TRUE.
966 : END DO
967 :
968 20 : DO la = 0, am_max
969 20 : IF (nf_per_l(la) == 1) THEN
970 4 : ipos = ipos + 1
971 4 : new_expo(ipos) = max_min_exp_per_l(1, la)
972 4 : IF (new_expo(ipos) >= exp_limits(1, ikind) .AND. new_expo(ipos) <= exp_limits(2, ikind)) THEN
973 4 : max_rel_dev(ipos) = (new_expo(ipos) - old_lower_B(ipos))/new_expo(ipos)
974 4 : IF (max_rel_dev(ipos) <= 0.1_dp) max_rel_dev(ipos) = 0.8_dp
975 : ELSE
976 0 : new_expo(ipos) = (exp_limits(2, ikind) + exp_limits(1, ikind))/2.0_dp
977 0 : max_rel_dev(ipos) = (new_expo(ipos) - exp_limits(1, ikind))/new_expo(ipos)
978 : END IF
979 4 : IF (has_to_be_changed(la)) change_expo(ipos) = 1
980 : ELSE
981 12 : IF (max_min_exp_per_l(2, la)/max_min_exp_per_l(1, la) < 2.0_dp) THEN
982 0 : max_min_exp_per_l(1, la) = max_min_exp_per_l(1, la)*0.5
983 0 : max_min_exp_per_l(2, la) = max_min_exp_per_l(2, la)*1.5
984 : END IF
985 12 : 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 64 : DO iexpo = 1, nf_per_l(la)
987 52 : ipos = ipos + 1
988 52 : new_expo(ipos) = max_min_exp_per_l(1, la)*(geom_fact**(iexpo - 1))
989 52 : max_rel_dev(ipos) = (geom_fact - 1.0_dp)/(geom_fact + 1.0_dp)*0.9_dp
990 64 : IF (has_to_be_changed(la)) change_expo(ipos) = 1
991 : END DO
992 : END IF
993 : END DO
994 :
995 4 : DEALLOCATE (has_to_be_changed)
996 :
997 4 : DEALLOCATE (nf_per_l)
998 8 : DEALLOCATE (max_min_exp_per_l)
999 : END DO
1000 :
1001 : ipos = 0
1002 8 : DO ikind = 1, nkind
1003 64 : DO iset = 1, RI_basis_parameter(ikind)%nset
1004 56 : ipos = ipos + 1
1005 60 : RI_basis_parameter(ikind)%zet(1, iset) = new_expo(ipos)
1006 : END DO
1007 : END DO
1008 :
1009 4 : CALL init_transf(nkind, RI_basis_parameter, lower_B, max_dev, max_rel_dev)
1010 :
1011 4 : ipos = 0
1012 8 : DO ikind = 1, nkind
1013 64 : DO iset = 1, RI_basis_parameter(ikind)%nset
1014 56 : ipos = ipos + 1
1015 60 : IF (change_expo(ipos) == 0) THEN
1016 : ! restore original
1017 52 : pnew(ipos) = old_pnew(ipos)
1018 52 : lower_B(ipos) = old_lower_B(ipos)
1019 52 : max_dev(ipos) = old_max_dev(ipos)
1020 52 : max_rel_dev(ipos) = old_max_rel_dev(ipos)
1021 52 : RI_basis_parameter(ikind)%zet(1, iset) = old_expo(ipos)
1022 : END IF
1023 : END DO
1024 : END DO
1025 :
1026 4 : 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 4362 : SUBROUTINE contract_integrals(DI, Emp2_RI, DRI, Emp2, homo, homo_beta, virtual, virtual_beta, &
1049 : fact, fact2, calc_ex, &
1050 4362 : MOenerg, MOenerg_beta, abij, &
1051 4362 : 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 4362 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_ab
1065 :
1066 17448 : ALLOCATE (mat_ab(virtual, virtual_beta))
1067 :
1068 4362 : DI = 0.0_dp
1069 4362 : Emp2_RI = 0.0_dp
1070 4362 : ij_counter = 0
1071 15994 : DO j = 1, homo_beta
1072 56706 : DO i = 1, homo
1073 40712 : ij_counter = ij_counter + 1
1074 40712 : IF (MOD(ij_counter, para_env%num_pe) /= para_env%mepos) CYCLE
1075 3908484 : mat_ab = 0.0_dp
1076 155386344 : mat_ab(:, :) = MATMUL(TRANSPOSE(Lai(:, :, i)), Lai_beta(:, :, j))
1077 424768 : DO b = 1, virtual_beta
1078 3911144 : DO a = 1, virtual
1079 3495348 : IF (calc_ex) THEN
1080 2419020 : t_iajb = fact*abij(a, b, i, j) - abij(b, a, i, j)
1081 2419020 : t_iajb_RI = fact*mat_ab(a, b) - mat_ab(b, a)
1082 : ELSE
1083 1076328 : t_iajb = fact*abij(a, b, i, j)
1084 1076328 : t_iajb_RI = fact*mat_ab(a, b)
1085 : END IF
1086 3495348 : t_iajb = t_iajb/(MOenerg(a + homo) + MOenerg_beta(b + homo_beta) - MOenerg(i) - MOenerg_beta(j))
1087 3495348 : t_iajb_RI = t_iajb_RI/(MOenerg(a + homo) + MOenerg_beta(b + homo_beta) - MOenerg(i) - MOenerg_beta(j))
1088 :
1089 3495348 : Emp2_RI = Emp2_RI - t_iajb_RI*mat_ab(a, b)*fact2
1090 :
1091 3870432 : DI = DI - t_iajb*mat_ab(a, b)*fact2
1092 :
1093 : END DO
1094 : END DO
1095 : END DO
1096 : END DO
1097 4362 : CALL para_env%sum(DI)
1098 4362 : CALL para_env%sum(Emp2_RI)
1099 :
1100 4362 : DRI = Emp2 - Emp2_RI
1101 4362 : DI = 2.0D+00*DI - Emp2 - Emp2_RI
1102 :
1103 4362 : DEALLOCATE (mat_ab)
1104 :
1105 4362 : 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 18 : 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 18 : elements_ij_proc = 0
1124 18 : ij_counter = -1
1125 78 : DO i = 1, homo
1126 246 : DO j = 1, homo_beta
1127 168 : ij_counter = ij_counter + 1
1128 228 : 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 54 : ALLOCATE (ij_list_proc(elements_ij_proc, 2))
1133 222 : ij_list_proc = 0
1134 18 : ij_counter = -1
1135 18 : elements_ij_proc = 0
1136 78 : DO i = 1, homo
1137 246 : DO j = 1, homo_beta
1138 168 : ij_counter = ij_counter + 1
1139 228 : IF (MOD(ij_counter, para_env%num_pe) == para_env%mepos) THEN
1140 84 : elements_ij_proc = elements_ij_proc + 1
1141 84 : ij_list_proc(elements_ij_proc, 1) = i
1142 84 : ij_list_proc(elements_ij_proc, 2) = j
1143 : END IF
1144 : END DO
1145 : END DO
1146 :
1147 18 : 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 6320 : 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 6320 : valout = 0.0_dp
1163 6320 : valout = lower_B + max_dev/(1.0_dp + EXP(-alpha*valin))
1164 :
1165 6320 : 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 6 : 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 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: l_expo, num_sgf_per_l, ordered_pos, &
1189 6 : RI_l_expo, RI_num_sgf_per_l, &
1190 6 : tot_num_exp_per_l
1191 6 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: l_tab
1192 6 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell
1193 6 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl
1194 : LOGICAL :: external_num_of_func
1195 : REAL(dp) :: geom_fact
1196 6 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: exponents, RI_exponents
1197 6 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: exp_tab, max_min_exp_l
1198 6 : REAL(dp), DIMENSION(:, :), POINTER :: sphi_a, zet
1199 6 : REAL(dp), DIMENSION(:, :, :), POINTER :: gcc
1200 6 : 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 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1203 : TYPE(qs_kind_type), POINTER :: atom_kind
1204 :
1205 6 : CALL timeset(routineN, handle)
1206 :
1207 6 : basis_quality = mp2_env%ri_opt_param%basis_quality
1208 6 : external_num_of_func = .FALSE.
1209 6 : IF (ALLOCATED(mp2_env%ri_opt_param%RI_nset_per_l)) external_num_of_func = .TRUE.
1210 :
1211 6 : NULLIFY (qs_kind_set)
1212 6 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
1213 6 : nkind = SIZE(qs_kind_set, 1)
1214 :
1215 18 : ALLOCATE (basis_was_assoc(nkind))
1216 12 : basis_was_assoc = .FALSE.
1217 :
1218 6 : 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 0 : "as specified by NUM_FUNC, will be used for all kinds.")
1222 : END IF
1223 :
1224 12 : DO ikind = 1, nkind
1225 6 : NULLIFY (atom_kind)
1226 6 : atom_kind => qs_kind_set(ikind)
1227 :
1228 6 : 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 6 : 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 6 : IF (basis_was_assoc(ikind)) THEN
1235 2 : nseta = orb_basis_a%nset
1236 2 : la_max => orb_basis_a%lmax
1237 2 : la_min => orb_basis_a%lmin
1238 2 : npgfa => orb_basis_a%npgf
1239 2 : nshell => orb_basis_a%nshell
1240 2 : zet => orb_basis_a%zet
1241 2 : nl => orb_basis_a%l
1242 :
1243 32 : max_am = MAXVAL(la_max)
1244 :
1245 2 : RI_max_am = max_am
1246 6 : ALLOCATE (RI_num_sgf_per_l(0:RI_max_am))
1247 10 : RI_num_sgf_per_l = 0
1248 2 : RI_nset = 0
1249 32 : DO iset = 1, nseta
1250 62 : DO la = la_min(iset), la_max(iset)
1251 30 : RI_nset = RI_nset + 1
1252 30 : RI_num_sgf_per_l(la) = RI_num_sgf_per_l(la) + 1
1253 60 : IF (npgfa(iset) > 1) THEN
1254 : CALL cp_warn(__LOCATION__, &
1255 : "The RI basis set optimizer can not handle contracted Gaussian. "// &
1256 0 : "Calculation continue with only uncontracted functions.")
1257 : END IF
1258 : END DO
1259 : END DO
1260 :
1261 16 : ALLOCATE (exp_tab(MAXVAL(RI_num_sgf_per_l), 0:RI_max_am))
1262 58 : exp_tab = 0.0_dp
1263 10 : DO iii = 0, RI_max_am
1264 : iexpo = 0
1265 130 : DO iset = 1, nseta
1266 248 : DO la = la_min(iset), la_max(iset)
1267 120 : IF (la /= iii) CYCLE
1268 30 : iexpo = iexpo + 1
1269 240 : exp_tab(iexpo, iii) = zet(1, iset)
1270 : END DO
1271 : END DO
1272 : END DO
1273 :
1274 : ! sort exponents
1275 10 : DO iii = 0, RI_max_am
1276 24 : ALLOCATE (ordered_pos(RI_num_sgf_per_l(iii)))
1277 38 : ordered_pos = 0
1278 8 : CALL sort(exp_tab(1:RI_num_sgf_per_l(iii), iii), RI_num_sgf_per_l(iii), ordered_pos)
1279 10 : DEALLOCATE (ordered_pos)
1280 : END DO
1281 :
1282 6 : ALLOCATE (RI_l_expo(RI_nset))
1283 32 : RI_l_expo = 0
1284 6 : ALLOCATE (RI_exponents(RI_nset))
1285 32 : RI_exponents = 0.0_dp
1286 :
1287 : iset = 0
1288 10 : DO iii = 0, RI_max_am
1289 40 : DO iexpo = 1, RI_num_sgf_per_l(iii)
1290 30 : iset = iset + 1
1291 30 : RI_l_expo(iset) = iii
1292 38 : RI_exponents(iset) = exp_tab(iexpo, iii)
1293 : END DO
1294 : END DO
1295 2 : DEALLOCATE (exp_tab)
1296 :
1297 6 : ALLOCATE (max_rel_dev(RI_nset))
1298 32 : max_rel_dev = 1.0_dp
1299 : iset = 0
1300 10 : DO iii = 0, RI_max_am
1301 10 : IF (RI_num_sgf_per_l(iii) == 1) THEN
1302 0 : iset = iset + 1
1303 0 : max_rel_dev(iset) = 0.35_dp
1304 : ELSE
1305 8 : iset = iset + 1
1306 8 : max_rel_dev(iset) = (RI_exponents(iset + 1) + RI_exponents(iset))/2.0_dp
1307 8 : max_rel_dev(iset) = max_rel_dev(iset)/RI_exponents(iset) - 1.0_dp
1308 30 : DO iexpo = 2, RI_num_sgf_per_l(iii)
1309 22 : iset = iset + 1
1310 22 : max_rel_dev(iset) = (RI_exponents(iset) + RI_exponents(iset - 1))/2.0_dp
1311 30 : max_rel_dev(iset) = 1.0_dp - max_rel_dev(iset)/RI_exponents(iset)
1312 : END DO
1313 : END IF
1314 : END DO
1315 32 : max_rel_dev(:) = max_rel_dev(:)*0.9_dp
1316 2 : DEALLOCATE (RI_num_sgf_per_l)
1317 :
1318 : ! deallocate the old basis before moving on
1319 2 : CALL remove_basis_from_container(qs_kind_set(ikind)%basis_sets, basis_type="RI_AUX")
1320 :
1321 : ELSE
1322 :
1323 4 : CALL get_qs_kind(qs_kind=atom_kind, basis_set=orb_basis_a)
1324 :
1325 4 : sphi_a => orb_basis_a%sphi
1326 4 : nseta = orb_basis_a%nset
1327 4 : la_max => orb_basis_a%lmax
1328 4 : la_min => orb_basis_a%lmin
1329 4 : npgfa => orb_basis_a%npgf
1330 4 : first_sgfa => orb_basis_a%first_sgf
1331 4 : nshell => orb_basis_a%nshell
1332 4 : zet => orb_basis_a%zet
1333 4 : gcc => orb_basis_a%gcc
1334 4 : nl => orb_basis_a%l
1335 4 : nsgfa = orb_basis_a%nsgf
1336 :
1337 16 : max_am = MAXVAL(la_max)
1338 :
1339 : nexpo_shell = 0
1340 16 : DO iset = 1, nseta
1341 36 : DO ishell = 1, nshell(iset)
1342 68 : DO la = la_min(iset), la_max(iset)
1343 36 : IF (ishell > 1) THEN
1344 16 : IF (nl(ishell, iset) == nl(ishell - 1, iset)) CYCLE
1345 : END IF
1346 36 : IF (la /= nl(ishell, iset)) CYCLE
1347 56 : nexpo_shell = nexpo_shell + npgfa(iset)
1348 : END DO
1349 : END DO
1350 : END DO
1351 :
1352 12 : ALLOCATE (exponents(nexpo_shell))
1353 40 : exponents = 0.0_dp
1354 12 : ALLOCATE (l_expo(nexpo_shell))
1355 40 : l_expo = 0
1356 12 : ALLOCATE (num_sgf_per_l(0:max_am))
1357 16 : num_sgf_per_l = 0
1358 : iexpo = 0
1359 16 : DO iset = 1, nseta
1360 36 : DO ishell = 1, nshell(iset)
1361 68 : DO la = la_min(iset), la_max(iset)
1362 36 : IF (ishell > 1) THEN
1363 16 : IF (nl(ishell, iset) == nl(ishell - 1, iset)) CYCLE
1364 : END IF
1365 36 : IF (la /= nl(ishell, iset)) CYCLE
1366 56 : DO ipgf = 1, npgfa(iset)
1367 36 : iexpo = iexpo + 1
1368 36 : exponents(iexpo) = zet(ipgf, iset)
1369 56 : l_expo(iexpo) = la
1370 : END DO
1371 56 : num_sgf_per_l(la) = num_sgf_per_l(la) + 1
1372 : END DO
1373 : END DO
1374 : END DO
1375 :
1376 16 : ALLOCATE (exp_tab(nexpo_shell, nexpo_shell))
1377 364 : exp_tab = 0.0_dp
1378 16 : ALLOCATE (l_tab(nexpo_shell, nexpo_shell))
1379 364 : l_tab = 0
1380 12 : ALLOCATE (tot_num_exp_per_l(0:max_am*2))
1381 24 : tot_num_exp_per_l = 0
1382 40 : DO iexpo = 1, nexpo_shell
1383 220 : DO jexpo = iexpo, nexpo_shell
1384 180 : exp_tab(jexpo, iexpo) = exponents(jexpo) + exponents(iexpo)
1385 180 : exp_tab(iexpo, jexpo) = exp_tab(jexpo, iexpo)
1386 180 : l_tab(jexpo, iexpo) = l_expo(jexpo) + l_expo(iexpo)
1387 180 : l_tab(iexpo, jexpo) = l_tab(jexpo, iexpo)
1388 216 : 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 4 : DEALLOCATE (l_expo)
1392 4 : DEALLOCATE (exponents)
1393 :
1394 12 : ALLOCATE (max_min_exp_l(2, 0:max_am*2))
1395 24 : max_min_exp_l(1, :) = HUGE(0)
1396 24 : max_min_exp_l(2, :) = -HUGE(0)
1397 :
1398 24 : DO la = 0, max_am*2
1399 204 : DO iexpo = 1, nexpo_shell
1400 1100 : DO jexpo = iexpo, nexpo_shell
1401 900 : IF (la /= l_tab(jexpo, iexpo)) CYCLE
1402 180 : IF (exp_tab(jexpo, iexpo) <= max_min_exp_l(1, la)) max_min_exp_l(1, la) = exp_tab(jexpo, iexpo)
1403 360 : 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 4 : DEALLOCATE (exp_tab)
1408 4 : 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 24 : max_min_exp_l(1, :) = max_min_exp_l(1, :)/0.80_dp
1413 24 : max_min_exp_l(2, :) = max_min_exp_l(2, :)/1.20_dp
1414 :
1415 12 : ALLOCATE (RI_num_sgf_per_l(0:max_am*2))
1416 24 : RI_num_sgf_per_l = 0
1417 :
1418 : SELECT CASE (basis_quality)
1419 : CASE (0)
1420 : ! normal
1421 : prog_func = 0
1422 2 : prog_l = 2
1423 : CASE (1)
1424 : ! large
1425 2 : prog_func = 1
1426 2 : prog_l = 2
1427 : CASE (2)
1428 0 : prog_func = 2
1429 0 : prog_l = 2
1430 : CASE DEFAULT
1431 : prog_func = 0
1432 4 : prog_l = 2
1433 : END SELECT
1434 :
1435 4 : IF (external_num_of_func) THEN
1436 : ! cp2k can not exceed angular momentum 7
1437 2 : RI_max_am = MIN(SIZE(mp2_env%ri_opt_param%RI_nset_per_l) - 1, 7)
1438 2 : IF (RI_max_am > max_am*2) THEN
1439 0 : DEALLOCATE (RI_num_sgf_per_l)
1440 0 : ALLOCATE (RI_num_sgf_per_l(0:RI_max_am))
1441 0 : RI_num_sgf_per_l = 0
1442 : END IF
1443 10 : DO la = 0, RI_max_am
1444 10 : RI_num_sgf_per_l(la) = mp2_env%ri_opt_param%RI_nset_per_l(la)
1445 : END DO
1446 : ELSE
1447 2 : RI_num_sgf_per_l(0) = num_sgf_per_l(0)*2 + prog_func
1448 6 : DO la = 1, max_am
1449 4 : RI_num_sgf_per_l(la) = RI_num_sgf_per_l(la - 1) - 1
1450 6 : IF (RI_num_sgf_per_l(la) == 0) THEN
1451 0 : RI_num_sgf_per_l(la) = 1
1452 0 : EXIT
1453 : END IF
1454 : END DO
1455 2 : DO la = max_am + 1, max_am*2
1456 2 : RI_num_sgf_per_l(la) = RI_num_sgf_per_l(la - 1) - prog_l
1457 2 : IF (RI_num_sgf_per_l(la) == 0) THEN
1458 0 : RI_num_sgf_per_l(la) = 1
1459 : END IF
1460 2 : IF (RI_num_sgf_per_l(la) == 1) EXIT
1461 : END DO
1462 2 : RI_max_am = MIN(max_am*2, 7)
1463 10 : DO la = 0, MIN(max_am*2, 7)
1464 10 : IF (RI_num_sgf_per_l(la) == 0) THEN
1465 2 : RI_max_am = la - 1
1466 2 : EXIT
1467 : END IF
1468 : END DO
1469 :
1470 2 : iii = 0
1471 2 : jjj = 0
1472 2 : nsgfa_RI = 0
1473 10 : DO la = 1, max_am*2
1474 8 : IF (tot_num_exp_per_l(la) >= iii) THEN
1475 2 : iii = tot_num_exp_per_l(la)
1476 2 : jjj = la
1477 : END IF
1478 10 : nsgfa_RI = nsgfa_RI + RI_num_sgf_per_l(la)*(la*2 + 1)
1479 : END DO
1480 2 : DEALLOCATE (tot_num_exp_per_l)
1481 2 : IF (REAL(nsgfa_RI, KIND=dp)/REAL(nsgfa, KIND=dp) <= 2.5_dp) THEN
1482 0 : RI_num_sgf_per_l(jjj) = RI_num_sgf_per_l(jjj) + 1
1483 : END IF
1484 : END IF
1485 :
1486 24 : RI_nset = SUM(RI_num_sgf_per_l)
1487 :
1488 12 : ALLOCATE (RI_exponents(RI_nset))
1489 60 : RI_exponents = 0.0_dp
1490 :
1491 12 : ALLOCATE (RI_l_expo(RI_nset))
1492 60 : RI_l_expo = 0
1493 :
1494 12 : ALLOCATE (max_rel_dev(RI_nset))
1495 60 : max_rel_dev = 1.0_dp
1496 :
1497 : iset = 0
1498 20 : DO la = 0, RI_max_am
1499 20 : IF (RI_num_sgf_per_l(la) == 1) THEN
1500 4 : iset = iset + 1
1501 4 : RI_exponents(iset) = (max_min_exp_l(2, la) + max_min_exp_l(1, la))/2.0_dp
1502 4 : RI_l_expo(iset) = la
1503 4 : geom_fact = max_min_exp_l(2, la)/max_min_exp_l(1, la)
1504 :
1505 4 : max_rel_dev(iset) = (geom_fact - 1.0_dp)/(geom_fact + 1.0_dp)*0.9_dp
1506 : ELSE
1507 12 : 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 64 : DO iexpo = 1, RI_num_sgf_per_l(la)
1509 52 : iset = iset + 1
1510 52 : RI_exponents(iset) = max_min_exp_l(1, la)*(geom_fact**(iexpo - 1))
1511 52 : RI_l_expo(iset) = la
1512 64 : 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 4 : DEALLOCATE (num_sgf_per_l)
1517 4 : DEALLOCATE (max_min_exp_l)
1518 4 : DEALLOCATE (RI_num_sgf_per_l)
1519 :
1520 : END IF ! RI basis not associated
1521 :
1522 : ! create the new basis
1523 6 : NULLIFY (tmp_basis)
1524 6 : CALL create_ri_basis(tmp_basis, RI_nset, RI_l_expo, RI_exponents)
1525 6 : CALL add_basis_set_to_container(qs_kind_set(ikind)%basis_sets, tmp_basis, "RI_AUX")
1526 :
1527 6 : DEALLOCATE (RI_exponents)
1528 6 : DEALLOCATE (RI_l_expo)
1529 :
1530 6 : IF (.NOT. ALLOCATED(max_rel_dev_output)) THEN
1531 18 : ALLOCATE (max_rel_dev_output(RI_nset))
1532 92 : max_rel_dev_output(:) = max_rel_dev
1533 : ELSE
1534 : ! make a copy
1535 0 : RI_prev_size = SIZE(max_rel_dev_output)
1536 0 : ALLOCATE (max_rel_dev_prev(RI_prev_size))
1537 0 : max_rel_dev_prev(:) = max_rel_dev_output
1538 0 : DEALLOCATE (max_rel_dev_output)
1539 : ! reallocate and copy
1540 0 : ALLOCATE (max_rel_dev_output(RI_prev_size + RI_nset))
1541 0 : max_rel_dev_output(1:RI_prev_size) = max_rel_dev_prev
1542 0 : max_rel_dev_output(RI_prev_size + 1:RI_prev_size + RI_nset) = max_rel_dev
1543 0 : DEALLOCATE (max_rel_dev_prev)
1544 : END IF
1545 12 : DEALLOCATE (max_rel_dev)
1546 :
1547 : END DO ! ikind
1548 :
1549 6 : IF (external_num_of_func) DEALLOCATE (mp2_env%ri_opt_param%RI_nset_per_l)
1550 :
1551 6 : CALL timestop(handle)
1552 :
1553 12 : 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 6 : 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 6 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
1574 6 : INTEGER, DIMENSION(:, :), POINTER :: l, n
1575 6 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
1576 6 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
1577 :
1578 6 : CALL timeset(routineN, handle)
1579 6 : NULLIFY (lmax, lmin, npgf, nshell, l, n, zet, gcc)
1580 :
1581 : ! allocate the basis
1582 6 : CALL allocate_gto_basis_set(gto_basis_set)
1583 :
1584 : ! brute force
1585 6 : nset = 0
1586 6 : maxshell = 0
1587 6 : maxpgf = 0
1588 6 : maxco = 0
1589 6 : ncgf = 0
1590 6 : nsgf = 0
1591 6 : gto_basis_set%nset = nset
1592 6 : gto_basis_set%ncgf = ncgf
1593 6 : gto_basis_set%nsgf = nsgf
1594 6 : CALL reallocate(gto_basis_set%lmax, 1, nset)
1595 6 : CALL reallocate(gto_basis_set%lmin, 1, nset)
1596 6 : CALL reallocate(gto_basis_set%npgf, 1, nset)
1597 6 : CALL reallocate(gto_basis_set%nshell, 1, nset)
1598 6 : CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
1599 6 : CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
1600 6 : CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
1601 6 : CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
1602 6 : CALL reallocate(gto_basis_set%set_radius, 1, nset)
1603 6 : CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
1604 6 : CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
1605 6 : CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
1606 6 : CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
1607 6 : CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
1608 6 : CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
1609 6 : CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
1610 6 : CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
1611 6 : CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
1612 6 : CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
1613 6 : CALL reallocate(gto_basis_set%lx, 1, ncgf)
1614 6 : CALL reallocate(gto_basis_set%ly, 1, ncgf)
1615 6 : CALL reallocate(gto_basis_set%lz, 1, ncgf)
1616 6 : CALL reallocate(gto_basis_set%m, 1, nsgf)
1617 6 : CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
1618 :
1619 6 : nset = RI_nset
1620 :
1621 : ! locals
1622 6 : CALL reallocate(npgf, 1, nset)
1623 6 : CALL reallocate(nshell, 1, nset)
1624 6 : CALL reallocate(lmax, 1, nset)
1625 6 : CALL reallocate(lmin, 1, nset)
1626 6 : CALL reallocate(n, 1, 1, 1, nset)
1627 :
1628 6 : maxl = 0
1629 : maxpgf = 0
1630 : maxshell = 0
1631 92 : DO iset = 1, nset
1632 86 : n(1, iset) = iset
1633 86 : lmin(iset) = RI_l_expo(iset)
1634 86 : lmax(iset) = RI_l_expo(iset)
1635 86 : npgf(iset) = 1
1636 :
1637 86 : maxl = MAX(maxl, lmax(iset))
1638 :
1639 86 : IF (npgf(iset) > maxpgf) THEN
1640 6 : maxpgf = npgf(iset)
1641 6 : CALL reallocate(zet, 1, maxpgf, 1, nset)
1642 6 : CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
1643 : END IF
1644 86 : nshell(iset) = 0
1645 172 : DO lshell = lmin(iset), lmax(iset)
1646 86 : nmin = n(1, iset) + lshell - lmin(iset)
1647 86 : ishell = 1
1648 86 : nshell(iset) = nshell(iset) + ishell
1649 86 : IF (nshell(iset) > maxshell) THEN
1650 6 : maxshell = nshell(iset)
1651 6 : CALL reallocate(n, 1, maxshell, 1, nset)
1652 6 : CALL reallocate(l, 1, maxshell, 1, nset)
1653 6 : CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
1654 : END IF
1655 258 : DO i = 1, ishell
1656 86 : n(nshell(iset) - ishell + i, iset) = nmin + i - 1
1657 172 : l(nshell(iset) - ishell + i, iset) = lshell
1658 : END DO
1659 : END DO
1660 :
1661 178 : DO ipgf = 1, npgf(iset)
1662 86 : zet(ipgf, iset) = RI_exponents(iset)
1663 258 : DO ishell = 1, nshell(iset)
1664 172 : gcc(ipgf, ishell, iset) = 1.0_dp
1665 : END DO
1666 : END DO
1667 : END DO
1668 :
1669 6 : CALL init_orbital_pointers(maxl)
1670 :
1671 6 : gto_basis_set%nset = nset
1672 6 : CALL reallocate(gto_basis_set%lmax, 1, nset)
1673 6 : CALL reallocate(gto_basis_set%lmin, 1, nset)
1674 6 : CALL reallocate(gto_basis_set%npgf, 1, nset)
1675 6 : CALL reallocate(gto_basis_set%nshell, 1, nset)
1676 6 : CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
1677 6 : CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
1678 6 : CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
1679 6 : 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 92 : DO iset = 1, nset
1684 86 : gto_basis_set%lmax(iset) = lmax(iset)
1685 86 : gto_basis_set%lmin(iset) = lmin(iset)
1686 86 : gto_basis_set%npgf(iset) = npgf(iset)
1687 86 : gto_basis_set%nshell(iset) = nshell(iset)
1688 172 : DO ishell = 1, nshell(iset)
1689 86 : gto_basis_set%n(ishell, iset) = n(ishell, iset)
1690 86 : gto_basis_set%l(ishell, iset) = l(ishell, iset)
1691 258 : DO ipgf = 1, npgf(iset)
1692 172 : gto_basis_set%gcc(ipgf, ishell, iset) = gcc(ipgf, ishell, iset)
1693 : END DO
1694 : END DO
1695 178 : DO ipgf = 1, npgf(iset)
1696 172 : 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 6 : CALL reallocate(gto_basis_set%set_radius, 1, nset)
1703 6 : CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
1704 6 : CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
1705 6 : CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
1706 6 : CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
1707 6 : CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
1708 6 : CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
1709 6 : CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
1710 :
1711 : maxco = 0
1712 : ncgf = 0
1713 : nsgf = 0
1714 :
1715 92 : DO iset = 1, nset
1716 86 : gto_basis_set%ncgf_set(iset) = 0
1717 86 : gto_basis_set%nsgf_set(iset) = 0
1718 172 : DO ishell = 1, nshell(iset)
1719 86 : lshell = gto_basis_set%l(ishell, iset)
1720 86 : gto_basis_set%first_cgf(ishell, iset) = ncgf + 1
1721 86 : ncgf = ncgf + nco(lshell)
1722 86 : gto_basis_set%last_cgf(ishell, iset) = ncgf
1723 : gto_basis_set%ncgf_set(iset) = &
1724 86 : gto_basis_set%ncgf_set(iset) + nco(lshell)
1725 86 : gto_basis_set%first_sgf(ishell, iset) = nsgf + 1
1726 86 : nsgf = nsgf + nso(lshell)
1727 86 : gto_basis_set%last_sgf(ishell, iset) = nsgf
1728 : gto_basis_set%nsgf_set(iset) = &
1729 172 : gto_basis_set%nsgf_set(iset) + nso(lshell)
1730 : END DO
1731 92 : maxco = MAX(maxco, npgf(iset)*ncoset(lmax(iset)))
1732 : END DO
1733 :
1734 6 : gto_basis_set%ncgf = ncgf
1735 6 : gto_basis_set%nsgf = nsgf
1736 :
1737 6 : CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
1738 6 : CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
1739 6 : CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
1740 6 : CALL reallocate(gto_basis_set%lx, 1, ncgf)
1741 6 : CALL reallocate(gto_basis_set%ly, 1, ncgf)
1742 6 : CALL reallocate(gto_basis_set%lz, 1, ncgf)
1743 6 : CALL reallocate(gto_basis_set%m, 1, nsgf)
1744 6 : CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
1745 18 : ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
1746 :
1747 18 : ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
1748 :
1749 6 : ncgf = 0
1750 6 : nsgf = 0
1751 :
1752 92 : DO iset = 1, nset
1753 178 : DO ishell = 1, nshell(iset)
1754 86 : lshell = gto_basis_set%l(ishell, iset)
1755 396 : DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
1756 310 : ncgf = ncgf + 1
1757 310 : gto_basis_set%lx(ncgf) = indco(1, ico)
1758 310 : gto_basis_set%ly(ncgf) = indco(2, ico)
1759 310 : 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 1326 : gto_basis_set%lz(ncgf)/))
1764 : END DO
1765 438 : DO m = -lshell, lshell
1766 266 : nsgf = nsgf + 1
1767 266 : gto_basis_set%m(nsgf) = m
1768 : gto_basis_set%sgf_symbol(nsgf) = &
1769 352 : sgf_symbol(n(ishell, iset), lshell, m)
1770 : END DO
1771 : END DO
1772 : END DO
1773 :
1774 6 : DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet)
1775 :
1776 6 : CALL timestop(handle)
1777 :
1778 6 : END SUBROUTINE
1779 :
1780 174 : END MODULE mp2_optimize_ri_basis
1781 :
|
