(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
33 USE kinds, ONLY: dp
35 USE machine, ONLY: default_output_unit,&
41 USE mp2_ri_libint, ONLY: libint_ri_mp2,&
44 USE mp2_types, ONLY: mp2_biel_type,&
46 USE orbital_pointers, ONLY: indco,&
48 nco,&
49 ncoset,&
50 nso
51 USE orbital_symbols, ONLY: cgf_symbol,&
55 USE qs_kind_types, ONLY: get_qs_kind,&
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
67
68CONTAINS
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)
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
1780END MODULE mp2_optimize_ri_basis
1781
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:48
Define the atomic kind types and their sub types.
subroutine, public remove_basis_from_container(container, inum, basis_type)
...
subroutine, public add_basis_set_to_container(container, basis_set, basis_set_type)
...
subroutine, public allocate_gto_basis_set(gto_basis_set)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
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
subroutine, public cp_logger_set(logger, local_filename, global_filename)
sets various attributes of the given 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 ...
subroutine, public cp_logger_release(logger)
releases this logger
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
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Types and set/get functions for HFX.
Definition hfx_types.F:15
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the Libint-Library or a c++ wrapper.
integer, parameter, public build_eri_size
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
integer, parameter, public default_output_unit
Definition machine.F:45
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:106
Utility routines for the memory handling.
Interface to the message passing library MPI.
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)
Routines to calculate MP2 energy.
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)
...
Routines to optimize the RI-MP2 basis. Only exponents of non-contracted auxiliary basis basis are opt...
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
Routines to calculate the 3 and 2 center ERI's needed in the RI approximation using libint.
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)
...
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.
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...
Types needed for MP2 calculations.
Definition mp2_types.F:14
Provides Cartesian and spherical orbital pointers and indices.
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
integer, dimension(:), allocatable, public nco
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :), allocatable, public indco
integer, dimension(:), allocatable, public nso
orbital_symbols
character(len=12) function, public cgf_symbol(n, lxyz)
Build a Cartesian orbital symbol (orbital labels for printing).
character(len=6) function, public sgf_symbol(n, l, m)
Build a spherical orbital symbol (orbital labels for printing).
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.
Define the quickstep kind type and their sub types.
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.
All kind of helpful little routines.
Definition util.F:14
Provides all information about an atomic kind.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.