1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE qs_gapw_densities 8 USE atomic_kind_types, ONLY: atomic_kind_type,& 9 get_atomic_kind 10 USE cp_control_types, ONLY: dft_control_type,& 11 gapw_control_type 12 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit 13 USE cp_para_types, ONLY: cp_para_env_type 14 USE kinds, ONLY: dp 15 USE message_passing, ONLY: mp_sum 16 USE qs_charges_types, ONLY: qs_charges_type 17 USE qs_environment_types, ONLY: get_qs_env,& 18 qs_environment_type 19 USE qs_kind_types, ONLY: get_qs_kind,& 20 qs_kind_type 21 USE qs_local_rho_types, ONLY: local_rho_type 22 USE qs_rho0_ggrid, ONLY: put_rho0_on_grid 23 USE qs_rho0_methods, ONLY: calculate_rho0_atom 24 USE qs_rho0_types, ONLY: rho0_atom_type,& 25 rho0_mpole_type 26 USE qs_rho_atom_methods, ONLY: calculate_rho_atom 27 USE qs_rho_atom_types, ONLY: rho_atom_type 28#include "./base/base_uses.f90" 29 30 IMPLICIT NONE 31 32 PRIVATE 33 34 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gapw_densities' 35 36 PUBLIC :: prepare_gapw_den 37 38CONTAINS 39 40! ************************************************************************************************** 41!> \brief ... 42!> \param qs_env ... 43!> \param local_rho_set ... 44!> \param do_rho0 ... 45!> \param kind_set_external can be provided to use different projectors/grids/basis than the default 46! ************************************************************************************************** 47 SUBROUTINE prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external) 48 49 TYPE(qs_environment_type), POINTER :: qs_env 50 TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set 51 LOGICAL, INTENT(IN), OPTIONAL :: do_rho0 52 TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, & 53 POINTER :: kind_set_external 54 55 CHARACTER(len=*), PARAMETER :: routineN = 'prepare_gapw_den', & 56 routineP = moduleN//':'//routineN 57 58 INTEGER :: handle, ikind, ispin, natom, nspins, & 59 output_unit 60 INTEGER, DIMENSION(:), POINTER :: atom_list 61 LOGICAL :: extern, my_do_rho0, paw_atom 62 REAL(dp) :: rho0_h_tot, tot_rs_int 63 REAL(dp), DIMENSION(:), POINTER :: rho1_h_tot, rho1_s_tot 64 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 65 TYPE(cp_para_env_type), POINTER :: para_env 66 TYPE(dft_control_type), POINTER :: dft_control 67 TYPE(gapw_control_type), POINTER :: gapw_control 68 TYPE(qs_charges_type), POINTER :: qs_charges 69 TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set 70 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set 71 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 72 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set 73 74 CALL timeset(routineN, handle) 75 76 NULLIFY (atomic_kind_set) 77 NULLIFY (my_kind_set) 78 NULLIFY (dft_control) 79 NULLIFY (gapw_control) 80 NULLIFY (para_env) 81 NULLIFY (atom_list) 82 NULLIFY (rho0_mpole) 83 NULLIFY (qs_charges) 84 NULLIFY (rho1_h_tot, rho1_s_tot) 85 NULLIFY (rho_atom_set) 86 NULLIFY (rho0_atom_set) 87 88 my_do_rho0 = .TRUE. 89 IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0 90 91 output_unit = cp_logger_get_default_io_unit() 92 93 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 94 para_env=para_env, & 95 qs_charges=qs_charges, & 96 qs_kind_set=my_kind_set, & 97 atomic_kind_set=atomic_kind_set, & 98 rho0_mpole=rho0_mpole, & 99 rho_atom_set=rho_atom_set, & 100 rho0_atom_set=rho0_atom_set) 101 102 gapw_control => dft_control%qs_control%gapw_control 103 104 IF (PRESENT(local_rho_set)) THEN 105 rho_atom_set => local_rho_set%rho_atom_set 106 IF (my_do_rho0) THEN 107 rho0_mpole => local_rho_set%rho0_mpole 108 rho0_atom_set => local_rho_set%rho0_atom_set 109 END IF 110 END IF 111 112 extern = .FALSE. 113 IF (PRESENT(kind_set_external)) THEN 114 CPASSERT(ASSOCIATED(kind_set_external)) 115 my_kind_set => kind_set_external 116 extern = .TRUE. 117 END IF 118 119 nspins = dft_control%nspins 120 121 rho0_h_tot = 0.0_dp 122 ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins)) 123 rho1_h_tot = 0.0_dp 124 rho1_s_tot = 0.0_dp 125 126 DO ikind = 1, SIZE(atomic_kind_set) 127 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom) 128 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom) 129 130! Calculate rho1_h and rho1_s on the radial grids centered on the atomic position 131 IF (paw_atom) & 132 CALL calculate_rho_atom(para_env, rho_atom_set, my_kind_set(ikind), & 133 atom_list, natom, nspins, rho1_h_tot, rho1_s_tot) 134 135! Calculate rho0_h and rho0_s on the radial grids centered on the atomic position 136 IF (my_do_rho0) & 137 CALL calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, rho0_mpole, & 138 atom_list, natom, ikind, my_kind_set(ikind), rho0_h_tot) 139 140 ENDDO 141 142 !Do not mess with charges if using a non-default kind_set 143 IF (.NOT. extern) THEN 144 CALL mp_sum(rho1_h_tot, para_env%group) 145 CALL mp_sum(rho1_s_tot, para_env%group) 146 DO ispin = 1, nspins 147 qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin) 148 qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin) 149 END DO 150 151 IF (my_do_rho0) THEN 152 rho0_mpole%total_rho0_h = -rho0_h_tot 153! Put the rho0_soft on the global grid 154 CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int) 155 IF (ABS(rho0_h_tot) .GE. 1.0E-5_dp) THEN 156 IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) .GT. 1.0E-3_dp) THEN 157 IF (output_unit > 0) THEN 158 WRITE (output_unit, '(/,72("*"))') 159 WRITE (output_unit, '(T2,A,T66,1E20.8)') & 160 "WARNING: rho0 calculated on the local grid is :", -rho0_h_tot, & 161 " rho0 calculated on the global grid is :", tot_rs_int 162 WRITE (output_unit, '(T2,A)') & 163 " bad integration" 164 WRITE (output_unit, '(72("*"),/)') 165 END IF 166 END IF 167 END IF 168 qs_charges%total_rho0_soft_rspace = tot_rs_int 169 qs_charges%total_rho0_hard_lebedev = rho0_h_tot 170 ELSE 171 qs_charges%total_rho0_hard_lebedev = 0.0_dp 172 END IF 173 END IF 174 175 DEALLOCATE (rho1_h_tot, rho1_s_tot) 176 177 CALL timestop(handle) 178 179 END SUBROUTINE prepare_gapw_den 180 181END MODULE qs_gapw_densities 182