Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : MODULE qs_gapw_densities
10 : USE atomic_kind_types, ONLY: atomic_kind_type,&
11 : get_atomic_kind
12 : USE cp_control_types, ONLY: dft_control_type,&
13 : gapw_control_type
14 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
15 : USE kinds, ONLY: dp
16 : USE message_passing, ONLY: mp_para_env_type
17 : USE pw_env_types, ONLY: pw_env_get,&
18 : pw_env_type
19 : USE pw_pool_types, ONLY: pw_pool_p_type
20 : USE qs_charges_types, ONLY: qs_charges_type
21 : USE qs_cneo_ggrid, ONLY: put_rhoz_cneo_s_on_grid
22 : USE qs_cneo_types, ONLY: rhoz_cneo_type
23 : USE qs_environment_types, ONLY: get_qs_env,&
24 : qs_environment_type
25 : USE qs_kind_types, ONLY: get_qs_kind,&
26 : qs_kind_type
27 : USE qs_local_rho_types, ONLY: local_rho_type
28 : USE qs_rho0_ggrid, ONLY: put_rho0_on_grid
29 : USE qs_rho0_methods, ONLY: calculate_rho0_atom
30 : USE qs_rho0_types, ONLY: rho0_atom_type,&
31 : rho0_mpole_type
32 : USE qs_rho_atom_methods, ONLY: calculate_rho_atom
33 : USE qs_rho_atom_types, ONLY: rho_atom_type
34 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
35 : realspace_grid_type
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gapw_densities'
43 :
44 : PUBLIC :: prepare_gapw_den
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief ...
50 : !> \param qs_env ...
51 : !> \param local_rho_set ...
52 : !> \param do_rho0 ...
53 : !> \param kind_set_external can be provided to use different projectors/grids/basis than the default
54 : !> \param pw_env_sub ...
55 : ! **************************************************************************************************
56 41028 : SUBROUTINE prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external, pw_env_sub)
57 :
58 : TYPE(qs_environment_type), POINTER :: qs_env
59 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
60 : LOGICAL, INTENT(IN), OPTIONAL :: do_rho0
61 : TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
62 : POINTER :: kind_set_external
63 : TYPE(pw_env_type), OPTIONAL :: pw_env_sub
64 :
65 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_gapw_den'
66 :
67 : INTEGER :: handle, ikind, ispin, natom, nspins, &
68 : output_unit
69 41028 : INTEGER, DIMENSION(:), POINTER :: atom_list
70 : LOGICAL :: extern, my_do_rho0, paw_atom
71 : REAL(dp) :: rho0_h_tot, rho1_h_aspin, rho1_h_spin, &
72 : rho1_s_aspin, rho1_s_spin, tot_rs_int
73 41028 : REAL(dp), DIMENSION(:), POINTER :: rho1_h_tot, rho1_s_tot
74 41028 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
75 : TYPE(dft_control_type), POINTER :: dft_control
76 : TYPE(gapw_control_type), POINTER :: gapw_control
77 : TYPE(mp_para_env_type), POINTER :: para_env
78 41028 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: my_pools
79 : TYPE(qs_charges_type), POINTER :: qs_charges
80 41028 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
81 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
82 41028 : POINTER :: my_rs_descs
83 41028 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: my_rs_grids
84 41028 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
85 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
86 41028 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
87 41028 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
88 :
89 41028 : CALL timeset(routineN, handle)
90 :
91 41028 : NULLIFY (atomic_kind_set)
92 41028 : NULLIFY (my_kind_set)
93 41028 : NULLIFY (dft_control)
94 41028 : NULLIFY (gapw_control)
95 41028 : NULLIFY (para_env)
96 41028 : NULLIFY (atom_list)
97 41028 : NULLIFY (rho0_mpole)
98 41028 : NULLIFY (qs_charges)
99 : NULLIFY (rho1_h_tot, rho1_s_tot)
100 41028 : NULLIFY (rho_atom_set)
101 41028 : NULLIFY (rho0_atom_set)
102 :
103 41028 : my_do_rho0 = .TRUE.
104 41028 : IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0
105 :
106 41028 : output_unit = cp_logger_get_default_io_unit()
107 :
108 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
109 : para_env=para_env, &
110 : qs_charges=qs_charges, &
111 : qs_kind_set=my_kind_set, &
112 : atomic_kind_set=atomic_kind_set, &
113 : rho0_mpole=rho0_mpole, &
114 : rho_atom_set=rho_atom_set, &
115 : rho0_atom_set=rho0_atom_set, &
116 41028 : rhoz_cneo_set=rhoz_cneo_set)
117 :
118 41028 : gapw_control => dft_control%qs_control%gapw_control
119 :
120 : ! If TDDFPT%MGRID is defined, overwrite QS grid info accordingly
121 41028 : IF (PRESENT(local_rho_set)) THEN
122 13286 : rho_atom_set => local_rho_set%rho_atom_set
123 13286 : rhoz_cneo_set => local_rho_set%rhoz_cneo_set
124 13286 : IF (my_do_rho0) THEN
125 5380 : rho0_mpole => local_rho_set%rho0_mpole
126 5380 : rho0_atom_set => local_rho_set%rho0_atom_set
127 : END IF
128 : END IF
129 :
130 41028 : extern = .FALSE.
131 41028 : IF (PRESENT(kind_set_external)) THEN
132 5568 : CPASSERT(ASSOCIATED(kind_set_external))
133 5568 : my_kind_set => kind_set_external
134 5568 : extern = .TRUE.
135 : END IF
136 :
137 41028 : nspins = dft_control%nspins
138 :
139 41028 : rho0_h_tot = 0.0_dp
140 164112 : ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins))
141 87394 : rho1_h_tot = 0.0_dp
142 87394 : rho1_s_tot = 0.0_dp
143 :
144 41028 : rho1_h_spin = 0.0_dp
145 41028 : rho1_s_spin = 0.0_dp
146 41028 : rho1_h_aspin = 0.0_dp
147 41028 : rho1_s_aspin = 0.0_dp
148 :
149 123730 : DO ikind = 1, SIZE(atomic_kind_set)
150 82702 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
151 82702 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
152 :
153 : !Calculate rho1_h and rho1_s on the radial grids centered on the atomic position
154 82702 : IF (paw_atom) THEN
155 : CALL calculate_rho_atom(para_env, rho_atom_set, my_kind_set(ikind), &
156 : atom_list, natom, nspins, rho1_h_tot, rho1_s_tot, &
157 75718 : rho1_h_spin, rho1_s_spin, rho1_h_aspin, rho1_s_aspin)
158 : END IF
159 :
160 : !Calculate rho0_h and rho0_s on the radial grids centered on the atomic position
161 206432 : IF (my_do_rho0) THEN
162 : CALL calculate_rho0_atom(gapw_control, rho_atom_set, rhoz_cneo_set, rho0_atom_set, &
163 : rho0_mpole, atom_list, natom, ikind, my_kind_set(ikind), &
164 56050 : rho0_h_tot)
165 : END IF
166 : END DO
167 :
168 : !Do not mess with charges if using a non-default kind_set
169 41028 : IF (.NOT. extern) THEN
170 115208 : CALL para_env%sum(rho1_h_tot)
171 115208 : CALL para_env%sum(rho1_s_tot)
172 75334 : DO ispin = 1, nspins
173 39874 : qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin)
174 75334 : qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin)
175 : END DO
176 : !spin
177 35460 : CALL para_env%sum(rho1_h_spin)
178 35460 : CALL para_env%sum(rho1_s_spin)
179 35460 : CALL para_env%sum(rho1_h_aspin)
180 35460 : CALL para_env%sum(rho1_s_aspin)
181 35460 : qs_charges%total_rho_hard_spin = rho1_h_spin
182 35460 : qs_charges%total_rho_soft_spin = rho1_s_spin
183 35460 : qs_charges%total_rho_hard_abs_spin = rho1_h_aspin
184 35460 : qs_charges%total_rho_soft_abs_spin = rho1_s_aspin
185 :
186 35460 : IF (my_do_rho0) THEN
187 28226 : rho0_mpole%total_rho0_h = -rho0_h_tot
188 :
189 : ! When MGRID is defined within TDDFPT
190 28226 : IF (PRESENT(pw_env_sub)) THEN
191 : ! Find pool
192 2282 : NULLIFY (my_pools, my_rs_grids, my_rs_descs)
193 : CALL pw_env_get(pw_env=pw_env_sub, rs_grids=my_rs_grids, &
194 2282 : rs_descs=my_rs_descs, pw_pools=my_pools)
195 : ! Put the rho0_soft on the global grid
196 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int, my_pools=my_pools, &
197 2282 : my_rs_grids=my_rs_grids, my_rs_descs=my_rs_descs)
198 : ELSE
199 : ! Put the rho0_soft on the global grid
200 25944 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int)
201 : END IF
202 :
203 28226 : IF (ABS(rho0_h_tot) >= 1.0E-5_dp) THEN
204 26088 : IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) > 1.0E-3_dp) THEN
205 1930 : IF (output_unit > 0) THEN
206 965 : WRITE (output_unit, '(/,72("*"))')
207 : WRITE (output_unit, '(T2,A,T66,1E20.8)') &
208 965 : "WARNING: rho0 calculated on the local grid is :", -rho0_h_tot, &
209 1930 : " rho0 calculated on the global grid is :", tot_rs_int
210 : WRITE (output_unit, '(T2,A)') &
211 965 : " bad integration"
212 965 : WRITE (output_unit, '(72("*"),/)')
213 : END IF
214 : END IF
215 : END IF
216 28226 : qs_charges%total_rho0_soft_rspace = tot_rs_int
217 28226 : qs_charges%total_rho0_hard_lebedev = rho0_h_tot
218 28226 : IF (rho0_mpole%do_cneo) THEN
219 : ! put soft tails of quantum nuclear charge densities on the global grid
220 48 : CALL put_rhoz_cneo_s_on_grid(qs_env, rho0_mpole, rhoz_cneo_set, tot_rs_int)
221 48 : IF (ABS(rho0_mpole%tot_rhoz_cneo_s) >= 1.0E-5_dp) THEN
222 40 : IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_mpole%tot_rhoz_cneo_s)) > 1.0E-3_dp) THEN
223 0 : IF (output_unit > 0) THEN
224 0 : WRITE (output_unit, '(/,72("*"))')
225 : WRITE (output_unit, '(T2,A,T66,1E20.8)') &
226 0 : "WARNING: rhoz_cneo_s calculated on the local grid is :", &
227 0 : rho0_mpole%tot_rhoz_cneo_s, &
228 0 : " rhoz_cneo_s calculated on the global grid is :", tot_rs_int
229 : WRITE (output_unit, '(T2,A)') &
230 0 : " bad integration"
231 0 : WRITE (output_unit, '(72("*"),/)')
232 : END IF
233 : END IF
234 : END IF
235 48 : qs_charges%total_rho1_soft_nuc_rspace = tot_rs_int
236 48 : qs_charges%total_rho1_soft_nuc_lebedev = rho0_mpole%tot_rhoz_cneo_s
237 : ELSE
238 28178 : qs_charges%total_rho1_soft_nuc_rspace = 0.0_dp
239 : END IF
240 : ELSE
241 7234 : qs_charges%total_rho0_hard_lebedev = 0.0_dp
242 : END IF
243 : END IF
244 :
245 41028 : DEALLOCATE (rho1_h_tot, rho1_s_tot)
246 :
247 41028 : CALL timestop(handle)
248 :
249 41028 : END SUBROUTINE prepare_gapw_den
250 :
251 : END MODULE qs_gapw_densities
|
