Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 26816 : 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 26816 : INTEGER, DIMENSION(:), POINTER :: atom_list
70 : LOGICAL :: extern, my_do_rho0, paw_atom
71 : REAL(dp) :: rho0_h_tot, tot_rs_int
72 26816 : REAL(dp), DIMENSION(:), POINTER :: rho1_h_tot, rho1_s_tot
73 26816 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
74 : TYPE(dft_control_type), POINTER :: dft_control
75 : TYPE(gapw_control_type), POINTER :: gapw_control
76 : TYPE(mp_para_env_type), POINTER :: para_env
77 26816 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: my_pools
78 : TYPE(qs_charges_type), POINTER :: qs_charges
79 26816 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
80 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
81 26816 : POINTER :: my_rs_descs
82 26816 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: my_rs_grids
83 26816 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
84 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
85 26816 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
86 26816 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
87 :
88 26816 : CALL timeset(routineN, handle)
89 :
90 26816 : NULLIFY (atomic_kind_set)
91 26816 : NULLIFY (my_kind_set)
92 26816 : NULLIFY (dft_control)
93 26816 : NULLIFY (gapw_control)
94 26816 : NULLIFY (para_env)
95 26816 : NULLIFY (atom_list)
96 26816 : NULLIFY (rho0_mpole)
97 26816 : NULLIFY (qs_charges)
98 26816 : NULLIFY (rho1_h_tot, rho1_s_tot)
99 26816 : NULLIFY (rho_atom_set)
100 26816 : NULLIFY (rho0_atom_set)
101 :
102 26816 : my_do_rho0 = .TRUE.
103 26816 : IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0
104 :
105 26816 : output_unit = cp_logger_get_default_io_unit()
106 :
107 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
108 : para_env=para_env, &
109 : qs_charges=qs_charges, &
110 : qs_kind_set=my_kind_set, &
111 : atomic_kind_set=atomic_kind_set, &
112 : rho0_mpole=rho0_mpole, &
113 : rho_atom_set=rho_atom_set, &
114 : rho0_atom_set=rho0_atom_set, &
115 26816 : rhoz_cneo_set=rhoz_cneo_set)
116 :
117 26816 : gapw_control => dft_control%qs_control%gapw_control
118 :
119 : ! If TDDFPT%MGRID is defined, overwrite QS grid info accordingly
120 26816 : IF (PRESENT(local_rho_set)) THEN
121 7760 : rho_atom_set => local_rho_set%rho_atom_set
122 7760 : rhoz_cneo_set => local_rho_set%rhoz_cneo_set
123 7760 : IF (my_do_rho0) THEN
124 2844 : rho0_mpole => local_rho_set%rho0_mpole
125 2844 : rho0_atom_set => local_rho_set%rho0_atom_set
126 : END IF
127 : END IF
128 :
129 26816 : extern = .FALSE.
130 26816 : IF (PRESENT(kind_set_external)) THEN
131 3522 : CPASSERT(ASSOCIATED(kind_set_external))
132 3522 : my_kind_set => kind_set_external
133 3522 : extern = .TRUE.
134 : END IF
135 :
136 26816 : nspins = dft_control%nspins
137 :
138 26816 : rho0_h_tot = 0.0_dp
139 134080 : ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins))
140 58346 : rho1_h_tot = 0.0_dp
141 58346 : rho1_s_tot = 0.0_dp
142 :
143 80430 : DO ikind = 1, SIZE(atomic_kind_set)
144 53614 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
145 53614 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
146 :
147 : !Calculate rho1_h and rho1_s on the radial grids centered on the atomic position
148 53614 : IF (paw_atom) THEN
149 : CALL calculate_rho_atom(para_env, rho_atom_set, my_kind_set(ikind), &
150 49082 : atom_list, natom, nspins, rho1_h_tot, rho1_s_tot)
151 : END IF
152 :
153 : !Calculate rho0_h and rho0_s on the radial grids centered on the atomic position
154 53614 : IF (my_do_rho0) &
155 : CALL calculate_rho0_atom(gapw_control, rho_atom_set, rhoz_cneo_set, rho0_atom_set, &
156 : rho0_mpole, atom_list, natom, ikind, my_kind_set(ikind), &
157 116738 : rho0_h_tot)
158 :
159 : END DO
160 :
161 : !Do not mess with charges if using a non-default kind_set
162 26816 : IF (.NOT. extern) THEN
163 77766 : CALL para_env%sum(rho1_h_tot)
164 77766 : CALL para_env%sum(rho1_s_tot)
165 50530 : DO ispin = 1, nspins
166 27236 : qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin)
167 50530 : qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin)
168 : END DO
169 :
170 23294 : IF (my_do_rho0) THEN
171 18394 : rho0_mpole%total_rho0_h = -rho0_h_tot
172 :
173 : ! When MGRID is defined within TDDFPT
174 18394 : IF (PRESENT(pw_env_sub)) THEN
175 : ! Find pool
176 1304 : NULLIFY (my_pools, my_rs_grids, my_rs_descs)
177 : CALL pw_env_get(pw_env=pw_env_sub, rs_grids=my_rs_grids, &
178 1304 : rs_descs=my_rs_descs, pw_pools=my_pools)
179 : ! Put the rho0_soft on the global grid
180 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int, my_pools=my_pools, &
181 1304 : my_rs_grids=my_rs_grids, my_rs_descs=my_rs_descs)
182 : ELSE
183 : ! Put the rho0_soft on the global grid
184 17090 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int)
185 : END IF
186 :
187 18394 : IF (ABS(rho0_h_tot) >= 1.0E-5_dp) THEN
188 17322 : IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) > 1.0E-3_dp) THEN
189 1778 : IF (output_unit > 0) THEN
190 889 : WRITE (output_unit, '(/,72("*"))')
191 : WRITE (output_unit, '(T2,A,T66,1E20.8)') &
192 889 : "WARNING: rho0 calculated on the local grid is :", -rho0_h_tot, &
193 1778 : " rho0 calculated on the global grid is :", tot_rs_int
194 : WRITE (output_unit, '(T2,A)') &
195 889 : " bad integration"
196 889 : WRITE (output_unit, '(72("*"),/)')
197 : END IF
198 : END IF
199 : END IF
200 18394 : qs_charges%total_rho0_soft_rspace = tot_rs_int
201 18394 : qs_charges%total_rho0_hard_lebedev = rho0_h_tot
202 18394 : IF (rho0_mpole%do_cneo) THEN
203 : ! put soft tails of quantum nuclear charge densities on the global grid
204 48 : CALL put_rhoz_cneo_s_on_grid(qs_env, rho0_mpole, rhoz_cneo_set, tot_rs_int)
205 48 : IF (ABS(rho0_mpole%tot_rhoz_cneo_s) >= 1.0E-5_dp) THEN
206 40 : IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_mpole%tot_rhoz_cneo_s)) > 1.0E-3_dp) THEN
207 0 : IF (output_unit > 0) THEN
208 0 : WRITE (output_unit, '(/,72("*"))')
209 : WRITE (output_unit, '(T2,A,T66,1E20.8)') &
210 0 : "WARNING: rhoz_cneo_s calculated on the local grid is :", &
211 0 : rho0_mpole%tot_rhoz_cneo_s, &
212 0 : " rhoz_cneo_s calculated on the global grid is :", tot_rs_int
213 : WRITE (output_unit, '(T2,A)') &
214 0 : " bad integration"
215 0 : WRITE (output_unit, '(72("*"),/)')
216 : END IF
217 : END IF
218 : END IF
219 48 : qs_charges%total_rho1_soft_nuc_rspace = tot_rs_int
220 48 : qs_charges%total_rho1_soft_nuc_lebedev = rho0_mpole%tot_rhoz_cneo_s
221 : ELSE
222 18346 : qs_charges%total_rho1_soft_nuc_rspace = 0.0_dp
223 : END IF
224 : ELSE
225 4900 : qs_charges%total_rho0_hard_lebedev = 0.0_dp
226 : END IF
227 : END IF
228 :
229 26816 : DEALLOCATE (rho1_h_tot, rho1_s_tot)
230 :
231 26816 : CALL timestop(handle)
232 :
233 26816 : END SUBROUTINE prepare_gapw_den
234 :
235 : END MODULE qs_gapw_densities
|
