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