!--------------------------------------------------------------------------------------------------! ! CP2K: A general program to perform molecular dynamics simulations ! ! Copyright (C) 2000 - 2019 CP2K developers group ! !--------------------------------------------------------------------------------------------------! ! ************************************************************************************************** !> \brief Calculate the interaction radii for the operator matrix calculation. !> \par History !> Joost VandeVondele : added exp_radius_very_extended !> 24.09.2002 overloaded init_interaction_radii for KG use (gt) !> 27.02.2004 integrated KG into QS version of init_interaction_radii (jgh) !> 07.03.2006 new routine to extend kind interaction radius for semi-empiric (jgh) !> \author MK (12.07.2000) ! ************************************************************************************************** MODULE qs_interactions USE ao_util, ONLY: exp_radius USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind USE basis_set_types, ONLY: get_gto_basis_set,& gto_basis_set_type,& set_gto_basis_set USE cp_control_types, ONLY: qs_control_type,& semi_empirical_control_type USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_type USE cp_output_handling, ONLY: cp_print_key_finished_output,& cp_print_key_unit_nr USE cp_units, ONLY: cp_unit_from_cp2k USE external_potential_types, ONLY: all_potential_type,& get_potential,& gth_potential_type,& set_potential,& sgp_potential_type USE input_section_types, ONLY: section_vals_type,& section_vals_val_get USE kinds, ONLY: default_string_length,& dp USE paw_proj_set_types, ONLY: get_paw_proj_set,& paw_proj_set_type,& set_paw_proj_set USE qs_kind_types, ONLY: get_qs_kind,& get_qs_kind_set,& qs_kind_type USE string_utilities, ONLY: uppercase #include "./base/base_uses.f90" IMPLICIT NONE PRIVATE ! *** Global parameters (only in this module) CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_interactions' ! *** Public subroutines *** PUBLIC :: init_interaction_radii, init_se_nlradius, init_interaction_radii_orb_basis PUBLIC :: write_paw_radii, write_ppnl_radii, write_ppl_radii, write_core_charge_radii, & write_pgf_orb_radii CONTAINS ! ************************************************************************************************** !> \brief Initialize all the atomic kind radii for a given threshold value. !> \param qs_control ... !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \date 24.09.2002 !> \author GT !> \version 1.0 ! ************************************************************************************************** SUBROUTINE init_interaction_radii(qs_control, atomic_kind_set, qs_kind_set) TYPE(qs_control_type), INTENT(IN) :: qs_control TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set CHARACTER(len=*), PARAMETER :: routineN = 'init_interaction_radii', & routineP = moduleN//':'//routineN INTEGER :: handle, i, iexp_ppl, ikind, ip, iprj_ppnl, j, l, lppl, lppnl, lprj, lprj_ppnl, & maxl, maxnset, n_local, nexp_lpot, nexp_lsd, nexp_ppl, nkind, nloc INTEGER, DIMENSION(1:10) :: nrloc INTEGER, DIMENSION(:), POINTER :: nct_lpot, nct_lsd, nprj, nprj_ppnl LOGICAL :: ecp_local, llsd, lpot, paw_atom LOGICAL, DIMENSION(0:5) :: is_nonlocal REAL(dp), DIMENSION(1:10) :: aloc, bloc REAL(KIND=dp) :: alpha_core_charge, alpha_ppl, ccore_charge, cerf_ppl, core_charge_radius, & cval, ppl_radius, ppnl_radius, rcprj, zeta REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, a_nonlocal, alpha_lpot, & alpha_lsd, alpha_ppnl, c_local, & cexp_ppl REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj_ppnl, cval_lpot, cval_lsd, rzetprj, & zet REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nonlocal TYPE(all_potential_type), POINTER :: all_potential TYPE(gth_potential_type), POINTER :: gth_potential TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, aux_gw_basis, & harris_basis, lri_basis, mao_basis, orb_basis_set, ri_aux_basis_set, ri_basis, & ri_xas_basis, soft_basis TYPE(paw_proj_set_type), POINTER :: paw_proj_set TYPE(sgp_potential_type), POINTER :: sgp_potential CALL timeset(routineN, handle) NULLIFY (all_potential, gth_potential, sgp_potential, orb_basis_set, aux_basis_set, soft_basis, & lri_basis, mao_basis, paw_proj_set, harris_basis, ri_xas_basis) NULLIFY (nprj_ppnl, nprj) NULLIFY (alpha_ppnl, cexp_ppl, cprj_ppnl, zet) nkind = SIZE(atomic_kind_set) CALL get_qs_kind_set(qs_kind_set, maxnset=maxnset) nexp_ppl = 0 DO ikind = 1, nkind CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB") CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT") CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis, basis_type="LRI_AUX") CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_HXC") CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_aux_basis_set, basis_type="RI_AUX") CALL get_qs_kind(qs_kind_set(ikind), basis_set=mao_basis, basis_type="MAO") CALL get_qs_kind(qs_kind_set(ikind), basis_set=harris_basis, basis_type="HARRIS") CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_gw_basis, basis_type="AUX_GW") CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_xas_basis, basis_type="RI_XAS") CALL get_qs_kind(qs_kind_set(ikind), soft_basis_set=soft_basis) CALL get_qs_kind(qs_kind_set(ikind), & paw_proj_set=paw_proj_set, & paw_atom=paw_atom, & all_potential=all_potential, & gth_potential=gth_potential, & sgp_potential=sgp_potential) ! Calculate the orbital basis function radii *** ! For ALL electron this has to come before the calculation of the ! radii used to calculate the 3c lists. ! Indeed, there the radii of the GTO primitives are needed IF (ASSOCIATED(orb_basis_set)) THEN CALL init_interaction_radii_orb_basis(orb_basis_set, qs_control%eps_pgf_orb, & qs_control%eps_kg_orb) END IF IF (ASSOCIATED(all_potential)) THEN CALL get_potential(potential=all_potential, & alpha_core_charge=alpha_core_charge, & ccore_charge=ccore_charge) ! Calculate the radius of the core charge distribution core_charge_radius = exp_radius(0, alpha_core_charge, & qs_control%eps_core_charge, & ccore_charge) CALL set_potential(potential=all_potential, & core_charge_radius=core_charge_radius) ELSE IF (ASSOCIATED(gth_potential)) THEN CALL get_potential(potential=gth_potential, & alpha_core_charge=alpha_core_charge, & ccore_charge=ccore_charge, & alpha_ppl=alpha_ppl, & cerf_ppl=cerf_ppl, & nexp_ppl=nexp_ppl, & cexp_ppl=cexp_ppl, & lpot_present=lpot, & lsd_present=llsd, & lppnl=lppnl, & alpha_ppnl=alpha_ppnl, & nprj_ppnl=nprj_ppnl, & cprj_ppnl=cprj_ppnl) ! Calculate the radius of the core charge distribution core_charge_radius = exp_radius(0, alpha_core_charge, & qs_control%eps_core_charge, & ccore_charge) ! Calculate the radii of the local part ! of the Goedecker pseudopotential (GTH) ppl_radius = exp_radius(0, alpha_ppl, qs_control%eps_ppl, cerf_ppl) DO iexp_ppl = 1, nexp_ppl lppl = 2*(iexp_ppl - 1) ppl_radius = MAX(ppl_radius, & exp_radius(lppl, alpha_ppl, & qs_control%eps_ppl, & cexp_ppl(iexp_ppl))) END DO IF (lpot) THEN ! extended local potential CALL get_potential(potential=gth_potential, & nexp_lpot=nexp_lpot, & alpha_lpot=alpha_lpot, & nct_lpot=nct_lpot, & cval_lpot=cval_lpot) DO j = 1, nexp_lpot DO i = 1, nct_lpot(j) lppl = 2*(i - 1) ppl_radius = MAX(ppl_radius, & exp_radius(lppl, alpha_lpot(j), qs_control%eps_ppl, cval_lpot(i, j))) END DO END DO END IF IF (llsd) THEN ! lsd local potential CALL get_potential(potential=gth_potential, & nexp_lsd=nexp_lsd, & alpha_lsd=alpha_lsd, & nct_lsd=nct_lsd, & cval_lsd=cval_lsd) DO j = 1, nexp_lsd DO i = 1, nct_lsd(j) lppl = 2*(i - 1) ppl_radius = MAX(ppl_radius, & exp_radius(lppl, alpha_lsd(j), qs_control%eps_ppl, cval_lsd(i, j))) END DO END DO END IF ! Calculate the radii of the non-local part ! of the Goedecker pseudopotential (GTH) ppnl_radius = 0.0_dp DO l = 0, lppnl DO iprj_ppnl = 1, nprj_ppnl(l) lprj_ppnl = l + 2*(iprj_ppnl - 1) ppnl_radius = MAX(ppnl_radius, & exp_radius(lprj_ppnl, alpha_ppnl(l), & qs_control%eps_pgf_orb, & cprj_ppnl(iprj_ppnl, l))) END DO END DO CALL set_potential(potential=gth_potential, & core_charge_radius=core_charge_radius, & ppl_radius=ppl_radius, & ppnl_radius=ppnl_radius) ELSE IF (ASSOCIATED(sgp_potential)) THEN ! Calculate the radius of the core charge distribution CALL get_potential(potential=sgp_potential, & alpha_core_charge=alpha_core_charge, & ccore_charge=ccore_charge) core_charge_radius = exp_radius(0, alpha_core_charge, & qs_control%eps_core_charge, & ccore_charge) ! Calculate the radii of the local part ppl_radius = core_charge_radius CALL get_potential(potential=sgp_potential, ecp_local=ecp_local) IF (ecp_local) THEN CALL get_potential(potential=sgp_potential, nloc=nloc, nrloc=nrloc, aloc=aloc, bloc=bloc) DO i = 1, nloc lppl = MAX(0, nrloc(i) - 2) ppl_radius = MAX(ppl_radius, & exp_radius(lppl, bloc(i), qs_control%eps_ppl, aloc(i))) END DO ELSE CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local) DO i = 1, n_local ppl_radius = MAX(ppl_radius, & exp_radius(0, a_local(i), qs_control%eps_ppl, c_local(i))) END DO END IF ! Calculate the radii of the non-local part ppnl_radius = 0.0_dp CALL get_potential(potential=sgp_potential, lmax=lppnl, n_nonlocal=nloc) CALL get_potential(potential=sgp_potential, is_nonlocal=is_nonlocal, & a_nonlocal=a_nonlocal, c_nonlocal=c_nonlocal) DO l = 0, lppnl IF (is_nonlocal(l)) THEN DO i = 1, nloc cval = MAXVAL(ABS(c_nonlocal(i, :, l))) ppnl_radius = MAX(ppnl_radius, & exp_radius(l, a_nonlocal(i), qs_control%eps_pgf_orb, cval)) END DO END IF END DO CALL set_potential(potential=sgp_potential, & core_charge_radius=core_charge_radius, & ppl_radius=ppl_radius, & ppnl_radius=ppnl_radius) END IF ! Calculate the aux fit orbital basis function radii IF (ASSOCIATED(aux_fit_basis_set)) THEN CALL init_interaction_radii_orb_basis(aux_fit_basis_set, qs_control%eps_pgf_orb, & qs_control%eps_kg_orb) END IF ! Calculate the aux orbital basis function radii IF (ASSOCIATED(aux_basis_set)) THEN CALL init_interaction_radii_orb_basis(aux_basis_set, qs_control%eps_pgf_orb) END IF ! Calculate the RI aux orbital basis function radii IF (ASSOCIATED(RI_aux_basis_set)) THEN CALL init_interaction_radii_orb_basis(RI_aux_basis_set, qs_control%eps_pgf_orb) END IF ! Calculate the MAO orbital basis function radii IF (ASSOCIATED(mao_basis)) THEN CALL init_interaction_radii_orb_basis(mao_basis, qs_control%eps_pgf_orb) END IF ! Calculate the HARRIS orbital basis function radii IF (ASSOCIATED(harris_basis)) THEN CALL init_interaction_radii_orb_basis(harris_basis, qs_control%eps_pgf_orb) END IF ! Calculate the AUX_GW orbital basis function radii IF (ASSOCIATED(aux_gw_basis)) THEN CALL init_interaction_radii_orb_basis(aux_gw_basis, qs_control%eps_pgf_orb) END IF ! Calculate the soft orbital basis function radii IF (ASSOCIATED(soft_basis)) THEN IF (paw_atom) THEN CALL init_interaction_radii_orb_basis(soft_basis, qs_control%eps_pgf_orb) END IF END IF ! Calculate the lri basis function radii IF (ASSOCIATED(lri_basis)) THEN CALL init_interaction_radii_orb_basis(lri_basis, qs_control%eps_pgf_orb) END IF IF (ASSOCIATED(ri_basis)) THEN CALL init_interaction_radii_orb_basis(ri_basis, qs_control%eps_pgf_orb) END IF IF (ASSOCIATED(ri_xas_basis)) THEN CALL init_interaction_radii_orb_basis(ri_xas_basis, qs_control%eps_pgf_orb) END IF ! Calculate the paw projector basis function radii IF (ASSOCIATED(paw_proj_set)) THEN CALL get_paw_proj_set(paw_proj_set=paw_proj_set, & nprj=nprj, & maxl=maxl, & zetprj=zet, & rzetprj=rzetprj) rcprj = 0.0_dp rzetprj = 0.0_dp DO lprj = 0, maxl DO ip = 1, nprj(lprj) zeta = zet(ip, lprj) rzetprj(ip, lprj) = MAX(rzetprj(ip, lprj), & exp_radius(lprj, zeta, & qs_control%eps_pgf_orb, 1.0_dp)) rcprj = MAX(rcprj, rzetprj(ip, lprj)) ENDDO ENDDO CALL set_paw_proj_set(paw_proj_set=paw_proj_set, & rzetprj=rzetprj, & rcprj=rcprj) END IF END DO ! ikind CALL timestop(handle) END SUBROUTINE init_interaction_radii ! ************************************************************************************************** !> \brief ... !> \param orb_basis_set ... !> \param eps_pgf_orb ... !> \param eps_pgf_short ... ! ************************************************************************************************** SUBROUTINE init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short) TYPE(gto_basis_set_type), POINTER :: orb_basis_set REAL(KIND=dp), INTENT(IN) :: eps_pgf_orb REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_pgf_short CHARACTER(len=*), PARAMETER :: routineN = 'init_interaction_radii_orb_basis', & routineP = moduleN//':'//routineN INTEGER :: ipgf, iset, ishell, l, nset INTEGER, DIMENSION(:), POINTER :: npgf, nshell INTEGER, DIMENSION(:, :), POINTER :: lshell REAL(KIND=dp) :: eps_short, gcca, kind_radius, & short_radius, zeta REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius, zet REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc IF (ASSOCIATED(orb_basis_set)) THEN IF (PRESENT(eps_pgf_short)) THEN eps_short = eps_pgf_short ELSE eps_short = eps_pgf_orb END IF CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & nset=nset, & nshell=nshell, & npgf=npgf, & l=lshell, & zet=zet, & gcc=gcc, & pgf_radius=pgf_radius, & set_radius=set_radius) kind_radius = 0.0_dp short_radius = 0.0_dp DO iset = 1, nset set_radius(iset) = 0.0_dp DO ipgf = 1, npgf(iset) pgf_radius(ipgf, iset) = 0.0_dp DO ishell = 1, nshell(iset) l = lshell(ishell, iset) gcca = gcc(ipgf, ishell, iset) zeta = zet(ipgf, iset) pgf_radius(ipgf, iset) = MAX(pgf_radius(ipgf, iset), & exp_radius(l, zeta, & eps_pgf_orb, & gcca)) short_radius = MAX(short_radius, & exp_radius(l, zeta, & eps_short, & gcca)) END DO set_radius(iset) = MAX(set_radius(iset), pgf_radius(ipgf, iset)) END DO kind_radius = MAX(kind_radius, set_radius(iset)) END DO CALL set_gto_basis_set(gto_basis_set=orb_basis_set, & pgf_radius=pgf_radius, & set_radius=set_radius, & kind_radius=kind_radius, & short_kind_radius=short_radius) END IF END SUBROUTINE init_interaction_radii_orb_basis ! ************************************************************************************************** !> \brief ... !> \param se_control ... !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \param subsys_section ... ! ************************************************************************************************** SUBROUTINE init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, & subsys_section) TYPE(semi_empirical_control_type), POINTER :: se_control TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: subsys_section CHARACTER(len=*), PARAMETER :: routineN = 'init_se_nlradius', & routineP = moduleN//':'//routineN INTEGER :: ikind, nkind REAL(KIND=dp) :: kind_radius TYPE(gto_basis_set_type), POINTER :: orb_basis_set TYPE(qs_kind_type), POINTER :: qs_kind NULLIFY (orb_basis_set, qs_kind) nkind = SIZE(qs_kind_set) DO ikind = 1, nkind qs_kind => qs_kind_set(ikind) CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set) IF (ASSOCIATED(orb_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & kind_radius=kind_radius) kind_radius = MAX(kind_radius, se_control%cutoff_exc) CALL set_gto_basis_set(gto_basis_set=orb_basis_set, & kind_radius=kind_radius) END IF END DO CALL write_kind_radii(atomic_kind_set, qs_kind_set, subsys_section) END SUBROUTINE init_se_nlradius ! ************************************************************************************************** !> \brief ... !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \param subsys_section ... ! ************************************************************************************************** SUBROUTINE write_kind_radii(atomic_kind_set, qs_kind_set, subsys_section) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: subsys_section CHARACTER(LEN=10) :: bas CHARACTER(LEN=default_string_length) :: name, unit_str INTEGER :: ikind, nkind, output_unit REAL(KIND=dp) :: conv, kind_radius TYPE(cp_logger_type), POINTER :: logger TYPE(gto_basis_set_type), POINTER :: orb_basis_set NULLIFY (logger) logger => cp_get_default_logger() NULLIFY (orb_basis_set) bas = "ORBITAL " nkind = SIZE(atomic_kind_set) ! *** Print the kind radii *** output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/KIND_RADII", extension=".Log") CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) IF (output_unit > 0) THEN WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & "RADII: "//TRIM(bas)//" BASIS in "//TRIM(unit_str), & "Kind", "Label", "Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) IF (ASSOCIATED(orb_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & kind_radius=kind_radius) WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & ikind, name, kind_radius*conv ELSE WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T72,A9)") & ikind, name, "no basis" END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/KIND_RADII") END SUBROUTINE write_kind_radii ! ************************************************************************************************** !> \brief Write the radii of the core charge distributions to the output !> unit. !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \param subsys_section ... !> \date 15.09.2000 !> \author MK !> \version 1.0 ! ************************************************************************************************** SUBROUTINE write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: subsys_section CHARACTER(LEN=default_string_length) :: name, unit_str INTEGER :: ikind, nkind, output_unit REAL(KIND=dp) :: conv, core_charge_radius TYPE(all_potential_type), POINTER :: all_potential TYPE(cp_logger_type), POINTER :: logger TYPE(gth_potential_type), POINTER :: gth_potential NULLIFY (logger) logger => cp_get_default_logger() output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/CORE_CHARGE_RADII", extension=".Log") CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) IF (output_unit > 0) THEN nkind = SIZE(atomic_kind_set) WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & "RADII: CORE CHARGE DISTRIBUTIONS in "// & TRIM(unit_str), "Kind", "Label", "Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) CALL get_qs_kind(qs_kind_set(ikind), & all_potential=all_potential, gth_potential=gth_potential) IF (ASSOCIATED(all_potential)) THEN CALL get_potential(potential=all_potential, & core_charge_radius=core_charge_radius) WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & ikind, name, core_charge_radius*conv ELSE IF (ASSOCIATED(gth_potential)) THEN CALL get_potential(potential=gth_potential, & core_charge_radius=core_charge_radius) WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & ikind, name, core_charge_radius*conv END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/CORE_CHARGE_RADII") END SUBROUTINE write_core_charge_radii ! ************************************************************************************************** !> \brief Write the orbital basis function radii to the output unit. !> \param basis ... !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \param subsys_section ... !> \date 05.06.2000 !> \author MK !> \version 1.0 ! ************************************************************************************************** SUBROUTINE write_pgf_orb_radii(basis, atomic_kind_set, qs_kind_set, subsys_section) CHARACTER(len=*) :: basis TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: subsys_section CHARACTER(len=*), PARAMETER :: routineN = 'write_pgf_orb_radii', & routineP = moduleN//':'//routineN CHARACTER(LEN=10) :: bas CHARACTER(LEN=default_string_length) :: name, unit_str INTEGER :: ikind, ipgf, iset, nkind, nset, & output_unit INTEGER, DIMENSION(:), POINTER :: npgf REAL(KIND=dp) :: conv, kind_radius, short_kind_radius REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius TYPE(cp_logger_type), POINTER :: logger TYPE(gto_basis_set_type), POINTER :: aux_basis_set, lri_basis_set, & orb_basis_set NULLIFY (logger) logger => cp_get_default_logger() NULLIFY (aux_basis_set, orb_basis_set, lri_basis_set) bas = " " bas(1:3) = basis(1:3) CALL uppercase(bas) IF (bas(1:3) == "ORB") THEN bas = "ORBITAL " ELSE IF (bas(1:3) == "AUX") THEN bas = "AUXILLIARY" ELSE IF (bas(1:3) == "LRI") THEN bas = "LOCAL RI" ELSE CPABORT("Undefined basis set type") END IF nkind = SIZE(qs_kind_set) ! *** Print the kind radii *** output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/KIND_RADII", extension=".Log") CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) IF (output_unit > 0) THEN WRITE (UNIT=output_unit, FMT="(/,T2,A,T46,A,T53,A,T63,A,T71,A)") & "RADII: "//TRIM(bas)//" BASIS in "//TRIM(unit_str), & "Kind", "Label", "Radius", "OCE Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) IF (bas(1:3) == "ORB") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) IF (ASSOCIATED(orb_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & kind_radius=kind_radius, & short_kind_radius=short_kind_radius) END IF ELSE IF (bas(1:3) == "AUX") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") IF (ASSOCIATED(aux_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=aux_basis_set, & kind_radius=kind_radius) END IF ELSE IF (bas(1:3) == "LOC") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX") IF (ASSOCIATED(lri_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=lri_basis_set, & kind_radius=kind_radius) END IF ELSE CPABORT("Undefined basis set type") END IF IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN WRITE (UNIT=output_unit, FMT="(T45,I5,T53,A5,T57,F12.6,T69,F12.6)") & ikind, name, kind_radius*conv, short_kind_radius*conv ELSE WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T72,A9)") & ikind, name, "no basis" END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/KIND_RADII") ! *** Print the shell set radii *** output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/SET_RADII", extension=".Log") IF (output_unit > 0) THEN WRITE (UNIT=output_unit, FMT="(/,T2,A,T51,A,T57,A,T65,A,T75,A)") & "RADII: SHELL SETS OF "//TRIM(bas)//" BASIS in "// & TRIM(unit_str), "Kind", "Label", "Set", "Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) IF (bas(1:3) == "ORB") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) IF (ASSOCIATED(orb_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & nset=nset, & set_radius=set_radius) END IF ELSE IF (bas(1:3) == "AUX") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") IF (ASSOCIATED(aux_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=aux_basis_set, & nset=nset, & set_radius=set_radius) END IF ELSE IF (bas(1:3) == "LOC") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX") IF (ASSOCIATED(lri_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=lri_basis_set, & nset=nset, & set_radius=set_radius) END IF ELSE CPABORT("Undefined basis set type") END IF IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,(T63,I5,T69,F12.6))") & ikind, name, (iset, set_radius(iset)*conv, iset=1, nset) ELSE WRITE (UNIT=output_unit, FMT="(T50,I5,T58,A5,T73,A8)") & ikind, name, "no basis" END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/SET_RADII") ! *** Print the primitive Gaussian function radii *** output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/PGF_RADII", extension=".Log") IF (output_unit > 0) THEN WRITE (UNIT=output_unit, FMT="(/,T2,A,T51,A,T57,A,T65,A,T75,A)") & "RADII: PRIMITIVE GAUSSIANS OF "//TRIM(bas)//" BASIS in "// & TRIM(unit_str), "Kind", "Label", "Set", "Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) IF (bas(1:3) == "ORB") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) IF (ASSOCIATED(orb_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & nset=nset, & npgf=npgf, & pgf_radius=pgf_radius) END IF ELSE IF (bas(1:3) == "AUX") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX") IF (ASSOCIATED(aux_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=aux_basis_set, & nset=nset, & npgf=npgf, & pgf_radius=pgf_radius) END IF ELSE IF (bas(1:3) == "LOC") THEN CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type="LRI_AUX") IF (ASSOCIATED(lri_basis_set)) THEN CALL get_gto_basis_set(gto_basis_set=lri_basis_set, & nset=nset, & npgf=npgf, & pgf_radius=pgf_radius) END IF ELSE CPABORT("Undefined basis type") END IF IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set) .OR. & ASSOCIATED(lri_basis_set)) THEN DO iset = 1, nset WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,T63,I5,(T69,F12.6))") & ikind, name, iset, & (pgf_radius(ipgf, iset)*conv, ipgf=1, npgf(iset)) END DO ELSE WRITE (UNIT=output_unit, FMT="(T50,I5,T58,A5,T73,A8)") & ikind, name, "no basis" END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/PGF_RADII") END SUBROUTINE write_pgf_orb_radii ! ************************************************************************************************** !> \brief Write the radii of the exponential functions of the Goedecker !> pseudopotential (GTH, local part) to the logical unit number !> "output_unit". !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \param subsys_section ... !> \date 06.10.2000 !> \author MK !> \version 1.0 ! ************************************************************************************************** SUBROUTINE write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: subsys_section CHARACTER(LEN=default_string_length) :: name, unit_str INTEGER :: ikind, nkind, output_unit REAL(KIND=dp) :: conv, ppl_radius TYPE(cp_logger_type), POINTER :: logger TYPE(gth_potential_type), POINTER :: gth_potential NULLIFY (logger) logger => cp_get_default_logger() output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/PPL_RADII", extension=".Log") CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) IF (output_unit > 0) THEN nkind = SIZE(atomic_kind_set) WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & "RADII: LOCAL PART OF GTH/ELP PP in "// & TRIM(unit_str), "Kind", "Label", "Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential) IF (ASSOCIATED(gth_potential)) THEN CALL get_potential(potential=gth_potential, & ppl_radius=ppl_radius) WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & ikind, name, ppl_radius*conv END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/PPL_RADII") END SUBROUTINE write_ppl_radii ! ************************************************************************************************** !> \brief Write the radii of the projector functions of the Goedecker !> pseudopotential (GTH, non-local part) to the logical unit number !> "output_unit". !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \param subsys_section ... !> \date 06.10.2000 !> \author MK !> \version 1.0 ! ************************************************************************************************** SUBROUTINE write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: subsys_section CHARACTER(LEN=default_string_length) :: name, unit_str INTEGER :: ikind, nkind, output_unit REAL(KIND=dp) :: conv, ppnl_radius TYPE(cp_logger_type), POINTER :: logger TYPE(gth_potential_type), POINTER :: gth_potential NULLIFY (logger) logger => cp_get_default_logger() output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/PPNL_RADII", extension=".Log") CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) IF (output_unit > 0) THEN nkind = SIZE(atomic_kind_set) WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & "RADII: NON-LOCAL PART OF GTH PP in "// & TRIM(unit_str), "Kind", "Label", "Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential) IF (ASSOCIATED(gth_potential)) THEN CALL get_potential(potential=gth_potential, & ppnl_radius=ppnl_radius) WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & ikind, name, ppnl_radius*conv END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/PPNL_RADII") END SUBROUTINE write_ppnl_radii ! ************************************************************************************************** !> \brief Write the radii of the one center projector !> \param atomic_kind_set ... !> \param qs_kind_set ... !> \param subsys_section ... !> \version 1.0 ! ************************************************************************************************** SUBROUTINE write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section) TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(section_vals_type), POINTER :: subsys_section CHARACTER(LEN=default_string_length) :: name, unit_str INTEGER :: ikind, nkind, output_unit LOGICAL :: paw_atom REAL(KIND=dp) :: conv, rcprj TYPE(cp_logger_type), POINTER :: logger TYPE(paw_proj_set_type), POINTER :: paw_proj_set NULLIFY (logger) logger => cp_get_default_logger() output_unit = cp_print_key_unit_nr(logger, subsys_section, & "PRINT%RADII/GAPW_PRJ_RADII", extension=".Log") CALL section_vals_val_get(subsys_section, "PRINT%RADII%UNIT", c_val=unit_str) conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) IF (output_unit > 0) THEN nkind = SIZE(qs_kind_set) WRITE (UNIT=output_unit, FMT="(/,T2,A,T56,A,T63,A,T75,A)") & "RADII: ONE CENTER PROJECTORS in "// & TRIM(unit_str), "Kind", "Label", "Radius" DO ikind = 1, nkind CALL get_atomic_kind(atomic_kind_set(ikind), name=name) CALL get_qs_kind(qs_kind_set(ikind), & paw_atom=paw_atom, paw_proj_set=paw_proj_set) IF (paw_atom .AND. ASSOCIATED(paw_proj_set)) THEN CALL get_paw_proj_set(paw_proj_set=paw_proj_set, & rcprj=rcprj) WRITE (UNIT=output_unit, FMT="(T55,I5,T63,A5,T69,F12.6)") & ikind, name, rcprj*conv END IF END DO END IF CALL cp_print_key_finished_output(output_unit, logger, subsys_section, & "PRINT%RADII/GAPW_PRJ_RADII") END SUBROUTINE write_paw_radii END MODULE qs_interactions