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 39484 : 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 39484 : 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 39484 : REAL(dp), DIMENSION(:), POINTER :: rho1_h_tot, rho1_s_tot
74 39484 : 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 39484 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: my_pools
79 : TYPE(qs_charges_type), POINTER :: qs_charges
80 39484 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
81 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
82 39484 : POINTER :: my_rs_descs
83 39484 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: my_rs_grids
84 39484 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
85 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
86 39484 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
87 39484 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
88 :
89 39484 : CALL timeset(routineN, handle)
90 :
91 39484 : NULLIFY (atomic_kind_set)
92 39484 : NULLIFY (my_kind_set)
93 39484 : NULLIFY (dft_control)
94 39484 : NULLIFY (gapw_control)
95 39484 : NULLIFY (para_env)
96 39484 : NULLIFY (atom_list)
97 39484 : NULLIFY (rho0_mpole)
98 39484 : NULLIFY (qs_charges)
99 : NULLIFY (rho1_h_tot, rho1_s_tot)
100 39484 : NULLIFY (rho_atom_set)
101 39484 : NULLIFY (rho0_atom_set)
102 :
103 39484 : my_do_rho0 = .TRUE.
104 39484 : IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0
105 :
106 39484 : 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 39484 : rhoz_cneo_set=rhoz_cneo_set)
117 :
118 39484 : gapw_control => dft_control%qs_control%gapw_control
119 :
120 : ! If TDDFPT%MGRID is defined, overwrite QS grid info accordingly
121 39484 : IF (PRESENT(local_rho_set)) THEN
122 12722 : rho_atom_set => local_rho_set%rho_atom_set
123 12722 : rhoz_cneo_set => local_rho_set%rhoz_cneo_set
124 12722 : IF (my_do_rho0) THEN
125 5270 : rho0_mpole => local_rho_set%rho0_mpole
126 5270 : rho0_atom_set => local_rho_set%rho0_atom_set
127 : END IF
128 : END IF
129 :
130 39484 : extern = .FALSE.
131 39484 : IF (PRESENT(kind_set_external)) THEN
132 5178 : CPASSERT(ASSOCIATED(kind_set_external))
133 5178 : my_kind_set => kind_set_external
134 5178 : extern = .TRUE.
135 : END IF
136 :
137 39484 : nspins = dft_control%nspins
138 :
139 39484 : rho0_h_tot = 0.0_dp
140 157936 : ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins))
141 84094 : rho1_h_tot = 0.0_dp
142 84094 : rho1_s_tot = 0.0_dp
143 :
144 39484 : rho1_h_spin = 0.0_dp
145 39484 : rho1_s_spin = 0.0_dp
146 39484 : rho1_h_aspin = 0.0_dp
147 39484 : rho1_s_aspin = 0.0_dp
148 :
149 118714 : DO ikind = 1, SIZE(atomic_kind_set)
150 79230 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
151 79230 : 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 79230 : 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 72754 : 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 79230 : IF (my_do_rho0) &
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 173160 : rho0_h_tot)
165 : END DO
166 :
167 : !Do not mess with charges if using a non-default kind_set
168 39484 : IF (.NOT. extern) THEN
169 111506 : CALL para_env%sum(rho1_h_tot)
170 111506 : CALL para_env%sum(rho1_s_tot)
171 72906 : DO ispin = 1, nspins
172 38600 : qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin)
173 72906 : qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin)
174 : END DO
175 : !spin
176 34306 : CALL para_env%sum(rho1_h_spin)
177 34306 : CALL para_env%sum(rho1_s_spin)
178 34306 : CALL para_env%sum(rho1_h_aspin)
179 34306 : CALL para_env%sum(rho1_s_aspin)
180 34306 : qs_charges%total_rho_hard_spin = rho1_h_spin
181 34306 : qs_charges%total_rho_soft_spin = rho1_s_spin
182 34306 : qs_charges%total_rho_hard_abs_spin = rho1_h_aspin
183 34306 : qs_charges%total_rho_soft_abs_spin = rho1_s_aspin
184 :
185 34306 : IF (my_do_rho0) THEN
186 27496 : rho0_mpole%total_rho0_h = -rho0_h_tot
187 :
188 : ! When MGRID is defined within TDDFPT
189 27496 : IF (PRESENT(pw_env_sub)) THEN
190 : ! Find pool
191 2282 : NULLIFY (my_pools, my_rs_grids, my_rs_descs)
192 : CALL pw_env_get(pw_env=pw_env_sub, rs_grids=my_rs_grids, &
193 2282 : rs_descs=my_rs_descs, pw_pools=my_pools)
194 : ! Put the rho0_soft on the global grid
195 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int, my_pools=my_pools, &
196 2282 : my_rs_grids=my_rs_grids, my_rs_descs=my_rs_descs)
197 : ELSE
198 : ! Put the rho0_soft on the global grid
199 25214 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int)
200 : END IF
201 :
202 27496 : IF (ABS(rho0_h_tot) >= 1.0E-5_dp) THEN
203 25462 : IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) > 1.0E-3_dp) THEN
204 1918 : IF (output_unit > 0) THEN
205 959 : WRITE (output_unit, '(/,72("*"))')
206 : WRITE (output_unit, '(T2,A,T66,1E20.8)') &
207 959 : "WARNING: rho0 calculated on the local grid is :", -rho0_h_tot, &
208 1918 : " rho0 calculated on the global grid is :", tot_rs_int
209 : WRITE (output_unit, '(T2,A)') &
210 959 : " bad integration"
211 959 : WRITE (output_unit, '(72("*"),/)')
212 : END IF
213 : END IF
214 : END IF
215 27496 : qs_charges%total_rho0_soft_rspace = tot_rs_int
216 27496 : qs_charges%total_rho0_hard_lebedev = rho0_h_tot
217 27496 : IF (rho0_mpole%do_cneo) THEN
218 : ! put soft tails of quantum nuclear charge densities on the global grid
219 48 : CALL put_rhoz_cneo_s_on_grid(qs_env, rho0_mpole, rhoz_cneo_set, tot_rs_int)
220 48 : IF (ABS(rho0_mpole%tot_rhoz_cneo_s) >= 1.0E-5_dp) THEN
221 40 : IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_mpole%tot_rhoz_cneo_s)) > 1.0E-3_dp) THEN
222 0 : IF (output_unit > 0) THEN
223 0 : WRITE (output_unit, '(/,72("*"))')
224 : WRITE (output_unit, '(T2,A,T66,1E20.8)') &
225 0 : "WARNING: rhoz_cneo_s calculated on the local grid is :", &
226 0 : rho0_mpole%tot_rhoz_cneo_s, &
227 0 : " rhoz_cneo_s calculated on the global grid is :", tot_rs_int
228 : WRITE (output_unit, '(T2,A)') &
229 0 : " bad integration"
230 0 : WRITE (output_unit, '(72("*"),/)')
231 : END IF
232 : END IF
233 : END IF
234 48 : qs_charges%total_rho1_soft_nuc_rspace = tot_rs_int
235 48 : qs_charges%total_rho1_soft_nuc_lebedev = rho0_mpole%tot_rhoz_cneo_s
236 : ELSE
237 27448 : qs_charges%total_rho1_soft_nuc_rspace = 0.0_dp
238 : END IF
239 : ELSE
240 6810 : qs_charges%total_rho0_hard_lebedev = 0.0_dp
241 : END IF
242 : END IF
243 :
244 39484 : DEALLOCATE (rho1_h_tot, rho1_s_tot)
245 :
246 39484 : CALL timestop(handle)
247 :
248 39484 : END SUBROUTINE prepare_gapw_den
249 :
250 : END MODULE qs_gapw_densities
|
