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