1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate the interaction radii for the operator matrix calculation. 8!> \par History 9!> Joost VandeVondele : added exp_radius_very_extended 10!> 24.09.2002 overloaded init_interaction_radii for KG use (gt) 11!> 27.02.2004 integrated KG into QS version of init_interaction_radii (jgh) 12!> 07.03.2006 new routine to extend kind interaction radius for semi-empiric (jgh) 13!> \author MK (12.07.2000) 14! ************************************************************************************************** 15MODULE qs_interactions 16 17 USE ao_util, ONLY: exp_radius 18 USE atomic_kind_types, ONLY: atomic_kind_type,& 19 get_atomic_kind 20 USE basis_set_types, ONLY: get_gto_basis_set,& 21 gto_basis_set_type,& 22 set_gto_basis_set 23 USE cp_control_types, ONLY: qs_control_type,& 24 semi_empirical_control_type 25 USE cp_log_handling, ONLY: cp_get_default_logger,& 26 cp_logger_type 27 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 28 cp_print_key_unit_nr 29 USE cp_units, ONLY: cp_unit_from_cp2k 30 USE external_potential_types, ONLY: all_potential_type,& 31 get_potential,& 32 gth_potential_type,& 33 set_potential,& 34 sgp_potential_type 35 USE input_section_types, ONLY: section_vals_type,& 36 section_vals_val_get 37 USE kinds, ONLY: default_string_length,& 38 dp 39 USE paw_proj_set_types, ONLY: get_paw_proj_set,& 40 paw_proj_set_type,& 41 set_paw_proj_set 42 USE qs_kind_types, ONLY: get_qs_kind,& 43 qs_kind_type 44 USE string_utilities, ONLY: uppercase 45#include "./base/base_uses.f90" 46 47 IMPLICIT NONE 48 49 PRIVATE 50 51! *** Global parameters (only in this module) 52 53 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_interactions' 54 55! *** Public subroutines *** 56 57 PUBLIC :: init_interaction_radii, init_se_nlradius, init_interaction_radii_orb_basis 58 59 PUBLIC :: write_paw_radii, write_ppnl_radii, write_ppl_radii, write_core_charge_radii, & 60 write_pgf_orb_radii 61 62CONTAINS 63 64! ************************************************************************************************** 65!> \brief Initialize all the atomic kind radii for a given threshold value. 66!> \param qs_control ... 67!> \param qs_kind_set ... 68!> \date 24.09.2002 69!> \author GT 70!> \version 1.0 71! ************************************************************************************************** 72 SUBROUTINE init_interaction_radii(qs_control, qs_kind_set) 73 74 TYPE(qs_control_type), INTENT(IN) :: qs_control 75 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 76 77 CHARACTER(len=*), PARAMETER :: routineN = 'init_interaction_radii', & 78 routineP = moduleN//':'//routineN 79 80 INTEGER :: handle, i, iexp_ppl, ikind, ip, iprj_ppnl, j, l, lppl, lppnl, lprj, lprj_ppnl, & 81 maxl, n_local, nexp_lpot, nexp_lsd, nexp_ppl, nkind, nloc 82 INTEGER, DIMENSION(1:10) :: nrloc 83 INTEGER, DIMENSION(:), POINTER :: nct_lpot, nct_lsd, nprj, nprj_ppnl 84 LOGICAL :: ecp_local, llsd, lpot, paw_atom 85 LOGICAL, DIMENSION(0:5) :: is_nonlocal 86 REAL(dp), DIMENSION(1:10) :: aloc, bloc 87 REAL(KIND=dp) :: alpha_core_charge, alpha_ppl, ccore_charge, cerf_ppl, core_charge_radius, & 88 cval, ppl_radius, ppnl_radius, rcprj, zeta 89 REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, a_nonlocal, alpha_lpot, & 90 alpha_lsd, alpha_ppnl, c_local, & 91 cexp_ppl 92 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj_ppnl, cval_lpot, cval_lsd, rzetprj, & 93 zet 94 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nonlocal 95 TYPE(all_potential_type), POINTER :: all_potential 96 TYPE(gth_potential_type), POINTER :: gth_potential 97 TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, aux_gw_basis, & 98 harris_basis, lri_basis, mao_basis, orb_basis_set, ri_aux_basis_set, ri_basis, & 99 ri_xas_basis, soft_basis 100 TYPE(paw_proj_set_type), POINTER :: paw_proj_set 101 TYPE(sgp_potential_type), POINTER :: sgp_potential 102 103 CALL timeset(routineN, handle) 104 105 NULLIFY (all_potential, gth_potential, sgp_potential, orb_basis_set, aux_basis_set, soft_basis, & 106 lri_basis, mao_basis, paw_proj_set, harris_basis, ri_xas_basis) 107 NULLIFY (nprj_ppnl, nprj) 108 NULLIFY (alpha_ppnl, cexp_ppl, cprj_ppnl, zet) 109 110 nkind = SIZE(qs_kind_set) 111 nexp_ppl = 0 112 113 DO ikind = 1, nkind 114 115 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB") 116 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") 117 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT") 118 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis, basis_type="LRI_AUX") 119 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_HXC") 120 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_aux_basis_set, basis_type="RI_AUX") 121 CALL get_qs_kind(qs_kind_set(ikind), basis_set=mao_basis, basis_type="MAO") 122 CALL get_qs_kind(qs_kind_set(ikind), basis_set=harris_basis, basis_type="HARRIS") 123 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_gw_basis, basis_type="AUX_GW") 124 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_xas_basis, basis_type="RI_XAS") 125 CALL get_qs_kind(qs_kind_set(ikind), soft_basis_set=soft_basis) 126 CALL get_qs_kind(qs_kind_set(ikind), & 127 paw_proj_set=paw_proj_set, & 128 paw_atom=paw_atom, & 129 all_potential=all_potential, & 130 gth_potential=gth_potential, & 131 sgp_potential=sgp_potential) 132 133 ! Calculate the orbital basis function radii *** 134 ! For ALL electron this has to come before the calculation of the 135 ! radii used to calculate the 3c lists. 136 ! Indeed, there the radii of the GTO primitives are needed 137 IF (ASSOCIATED(orb_basis_set)) THEN 138 CALL init_interaction_radii_orb_basis(orb_basis_set, qs_control%eps_pgf_orb, & 139 qs_control%eps_kg_orb) 140 END IF 141 142 IF (ASSOCIATED(all_potential)) THEN 143 144 CALL get_potential(potential=all_potential, & 145 alpha_core_charge=alpha_core_charge, & 146 ccore_charge=ccore_charge) 147 148 ! Calculate the radius of the core charge distribution 149 150 core_charge_radius = exp_radius(0, alpha_core_charge, & 151 qs_control%eps_core_charge, & 152 ccore_charge) 153 154 CALL set_potential(potential=all_potential, & 155 core_charge_radius=core_charge_radius) 156 157 ELSE IF (ASSOCIATED(gth_potential)) THEN 158 159 CALL get_potential(potential=gth_potential, & 160 alpha_core_charge=alpha_core_charge, & 161 ccore_charge=ccore_charge, & 162 alpha_ppl=alpha_ppl, & 163 cerf_ppl=cerf_ppl, & 164 nexp_ppl=nexp_ppl, & 165 cexp_ppl=cexp_ppl, & 166 lpot_present=lpot, & 167 lsd_present=llsd, & 168 lppnl=lppnl, & 169 alpha_ppnl=alpha_ppnl, & 170 nprj_ppnl=nprj_ppnl, & 171 cprj_ppnl=cprj_ppnl) 172 173 ! Calculate the radius of the core charge distribution 174 core_charge_radius = exp_radius(0, alpha_core_charge, & 175 qs_control%eps_core_charge, & 176 ccore_charge) 177 178 ! Calculate the radii of the local part 179 ! of the Goedecker pseudopotential (GTH) 180 ppl_radius = exp_radius(0, alpha_ppl, qs_control%eps_ppl, cerf_ppl) 181 182 DO iexp_ppl = 1, nexp_ppl 183 lppl = 2*(iexp_ppl - 1) 184 ppl_radius = MAX(ppl_radius, & 185 exp_radius(lppl, alpha_ppl, & 186 qs_control%eps_ppl, & 187 cexp_ppl(iexp_ppl), & 188 rlow=ppl_radius)) 189 END DO 190 IF (lpot) THEN 191 ! extended local potential 192 CALL get_potential(potential=gth_potential, & 193 nexp_lpot=nexp_lpot, & 194 alpha_lpot=alpha_lpot, & 195 nct_lpot=nct_lpot, & 196 cval_lpot=cval_lpot) 197 DO j = 1, nexp_lpot 198 DO i = 1, nct_lpot(j) 199 lppl = 2*(i - 1) 200 ppl_radius = MAX(ppl_radius, & 201 exp_radius(lppl, alpha_lpot(j), qs_control%eps_ppl, & 202 cval_lpot(i, j), rlow=ppl_radius)) 203 END DO 204 END DO 205 END IF 206 IF (llsd) THEN 207 ! lsd local potential 208 CALL get_potential(potential=gth_potential, & 209 nexp_lsd=nexp_lsd, & 210 alpha_lsd=alpha_lsd, & 211 nct_lsd=nct_lsd, & 212 cval_lsd=cval_lsd) 213 DO j = 1, nexp_lsd 214 DO i = 1, nct_lsd(j) 215 lppl = 2*(i - 1) 216 ppl_radius = MAX(ppl_radius, & 217 exp_radius(lppl, alpha_lsd(j), qs_control%eps_ppl, & 218 cval_lsd(i, j), rlow=ppl_radius)) 219 END DO 220 END DO 221 END IF 222 223 ! Calculate the radii of the non-local part 224 ! of the Goedecker pseudopotential (GTH) 225 ppnl_radius = 0.0_dp 226 DO l = 0, lppnl 227 DO iprj_ppnl = 1, nprj_ppnl(l) 228 lprj_ppnl = l + 2*(iprj_ppnl - 1) 229 ppnl_radius = MAX(ppnl_radius, & 230 exp_radius(lprj_ppnl, alpha_ppnl(l), & 231 qs_control%eps_pgf_orb, & 232 cprj_ppnl(iprj_ppnl, l), & 233 rlow=ppnl_radius)) 234 END DO 235 END DO 236 CALL set_potential(potential=gth_potential, & 237 core_charge_radius=core_charge_radius, & 238 ppl_radius=ppl_radius, & 239 ppnl_radius=ppnl_radius) 240 241 ELSE IF (ASSOCIATED(sgp_potential)) THEN 242 243 ! Calculate the radius of the core charge distribution 244 CALL get_potential(potential=sgp_potential, & 245 alpha_core_charge=alpha_core_charge, & 246 ccore_charge=ccore_charge) 247 core_charge_radius = exp_radius(0, alpha_core_charge, & 248 qs_control%eps_core_charge, & 249 ccore_charge) 250 ! Calculate the radii of the local part 251 ppl_radius = core_charge_radius 252 CALL get_potential(potential=sgp_potential, ecp_local=ecp_local) 253 IF (ecp_local) THEN 254 CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc) 255 DO i = 1, nloc 256 lppl = MAX(0, nrloc(i) - 2) 257 ppl_radius = MAX(ppl_radius, & 258 exp_radius(lppl, bloc(i), qs_control%eps_ppl, aloc(i), rlow=ppl_radius)) 259 END DO 260 ELSE 261 CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local) 262 DO i = 1, n_local 263 ppl_radius = MAX(ppl_radius, & 264 exp_radius(0, a_local(i), qs_control%eps_ppl, c_local(i), rlow=ppl_radius)) 265 END DO 266 END IF 267 ! Calculate the radii of the non-local part 268 ppnl_radius = 0.0_dp 269 CALL get_potential(potential=sgp_potential, lmax=lppnl, n_nonlocal=nloc) 270 CALL get_potential(potential=sgp_potential, is_nonlocal=is_nonlocal, & 271 a_nonlocal=a_nonlocal, c_nonlocal=c_nonlocal) 272 DO l = 0, lppnl 273 IF (is_nonlocal(l)) THEN 274 DO i = 1, nloc 275 cval = MAXVAL(ABS(c_nonlocal(i, :, l))) 276 ppnl_radius = MAX(ppnl_radius, & 277 exp_radius(l, a_nonlocal(i), qs_control%eps_pgf_orb, cval, rlow=ppnl_radius)) 278 END DO 279 END IF 280 END DO 281 CALL set_potential(potential=sgp_potential, & 282 core_charge_radius=core_charge_radius, & 283 ppl_radius=ppl_radius, & 284 ppnl_radius=ppnl_radius) 285 END IF 286 287 ! Calculate the aux fit orbital basis function radii 288 IF (ASSOCIATED(aux_fit_basis_set)) THEN 289 CALL init_interaction_radii_orb_basis(aux_fit_basis_set, qs_control%eps_pgf_orb, & 290 qs_control%eps_kg_orb) 291 END IF 292 293 ! Calculate the aux orbital basis function radii 294 IF (ASSOCIATED(aux_basis_set)) THEN 295 CALL init_interaction_radii_orb_basis(aux_basis_set, qs_control%eps_pgf_orb) 296 END IF 297 298 ! Calculate the RI aux orbital basis function radii 299 IF (ASSOCIATED(RI_aux_basis_set)) THEN 300 CALL init_interaction_radii_orb_basis(RI_aux_basis_set, qs_control%eps_pgf_orb) 301 END IF 302 303 ! Calculate the MAO orbital basis function radii 304 IF (ASSOCIATED(mao_basis)) THEN 305 CALL init_interaction_radii_orb_basis(mao_basis, qs_control%eps_pgf_orb) 306 END IF 307 308 ! Calculate the HARRIS orbital basis function radii 309 IF (ASSOCIATED(harris_basis)) THEN 310 CALL init_interaction_radii_orb_basis(harris_basis, qs_control%eps_pgf_orb) 311 END IF 312 313 ! Calculate the AUX_GW orbital basis function radii 314 IF (ASSOCIATED(aux_gw_basis)) THEN 315 CALL init_interaction_radii_orb_basis(aux_gw_basis, qs_control%eps_pgf_orb) 316 END IF 317 318 ! Calculate the soft orbital basis function radii 319 IF (ASSOCIATED(soft_basis)) THEN 320 IF (paw_atom) THEN 321 CALL init_interaction_radii_orb_basis(soft_basis, qs_control%eps_pgf_orb) 322 END IF 323 END IF 324 325 ! Calculate the lri basis function radii 326 IF (ASSOCIATED(lri_basis)) THEN 327 CALL init_interaction_radii_orb_basis(lri_basis, qs_control%eps_pgf_orb) 328 END IF 329 IF (ASSOCIATED(ri_basis)) THEN 330 CALL init_interaction_radii_orb_basis(ri_basis, qs_control%eps_pgf_orb) 331 END IF 332 IF (ASSOCIATED(ri_xas_basis)) THEN 333 CALL init_interaction_radii_orb_basis(ri_xas_basis, qs_control%eps_pgf_orb) 334 END IF 335 336 ! Calculate the paw projector basis function radii 337 IF (ASSOCIATED(paw_proj_set)) THEN 338 CALL get_paw_proj_set(paw_proj_set=paw_proj_set, & 339 nprj=nprj, & 340 maxl=maxl, & 341 zetprj=zet, & 342 rzetprj=rzetprj) 343 rcprj = 0.0_dp 344 rzetprj = 0.0_dp 345 DO lprj = 0, maxl 346 DO ip = 1, nprj(lprj) 347 zeta = zet(ip, lprj) 348 rzetprj(ip, lprj) = MAX(rzetprj(ip, lprj), & 349 exp_radius(lprj, zeta, qs_control%eps_pgf_orb, & 350 1.0_dp, rlow=rzetprj(ip, lprj))) 351 rcprj = MAX(rcprj, rzetprj(ip, lprj)) 352 ENDDO 353 ENDDO 354 CALL set_paw_proj_set(paw_proj_set=paw_proj_set, & 355 rzetprj=rzetprj, & 356 rcprj=rcprj) 357 358 END IF 359 END DO ! ikind 360 361 CALL timestop(handle) 362 363 END SUBROUTINE init_interaction_radii 364 365! ************************************************************************************************** 366!> \brief ... 367!> \param orb_basis_set ... 368!> \param eps_pgf_orb ... 369!> \param eps_pgf_short ... 370! ************************************************************************************************** 371 SUBROUTINE init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short) 372 373 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 374 REAL(KIND=dp), INTENT(IN) :: eps_pgf_orb 375 REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_pgf_short 376 377 CHARACTER(len=*), PARAMETER :: routineN = 'init_interaction_radii_orb_basis', & 378 routineP = moduleN//':'//routineN 379 380 INTEGER :: ipgf, iset, ishell, l, nset 381 INTEGER, DIMENSION(:), POINTER :: npgf, nshell 382 INTEGER, DIMENSION(:, :), POINTER :: lshell 383 REAL(KIND=dp) :: eps_short, gcca, kind_radius, & 384 short_radius, zeta 385 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius 386 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius, zet 387 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc 388 389 IF (ASSOCIATED(orb_basis_set)) THEN 390 391 IF (PRESENT(eps_pgf_short)) THEN 392 eps_short = eps_pgf_short 393 ELSE 394 eps_short = eps_pgf_orb 395 END IF 396 397 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 398 nset=nset, & 399 nshell=nshell, & 400 npgf=npgf, & 401 l=lshell, & 402 zet=zet, & 403 gcc=gcc, & 404 pgf_radius=pgf_radius, & 405 set_radius=set_radius) 406 407 kind_radius = 0.0_dp 408 short_radius = 0.0_dp 409 410 DO iset = 1, nset 411 set_radius(iset) = 0.0_dp 412 DO ipgf = 1, npgf(iset) 413 pgf_radius(ipgf, iset) = 0.0_dp 414 DO ishell = 1, nshell(iset) 415 l = lshell(ishell, iset) 416 gcca = gcc(ipgf, ishell, iset) 417 zeta = zet(ipgf, iset) 418 pgf_radius(ipgf, iset) = MAX(pgf_radius(ipgf, iset), & 419 exp_radius(l, zeta, eps_pgf_orb, gcca, & 420 rlow=pgf_radius(ipgf, iset))) 421 short_radius = MAX(short_radius, & 422 exp_radius(l, zeta, eps_short, gcca, & 423 rlow=short_radius)) 424 END DO 425 set_radius(iset) = MAX(set_radius(iset), pgf_radius(ipgf, iset)) 426 END DO 427 kind_radius = MAX(kind_radius, set_radius(iset)) 428 END DO 429 430 CALL set_gto_basis_set(gto_basis_set=orb_basis_set, & 431 pgf_radius=pgf_radius, & 432 set_radius=set_radius, & 433 kind_radius=kind_radius, & 434 short_kind_radius=short_radius) 435 436 END IF 437 438 END SUBROUTINE init_interaction_radii_orb_basis 439 440! ************************************************************************************************** 441!> \brief ... 442!> \param se_control ... 443!> \param atomic_kind_set ... 444!> \param qs_kind_set ... 445!> \param subsys_section ... 446! ************************************************************************************************** 447 SUBROUTINE init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, & 448 subsys_section) 449 450 TYPE(semi_empirical_control_type), POINTER :: se_control 451 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 452 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 453 TYPE(section_vals_type), POINTER :: subsys_section 454 455 CHARACTER(len=*), PARAMETER :: routineN = 'init_se_nlradius', & 456 routineP = moduleN//':'//routineN 457 458 INTEGER :: ikind, nkind 459 REAL(KIND=dp) :: kind_radius 460 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 461 TYPE(qs_kind_type), POINTER :: qs_kind 462 463 NULLIFY (orb_basis_set, qs_kind) 464 465 nkind = SIZE(qs_kind_set) 466 467 DO ikind = 1, nkind 468 469 qs_kind => qs_kind_set(ikind) 470 471 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set) 472 473 IF (ASSOCIATED(orb_basis_set)) THEN 474 475 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 476 kind_radius=kind_radius) 477 478 kind_radius = MAX(kind_radius, se_control%cutoff_exc) 479 480 CALL set_gto_basis_set(gto_basis_set=orb_basis_set, & 481 kind_radius=kind_radius) 482 483 END IF 484 485 END DO 486 487 CALL write_kind_radii(atomic_kind_set, qs_kind_set, subsys_section) 488 489 END SUBROUTINE init_se_nlradius 490 491! ************************************************************************************************** 492!> \brief ... 493!> \param atomic_kind_set ... 494!> \param qs_kind_set ... 495!> \param subsys_section ... 496! ************************************************************************************************** 497 SUBROUTINE write_kind_radii(atomic_kind_set, qs_kind_set, subsys_section) 498 499 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 500 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 501 TYPE(section_vals_type), POINTER :: subsys_section 502 503 CHARACTER(LEN=10) :: bas 504 CHARACTER(LEN=default_string_length) :: name, unit_str 505 INTEGER :: ikind, nkind, output_unit 506 REAL(KIND=dp) :: conv, kind_radius 507 TYPE(cp_logger_type), POINTER :: logger 508 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 509 510 NULLIFY (logger) 511 logger => cp_get_default_logger() 512 NULLIFY (orb_basis_set) 513 bas = "ORBITAL " 514 515 nkind = SIZE(atomic_kind_set) 516! *** Print the kind radii *** 517 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 518 "PRINT%RADII/KIND_RADII", extension=".Log") 519 CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) 520 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 521 IF (output_unit > 0) THEN 522 WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & 523 "RADII: "//TRIM(bas)//" BASIS in "//TRIM(unit_str), & 524 "Kind", "Label", "Radius" 525 DO ikind = 1, nkind 526 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 527 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 528 IF (ASSOCIATED(orb_basis_set)) THEN 529 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 530 kind_radius=kind_radius) 531 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & 532 ikind, name, kind_radius*conv 533 ELSE 534 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T72,A9)") & 535 ikind, name, "no basis" 536 END IF 537 END DO 538 END IF 539 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 540 "PRINT%RADII/KIND_RADII") 541 542 END SUBROUTINE write_kind_radii 543 544! ************************************************************************************************** 545!> \brief Write the radii of the core charge distributions to the output 546!> unit. 547!> \param atomic_kind_set ... 548!> \param qs_kind_set ... 549!> \param subsys_section ... 550!> \date 15.09.2000 551!> \author MK 552!> \version 1.0 553! ************************************************************************************************** 554 SUBROUTINE write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section) 555 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 556 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 557 TYPE(section_vals_type), POINTER :: subsys_section 558 559 CHARACTER(LEN=default_string_length) :: name, unit_str 560 INTEGER :: ikind, nkind, output_unit 561 REAL(KIND=dp) :: conv, core_charge_radius 562 TYPE(all_potential_type), POINTER :: all_potential 563 TYPE(cp_logger_type), POINTER :: logger 564 TYPE(gth_potential_type), POINTER :: gth_potential 565 566 NULLIFY (logger) 567 logger => cp_get_default_logger() 568 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 569 "PRINT%RADII/CORE_CHARGE_RADII", extension=".Log") 570 CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) 571 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 572 IF (output_unit > 0) THEN 573 nkind = SIZE(atomic_kind_set) 574 WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & 575 "RADII: CORE CHARGE DISTRIBUTIONS in "// & 576 TRIM(unit_str), "Kind", "Label", "Radius" 577 DO ikind = 1, nkind 578 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 579 CALL get_qs_kind(qs_kind_set(ikind), & 580 all_potential=all_potential, gth_potential=gth_potential) 581 582 IF (ASSOCIATED(all_potential)) THEN 583 CALL get_potential(potential=all_potential, & 584 core_charge_radius=core_charge_radius) 585 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & 586 ikind, name, core_charge_radius*conv 587 ELSE IF (ASSOCIATED(gth_potential)) THEN 588 CALL get_potential(potential=gth_potential, & 589 core_charge_radius=core_charge_radius) 590 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & 591 ikind, name, core_charge_radius*conv 592 END IF 593 END DO 594 END IF 595 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 596 "PRINT%RADII/CORE_CHARGE_RADII") 597 END SUBROUTINE write_core_charge_radii 598 599! ************************************************************************************************** 600!> \brief Write the orbital basis function radii to the output unit. 601!> \param basis ... 602!> \param atomic_kind_set ... 603!> \param qs_kind_set ... 604!> \param subsys_section ... 605!> \date 05.06.2000 606!> \author MK 607!> \version 1.0 608! ************************************************************************************************** 609 SUBROUTINE write_pgf_orb_radii(basis, atomic_kind_set, qs_kind_set, subsys_section) 610 611 CHARACTER(len=*) :: basis 612 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 613 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 614 TYPE(section_vals_type), POINTER :: subsys_section 615 616 CHARACTER(len=*), PARAMETER :: routineN = 'write_pgf_orb_radii', & 617 routineP = moduleN//':'//routineN 618 619 CHARACTER(LEN=10) :: bas 620 CHARACTER(LEN=default_string_length) :: name, unit_str 621 INTEGER :: ikind, ipgf, iset, nkind, nset, & 622 output_unit 623 INTEGER, DIMENSION(:), POINTER :: npgf 624 REAL(KIND=dp) :: conv, kind_radius, short_kind_radius 625 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius 626 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius 627 TYPE(cp_logger_type), POINTER :: logger 628 TYPE(gto_basis_set_type), POINTER :: aux_basis_set, lri_basis_set, & 629 orb_basis_set 630 631 NULLIFY (logger) 632 logger => cp_get_default_logger() 633 NULLIFY (aux_basis_set, orb_basis_set, lri_basis_set) 634 bas = " " 635 bas(1:3) = basis(1:3) 636 CALL uppercase(bas) 637 IF (bas(1:3) == "ORB") THEN 638 bas = "ORBITAL " 639 ELSE IF (bas(1:3) == "AUX") THEN 640 bas = "AUXILLIARY" 641 ELSE IF (bas(1:3) == "LRI") THEN 642 bas = "LOCAL RI" 643 ELSE 644 CPABORT("Undefined basis set type") 645 END IF 646 647 nkind = SIZE(qs_kind_set) 648 649! *** Print the kind radii *** 650 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 651 "PRINT%RADII/KIND_RADII", extension=".Log") 652 CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) 653 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 654 IF (output_unit > 0) THEN 655 WRITE (UNIT=output_unit, FMT="(/,T2,A,T46,A,T53,A,T63,A,T71,A)") & 656 "RADII: "//TRIM(bas)//" BASIS in "//TRIM(unit_str), & 657 "Kind", "Label", "Radius", "OCE Radius" 658 DO ikind = 1, nkind 659 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 660 IF (bas(1:3) == "ORB") THEN 661 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 662 IF (ASSOCIATED(orb_basis_set)) THEN 663 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 664 kind_radius=kind_radius, & 665 short_kind_radius=short_kind_radius) 666 END IF 667 ELSE IF (bas(1:3) == "AUX") THEN 668 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") 669 IF (ASSOCIATED(aux_basis_set)) THEN 670 CALL get_gto_basis_set(gto_basis_set=aux_basis_set, & 671 kind_radius=kind_radius) 672 END IF 673 ELSE IF (bas(1:3) == "LOC") THEN 674 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX") 675 IF (ASSOCIATED(lri_basis_set)) THEN 676 CALL get_gto_basis_set(gto_basis_set=lri_basis_set, & 677 kind_radius=kind_radius) 678 END IF 679 ELSE 680 CPABORT("Undefined basis set type") 681 END IF 682 IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN 683 WRITE (UNIT=output_unit, FMT="(T45,I5,T53,A5,T57,F12.6,T69,F12.6)") & 684 ikind, name, kind_radius*conv, short_kind_radius*conv 685 ELSE 686 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T72,A9)") & 687 ikind, name, "no basis" 688 END IF 689 END DO 690 END IF 691 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 692 "PRINT%RADII/KIND_RADII") 693 694! *** Print the shell set radii *** 695 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 696 "PRINT%RADII/SET_RADII", extension=".Log") 697 IF (output_unit > 0) THEN 698 WRITE (UNIT=output_unit, FMT="(/,T2,A,T51,A,T57,A,T65,A,T75,A)") & 699 "RADII: SHELL SETS OF "//TRIM(bas)//" BASIS in "// & 700 TRIM(unit_str), "Kind", "Label", "Set", "Radius" 701 DO ikind = 1, nkind 702 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 703 IF (bas(1:3) == "ORB") THEN 704 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 705 IF (ASSOCIATED(orb_basis_set)) THEN 706 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 707 nset=nset, & 708 set_radius=set_radius) 709 END IF 710 ELSE IF (bas(1:3) == "AUX") THEN 711 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") 712 IF (ASSOCIATED(aux_basis_set)) THEN 713 CALL get_gto_basis_set(gto_basis_set=aux_basis_set, & 714 nset=nset, & 715 set_radius=set_radius) 716 END IF 717 ELSE IF (bas(1:3) == "LOC") THEN 718 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX") 719 IF (ASSOCIATED(lri_basis_set)) THEN 720 CALL get_gto_basis_set(gto_basis_set=lri_basis_set, & 721 nset=nset, & 722 set_radius=set_radius) 723 END IF 724 ELSE 725 CPABORT("Undefined basis set type") 726 END IF 727 IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN 728 WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,(T63,I5,T69,F12.6))") & 729 ikind, name, (iset, set_radius(iset)*conv, iset=1, nset) 730 ELSE 731 WRITE (UNIT=output_unit, FMT="(T50,I5,T58,A5,T73,A8)") & 732 ikind, name, "no basis" 733 END IF 734 END DO 735 END IF 736 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 737 "PRINT%RADII/SET_RADII") 738! *** Print the primitive Gaussian function radii *** 739 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 740 "PRINT%RADII/PGF_RADII", extension=".Log") 741 IF (output_unit > 0) THEN 742 WRITE (UNIT=output_unit, FMT="(/,T2,A,T51,A,T57,A,T65,A,T75,A)") & 743 "RADII: PRIMITIVE GAUSSIANS OF "//TRIM(bas)//" BASIS in "// & 744 TRIM(unit_str), "Kind", "Label", "Set", "Radius" 745 DO ikind = 1, nkind 746 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 747 IF (bas(1:3) == "ORB") THEN 748 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 749 IF (ASSOCIATED(orb_basis_set)) THEN 750 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 751 nset=nset, & 752 npgf=npgf, & 753 pgf_radius=pgf_radius) 754 END IF 755 ELSE IF (bas(1:3) == "AUX") THEN 756 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") 757 IF (ASSOCIATED(aux_basis_set)) THEN 758 CALL get_gto_basis_set(gto_basis_set=aux_basis_set, & 759 nset=nset, & 760 npgf=npgf, & 761 pgf_radius=pgf_radius) 762 END IF 763 ELSE IF (bas(1:3) == "LOC") THEN 764 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX") 765 IF (ASSOCIATED(lri_basis_set)) THEN 766 CALL get_gto_basis_set(gto_basis_set=lri_basis_set, & 767 nset=nset, & 768 npgf=npgf, & 769 pgf_radius=pgf_radius) 770 END IF 771 ELSE 772 CPABORT("Undefined basis type") 773 END IF 774 IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set) .OR. & 775 ASSOCIATED(lri_basis_set)) THEN 776 DO iset = 1, nset 777 WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,T63,I5,(T69,F12.6))") & 778 ikind, name, iset, & 779 (pgf_radius(ipgf, iset)*conv, ipgf=1, npgf(iset)) 780 END DO 781 ELSE 782 WRITE (UNIT=output_unit, FMT="(T50,I5,T58,A5,T73,A8)") & 783 ikind, name, "no basis" 784 END IF 785 END DO 786 END IF 787 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 788 "PRINT%RADII/PGF_RADII") 789 END SUBROUTINE write_pgf_orb_radii 790 791! ************************************************************************************************** 792!> \brief Write the radii of the exponential functions of the Goedecker 793!> pseudopotential (GTH, local part) to the logical unit number 794!> "output_unit". 795!> \param atomic_kind_set ... 796!> \param qs_kind_set ... 797!> \param subsys_section ... 798!> \date 06.10.2000 799!> \author MK 800!> \version 1.0 801! ************************************************************************************************** 802 SUBROUTINE write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section) 803 804 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 805 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 806 TYPE(section_vals_type), POINTER :: subsys_section 807 808 CHARACTER(LEN=default_string_length) :: name, unit_str 809 INTEGER :: ikind, nkind, output_unit 810 REAL(KIND=dp) :: conv, ppl_radius 811 TYPE(cp_logger_type), POINTER :: logger 812 TYPE(gth_potential_type), POINTER :: gth_potential 813 814 NULLIFY (logger) 815 logger => cp_get_default_logger() 816 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 817 "PRINT%RADII/PPL_RADII", extension=".Log") 818 CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) 819 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 820 IF (output_unit > 0) THEN 821 nkind = SIZE(atomic_kind_set) 822 WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & 823 "RADII: LOCAL PART OF GTH/ELP PP in "// & 824 TRIM(unit_str), "Kind", "Label", "Radius" 825 DO ikind = 1, nkind 826 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 827 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential) 828 IF (ASSOCIATED(gth_potential)) THEN 829 CALL get_potential(potential=gth_potential, & 830 ppl_radius=ppl_radius) 831 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & 832 ikind, name, ppl_radius*conv 833 END IF 834 END DO 835 END IF 836 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 837 "PRINT%RADII/PPL_RADII") 838 END SUBROUTINE write_ppl_radii 839 840! ************************************************************************************************** 841!> \brief Write the radii of the projector functions of the Goedecker 842!> pseudopotential (GTH, non-local part) to the logical unit number 843!> "output_unit". 844!> \param atomic_kind_set ... 845!> \param qs_kind_set ... 846!> \param subsys_section ... 847!> \date 06.10.2000 848!> \author MK 849!> \version 1.0 850! ************************************************************************************************** 851 SUBROUTINE write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section) 852 853 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 854 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 855 TYPE(section_vals_type), POINTER :: subsys_section 856 857 CHARACTER(LEN=default_string_length) :: name, unit_str 858 INTEGER :: ikind, nkind, output_unit 859 REAL(KIND=dp) :: conv, ppnl_radius 860 TYPE(cp_logger_type), POINTER :: logger 861 TYPE(gth_potential_type), POINTER :: gth_potential 862 863 NULLIFY (logger) 864 logger => cp_get_default_logger() 865 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 866 "PRINT%RADII/PPNL_RADII", extension=".Log") 867 CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) 868 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 869 IF (output_unit > 0) THEN 870 nkind = SIZE(atomic_kind_set) 871 WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & 872 "RADII: NON-LOCAL PART OF GTH PP in "// & 873 TRIM(unit_str), "Kind", "Label", "Radius" 874 DO ikind = 1, nkind 875 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 876 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential) 877 IF (ASSOCIATED(gth_potential)) THEN 878 CALL get_potential(potential=gth_potential, & 879 ppnl_radius=ppnl_radius) 880 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & 881 ikind, name, ppnl_radius*conv 882 END IF 883 END DO 884 END IF 885 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 886 "PRINT%RADII/PPNL_RADII") 887 END SUBROUTINE write_ppnl_radii 888 889! ************************************************************************************************** 890!> \brief Write the radii of the one center projector 891!> \param atomic_kind_set ... 892!> \param qs_kind_set ... 893!> \param subsys_section ... 894!> \version 1.0 895! ************************************************************************************************** 896 SUBROUTINE write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section) 897 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 898 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 899 TYPE(section_vals_type), POINTER :: subsys_section 900 901 CHARACTER(LEN=default_string_length) :: name, unit_str 902 INTEGER :: ikind, nkind, output_unit 903 LOGICAL :: paw_atom 904 REAL(KIND=dp) :: conv, rcprj 905 TYPE(cp_logger_type), POINTER :: logger 906 TYPE(paw_proj_set_type), POINTER :: paw_proj_set 907 908 NULLIFY (logger) 909 logger => cp_get_default_logger() 910 output_unit = cp_print_key_unit_nr(logger, subsys_section, & 911 "PRINT%RADII/GAPW_PRJ_RADII", extension=".Log") 912 CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) 913 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 914 IF (output_unit > 0) THEN 915 nkind = SIZE(qs_kind_set) 916 WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & 917 "RADII: ONE CENTER PROJECTORS in "// & 918 TRIM(unit_str), "Kind", "Label", "Radius" 919 DO ikind = 1, nkind 920 CALL get_atomic_kind(atomic_kind_set(ikind), name=name) 921 CALL get_qs_kind(qs_kind_set(ikind), & 922 paw_atom=paw_atom, paw_proj_set=paw_proj_set) 923 IF (paw_atom .AND. ASSOCIATED(paw_proj_set)) THEN 924 CALL get_paw_proj_set(paw_proj_set=paw_proj_set, & 925 rcprj=rcprj) 926 WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & 927 ikind, name, rcprj*conv 928 END IF 929 END DO 930 END IF 931 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & 932 "PRINT%RADII/GAPW_PRJ_RADII") 933 END SUBROUTINE write_paw_radii 934 935END MODULE qs_interactions 936