LCOV - code coverage report
Current view: top level - src - lri_optimize_ri_basis.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.3 % 299 297
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 12 12

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

Generated by: LCOV version 2.0-1