58#include "./base/base_uses.f90"
64 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'mp2_optimize_ri_basis'
94 mp2_biel, mp2_env, C, Auto, kind_of, &
96 unit_nr, homo_beta, C_beta, Auto_beta)
98 REAL(kind=
dp),
INTENT(OUT) :: emp2, emp2_cou, emp2_ex, emp2_s, emp2_t
99 INTEGER,
INTENT(IN) :: dimen, natom, homo
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
107 INTEGER,
INTENT(IN) :: unit_nr
108 INTEGER,
INTENT(IN),
OPTIONAL :: homo_beta
109 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
111 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: auto_beta
113 CHARACTER(len=*),
PARAMETER :: routinen =
'optimize_ri_basis_main'
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, &
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, &
133 TYPE(
hfx_basis_type),
DIMENSION(:),
POINTER :: basis_s0, ri_basis_parameter
135 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
137 CALL timeset(routinen, handle)
140 open_shell_case = .false.
141 IF (
PRESENT(homo_beta) .AND.
PRESENT(c_beta) .AND.
PRESENT(auto_beta)) open_shell_case = .true.
143 virtual = dimen - homo
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
156 IF (open_shell_case)
THEN
158 virtual_beta = dimen - homo_beta
164 CALL calc_elem_ij_proc(homo, homo, para_env, elements_ij_proc, ij_list_proc)
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)
178 CALL calc_elem_ij_proc(homo_beta, homo_beta, para_env, elements_ij_proc, ij_list_proc)
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)
192 CALL calc_elem_ij_proc(homo, homo_beta, para_env, elements_ij_proc, ij_list_proc)
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)
202 emp2 = emp2_aa + emp2_bb + emp2_ab*2.0_dp
203 emp2_cou = emp2_aa_cou + emp2_bb_cou + emp2_ab_cou*2.0_dp
204 emp2_ex = emp2_aa_ex + emp2_bb_ex + emp2_ab_ex*2.0_dp
206 emp2_s = emp2_ab*2.0_dp
207 emp2_t = emp2_aa + emp2_bb
210 CALL para_env%sum(integ_mp2_aa)
211 CALL para_env%sum(integ_mp2_bb)
212 CALL para_env%sum(integ_mp2_ab)
216 CALL calc_elem_ij_proc(homo, homo, para_env, elements_ij_proc, ij_list_proc)
218 mp2_biel, dimen, c, auto, 0, homo, homo, &
219 elements_ij_proc, ij_list_proc, homo, 0, &
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)
227 CALL para_env%sum(integ_mp2)
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)
237 IF (para_env%is_source())
THEN
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")
248 CALL generate_ri_init_basis(qs_env, mp2_env, nkind, max_rel_dev, basis_was_assoc)
251 natom, nkind, kind_of, index_table_ri, dimen_ri, &
257 DO iset = 1, ri_basis_parameter(ikind)%nset
259 max_l_am = max(max_l_am, maxval(ri_basis_parameter(ikind)%lmax))
263 ALLOCATE (exp_limits(2, nkind))
264 exp_limits(1, :) = huge(0)
265 exp_limits(2, :) = -huge(0)
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
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
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
280 DEALLOCATE (basis_was_assoc)
284 cpabort(
"The angular momentum needed exceeds the value assumed when configuring libint.")
298 ALLOCATE (pnew(ndof))
300 ALLOCATE (deriv(ndof))
302 ALLOCATE (hessin(ndof, ndof))
305 hessin(i, i) = 1.0_dp
309 ALLOCATE (lower_b(ndof))
311 ALLOCATE (max_dev(ndof))
315 CALL init_transf(nkind, ri_basis_parameter, lower_b, max_dev, max_rel_dev)
318 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
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, &
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, &
344 DO iiter = 1, max_num_iter
345 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,I5)')
'OPTIMIZATION STEP NUMBER', iiter
349 CALL p2basis(nkind, ri_basis_parameter, lower_b, max_dev, pnew)
352 reset_boundary = .false.
355 DO iset = 1, ri_basis_parameter(ikind)%nset
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.
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)
378 hessin(i, i) = 1.0_dp
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, &
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, &
403 CALL p2basis(nkind, ri_basis_parameter, lower_b, max_dev, pnew)
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, &
421 IF (unit_nr > 0)
THEN
424 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=ri_aux_basis, &
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))
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'
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, &
467 hdg(:) = matmul(hessin, dg)
475 IF (
fac**2 > sumdg*sumxi*3.0e-8_dp)
THEN
478 dg(:) =
fac*xi - fad*hdg
481 hessin(i, j) = hessin(i, j) +
fac*xi(i)*xi(j) &
482 - fad*hdg(i)*hdg(j) &
487 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A)')
'Skip Hessian Update'
492 xi(:) = -matmul(hessin, g)
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, *)
501 DEALLOCATE (max_rel_dev)
513 DEALLOCATE (exp_limits)
515 IF (open_shell_case)
THEN
516 DEALLOCATE (integ_mp2_aa)
517 DEALLOCATE (integ_mp2_bb)
518 DEALLOCATE (integ_mp2_ab)
520 DEALLOCATE (integ_mp2)
522 DEALLOCATE (index_table_ri)
531 CALL timestop(handle)
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, &
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, &
593 REAL(kind=
dp),
INTENT(IN) :: eps
595 INTEGER,
INTENT(IN) :: nkind, natom, dimen, dimen_ri, homo, &
597 INTEGER,
DIMENSION(:),
INTENT(IN) :: kind_of
598 INTEGER,
ALLOCATABLE,
DIMENSION(:, :), &
599 INTENT(INOUT) :: index_table_ri
601 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
602 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: auto
603 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: c
605 POINTER :: ri_basis_parameter
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
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
618 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_energy_func_der'
620 INTEGER :: handle, ideriv, ikind, iset, nseta
621 REAL(kind=
dp) :: di, dri, emp2_ri, new_basis_val, &
623 REAL(kind=
dp),
VOLATILE :: step, temp
625 CALL timeset(routinen, handle)
633 nseta = ri_basis_parameter(ikind)%nset
637 IF (mod(ideriv, number_groups) /= color_sub) cycle
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
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.)
659 ri_basis_parameter(ikind)%zet(1, iset) = orig_basis_val
661 IF (para_env_sub%mepos == 0)
THEN
663 temp = temp/exp(di_ref)
664 deriv(ideriv) = log(temp)/step
670 CALL para_env%sum(deriv)
672 CALL timestop(handle)
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, &
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, &
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
731 TYPE(
mp2_type),
INTENT(INOUT) :: mp2_env
732 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: auto
733 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: c
735 POINTER :: ri_basis_parameter
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
744 INTEGER,
INTENT(IN) :: unit_nr
745 LOGICAL,
INTENT(IN) :: no_write
747 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_energy_func'
750 REAL(kind=
dp) :: di_aa, di_ab, di_bb, dri_aa, dri_ab, &
751 dri_bb, emp2_ri_aa, emp2_ri_ab, &
753 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: lai, lai_beta
755 CALL timeset(routinen, handle)
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)
767 IF (open_shell_case)
THEN
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, &
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)
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)
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
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, &
796 IF (.NOT. no_write)
THEN
797 IF (unit_nr > 0)
THEN
798 WRITE (unit_nr,
'(/,(T3,A,T56,F25.14))') &
801 WRITE (unit_nr,
'(T3,A,T56,ES25.10)') &
804 'DI/|Emp2| =', di/abs(emp2)
810 IF (open_shell_case)
DEALLOCATE (lai_beta)
812 CALL timestop(handle)
824 PURE SUBROUTINE init_transf(nkind, RI_basis_parameter, lower_B, max_dev, max_rel_dev)
825 INTEGER,
INTENT(IN) :: nkind
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
831 INTEGER :: ikind, ipos, iset
835 DO iset = 1, ri_basis_parameter(ikind)%nset
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)
852 SUBROUTINE p2basis(nkind, RI_basis_parameter, Lower_B, max_dev, p)
853 INTEGER,
INTENT(IN) :: nkind
855 POINTER :: ri_basis_parameter
856 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
857 INTENT(IN) :: lower_b, max_dev, p
859 INTEGER :: ikind, ipos, iset
860 REAL(kind=
dp) :: valout
864 DO iset = 1, ri_basis_parameter(ikind)%nset
866 valout = transf_val(lower_b(ipos), max_dev(ipos), p(ipos))
867 ri_basis_parameter(ikind)%zet(1, iset) = valout
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
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
894 INTEGER :: am_max, iexpo, ikind, ipos, ipos_p, &
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
907 old_lower_b = lower_b
908 old_max_dev = max_dev
909 old_max_rel_dev = max_rel_dev
913 DO iset = 1, ri_basis_parameter(ikind)%nset
915 old_expo(ipos) = ri_basis_parameter(ikind)%zet(1, iset)
930 am_max = maxval(ri_basis_parameter(ikind)%lmax(:))
931 ALLOCATE (nf_per_l(0:am_max))
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)
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
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))
954 ALLOCATE (has_to_be_changed(0:am_max))
955 has_to_be_changed = .false.
958 DO iexpo = 1, nf_per_l(la)
960 IF (abs(old_pnew(ipos_p)) >= pmax) pmax = abs(old_pnew(ipos_p))
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.
965 IF (pmax > reset_edge) has_to_be_changed(la) = .true.
969 IF (nf_per_l(la) == 1)
THEN
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
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)
979 IF (has_to_be_changed(la)) change_expo(ipos) = 1
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
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)
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
995 DEALLOCATE (has_to_be_changed)
997 DEALLOCATE (nf_per_l)
998 DEALLOCATE (max_min_exp_per_l)
1003 DO iset = 1, ri_basis_parameter(ikind)%nset
1005 ri_basis_parameter(ikind)%zet(1, iset) = new_expo(ipos)
1009 CALL init_transf(nkind, ri_basis_parameter, lower_b, max_dev, max_rel_dev)
1013 DO iset = 1, ri_basis_parameter(ikind)%nset
1015 IF (change_expo(ipos) == 0)
THEN
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)
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
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
1066 ALLOCATE (mat_ab(virtual, virtual_beta))
1073 ij_counter = ij_counter + 1
1074 IF (mod(ij_counter, para_env%num_pe) /= para_env%mepos) cycle
1076 mat_ab(:, :) = matmul(transpose(lai(:, :, i)), lai_beta(:, :, j))
1077 DO b = 1, virtual_beta
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)
1083 t_iajb = fact*abij(a, b, i, j)
1084 t_iajb_ri = fact*mat_ab(a, b)
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))
1089 emp2_ri = emp2_ri - t_iajb_ri*mat_ab(a, b)*fact2
1091 di = di - t_iajb*mat_ab(a, b)*fact2
1097 CALL para_env%sum(di)
1098 CALL para_env%sum(emp2_ri)
1100 dri = emp2 - emp2_ri
1101 di = 2.0d+00*di - emp2 - emp2_ri
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
1118 INTEGER,
INTENT(OUT) :: elements_ij_proc
1119 INTEGER,
ALLOCATABLE,
DIMENSION(:, :),
INTENT(OUT) :: ij_list_proc
1121 INTEGER :: i, ij_counter, j
1123 elements_ij_proc = 0
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
1132 ALLOCATE (ij_list_proc(elements_ij_proc, 2))
1135 elements_ij_proc = 0
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
1147 END SUBROUTINE calc_elem_ij_proc
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
1160 REAL(kind=
dp),
PARAMETER :: alpha = 2.633915794_dp
1163 valout = lower_b + max_dev/(1.0_dp + exp(-alpha*valin))
1175 SUBROUTINE generate_ri_init_basis(qs_env, mp2_env, nkind, max_rel_dev_output, basis_was_assoc)
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
1183 CHARACTER(LEN=*),
PARAMETER :: routinen =
'generate_RI_init_basis'
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, &
1188 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: l_expo, num_sgf_per_l, ordered_pos, &
1189 ri_l_expo, ri_num_sgf_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
1202 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1205 CALL timeset(routinen, handle)
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.
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)
1215 ALLOCATE (basis_was_assoc(nkind))
1216 basis_was_assoc = .false.
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.")
1226 atom_kind => qs_kind_set(ikind)
1228 CALL get_qs_kind(qs_kind=atom_kind, basis_set=orb_basis_a, basis_type=
"RI_AUX")
1230 basis_was_assoc(ikind) =
ASSOCIATED(orb_basis_a)
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
1243 max_am = maxval(la_max)
1246 ALLOCATE (ri_num_sgf_per_l(0:ri_max_am))
1247 ri_num_sgf_per_l = 0
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.")
1261 ALLOCATE (exp_tab(maxval(ri_num_sgf_per_l), 0:ri_max_am))
1263 DO iii = 0, ri_max_am
1266 DO la = la_min(iset), la_max(iset)
1267 IF (la /= iii) cycle
1269 exp_tab(iexpo, iii) = zet(1, iset)
1275 DO iii = 0, ri_max_am
1276 ALLOCATE (ordered_pos(ri_num_sgf_per_l(iii)))
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)
1282 ALLOCATE (ri_l_expo(ri_nset))
1284 ALLOCATE (ri_exponents(ri_nset))
1285 ri_exponents = 0.0_dp
1288 DO iii = 0, ri_max_am
1289 DO iexpo = 1, ri_num_sgf_per_l(iii)
1291 ri_l_expo(iset) = iii
1292 ri_exponents(iset) = exp_tab(iexpo, iii)
1295 DEALLOCATE (exp_tab)
1297 ALLOCATE (max_rel_dev(ri_nset))
1298 max_rel_dev = 1.0_dp
1300 DO iii = 0, ri_max_am
1301 IF (ri_num_sgf_per_l(iii) == 1)
THEN
1303 max_rel_dev(iset) = 0.35_dp
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)
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)
1315 max_rel_dev(:) = max_rel_dev(:)*0.9_dp
1316 DEALLOCATE (ri_num_sgf_per_l)
1323 CALL get_qs_kind(qs_kind=atom_kind, basis_set=orb_basis_a)
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
1335 nsgfa = orb_basis_a%nsgf
1337 max_am = maxval(la_max)
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
1346 IF (la /= nl(ishell, iset)) cycle
1347 nexpo_shell = nexpo_shell + npgfa(iset)
1352 ALLOCATE (exponents(nexpo_shell))
1354 ALLOCATE (l_expo(nexpo_shell))
1356 ALLOCATE (num_sgf_per_l(0:max_am))
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
1365 IF (la /= nl(ishell, iset)) cycle
1366 DO ipgf = 1, npgfa(iset)
1368 exponents(iexpo) = zet(ipgf, iset)
1371 num_sgf_per_l(la) = num_sgf_per_l(la) + 1
1376 ALLOCATE (exp_tab(nexpo_shell, nexpo_shell))
1378 ALLOCATE (l_tab(nexpo_shell, nexpo_shell))
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
1392 DEALLOCATE (exponents)
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)
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)
1407 DEALLOCATE (exp_tab)
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
1415 ALLOCATE (ri_num_sgf_per_l(0:max_am*2))
1416 ri_num_sgf_per_l = 0
1418 SELECT CASE (basis_quality)
1435 IF (external_num_of_func)
THEN
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
1443 DO la = 0, ri_max_am
1444 ri_num_sgf_per_l(la) = mp2_env%ri_opt_param%RI_nset_per_l(la)
1447 ri_num_sgf_per_l(0) = num_sgf_per_l(0)*2 + prog_func
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
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
1460 IF (ri_num_sgf_per_l(la) == 1)
EXIT
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
1474 IF (tot_num_exp_per_l(la) >= iii)
THEN
1475 iii = tot_num_exp_per_l(la)
1478 nsgfa_ri = nsgfa_ri + ri_num_sgf_per_l(la)*(la*2 + 1)
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
1486 ri_nset = sum(ri_num_sgf_per_l)
1488 ALLOCATE (ri_exponents(ri_nset))
1489 ri_exponents = 0.0_dp
1491 ALLOCATE (ri_l_expo(ri_nset))
1494 ALLOCATE (max_rel_dev(ri_nset))
1495 max_rel_dev = 1.0_dp
1498 DO la = 0, ri_max_am
1499 IF (ri_num_sgf_per_l(la) == 1)
THEN
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)
1505 max_rel_dev(iset) = (geom_fact - 1.0_dp)/(geom_fact + 1.0_dp)*0.9_dp
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)
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
1516 DEALLOCATE (num_sgf_per_l)
1517 DEALLOCATE (max_min_exp_l)
1518 DEALLOCATE (ri_num_sgf_per_l)
1524 CALL create_ri_basis(tmp_basis, ri_nset, ri_l_expo, ri_exponents)
1527 DEALLOCATE (ri_exponents)
1528 DEALLOCATE (ri_l_expo)
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
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)
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)
1545 DEALLOCATE (max_rel_dev)
1549 IF (external_num_of_func)
DEALLOCATE (mp2_env%ri_opt_param%RI_nset_per_l)
1551 CALL timestop(handle)
1553 END SUBROUTINE generate_ri_init_basis
1562 SUBROUTINE create_ri_basis(gto_basis_set, RI_nset, RI_l_expo, RI_exponents)
1564 INTEGER,
INTENT(IN) :: ri_nset
1565 INTEGER,
DIMENSION(:),
INTENT(IN) :: ri_l_expo
1566 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: ri_exponents
1568 CHARACTER(LEN=*),
PARAMETER :: routinen =
'create_ri_basis'
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
1578 CALL timeset(routinen, handle)
1579 NULLIFY (lmax, lmin, npgf, nshell, l, n, zet, gcc)
1591 gto_basis_set%nset = nset
1592 gto_basis_set%ncgf = ncgf
1593 gto_basis_set%nsgf = nsgf
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)
1617 CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
1633 lmin(iset) = ri_l_expo(iset)
1634 lmax(iset) = ri_l_expo(iset)
1637 maxl = max(maxl, lmax(iset))
1639 IF (npgf(iset) > maxpgf)
THEN
1642 CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
1645 DO lshell = lmin(iset), lmax(iset)
1646 nmin = n(1, iset) + lshell - lmin(iset)
1648 nshell(iset) = nshell(iset) + ishell
1649 IF (nshell(iset) > maxshell)
THEN
1650 maxshell = nshell(iset)
1653 CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
1656 n(nshell(iset) - ishell + i, iset) = nmin + i - 1
1657 l(nshell(iset) - ishell + i, iset) = lshell
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
1671 gto_basis_set%nset = 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)
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)
1695 DO ipgf = 1, npgf(iset)
1696 gto_basis_set%zet(ipgf, iset) = zet(ipgf, iset)
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)
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)
1731 maxco = max(maxco, npgf(iset)*
ncoset(lmax(iset)))
1734 gto_basis_set%ncgf = ncgf
1735 gto_basis_set%nsgf = nsgf
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)
1744 CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
1745 ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
1747 ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
1753 DO ishell = 1, nshell(iset)
1754 lshell = gto_basis_set%l(ishell, iset)
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)/))
1765 DO m = -lshell, lshell
1767 gto_basis_set%m(nsgf) = m
1768 gto_basis_set%sgf_symbol(nsgf) = &
1774 DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet)
1776 CALL timestop(handle)
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
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.
Defines the basic variable types.
integer, parameter, public dp
Interface to the Libint-Library or a c++ wrapper.
integer, parameter, public build_eri_size
Machine interface based on Fortran 2003 and POSIX.
integer, parameter, public default_output_unit
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
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.
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
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_pp, 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, harris_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, eeq, 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, zatom, 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_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_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
All kind of helpful little routines.
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.