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