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_environment_types, ONLY: get_qs_env,&
22 : qs_environment_type
23 : USE qs_kind_types, ONLY: get_qs_kind,&
24 : qs_kind_type
25 : USE qs_local_rho_types, ONLY: local_rho_type
26 : USE qs_rho0_ggrid, ONLY: put_rho0_on_grid
27 : USE qs_rho0_methods, ONLY: calculate_rho0_atom
28 : USE qs_rho0_types, ONLY: rho0_atom_type,&
29 : rho0_mpole_type
30 : USE qs_rho_atom_methods, ONLY: calculate_rho_atom
31 : USE qs_rho_atom_types, ONLY: rho_atom_type
32 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
33 : realspace_grid_type
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gapw_densities'
41 :
42 : PUBLIC :: prepare_gapw_den
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief ...
48 : !> \param qs_env ...
49 : !> \param local_rho_set ...
50 : !> \param do_rho0 ...
51 : !> \param kind_set_external can be provided to use different projectors/grids/basis than the default
52 : !> \param pw_env_sub ...
53 : ! **************************************************************************************************
54 24156 : SUBROUTINE prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external, pw_env_sub)
55 :
56 : TYPE(qs_environment_type), POINTER :: qs_env
57 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
58 : LOGICAL, INTENT(IN), OPTIONAL :: do_rho0
59 : TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
60 : POINTER :: kind_set_external
61 : TYPE(pw_env_type), OPTIONAL :: pw_env_sub
62 :
63 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_gapw_den'
64 :
65 : INTEGER :: handle, ikind, ispin, natom, nspins, &
66 : output_unit
67 24156 : INTEGER, DIMENSION(:), POINTER :: atom_list
68 : LOGICAL :: extern, my_do_rho0, paw_atom
69 : REAL(dp) :: rho0_h_tot, tot_rs_int
70 24156 : REAL(dp), DIMENSION(:), POINTER :: rho1_h_tot, rho1_s_tot
71 24156 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
72 : TYPE(dft_control_type), POINTER :: dft_control
73 : TYPE(gapw_control_type), POINTER :: gapw_control
74 : TYPE(mp_para_env_type), POINTER :: para_env
75 24156 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: my_pools
76 : TYPE(qs_charges_type), POINTER :: qs_charges
77 24156 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
78 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
79 24156 : POINTER :: my_rs_descs
80 24156 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: my_rs_grids
81 24156 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
82 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
83 24156 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
84 :
85 24156 : CALL timeset(routineN, handle)
86 :
87 24156 : NULLIFY (atomic_kind_set)
88 24156 : NULLIFY (my_kind_set)
89 24156 : NULLIFY (dft_control)
90 24156 : NULLIFY (gapw_control)
91 24156 : NULLIFY (para_env)
92 24156 : NULLIFY (atom_list)
93 24156 : NULLIFY (rho0_mpole)
94 24156 : NULLIFY (qs_charges)
95 24156 : NULLIFY (rho1_h_tot, rho1_s_tot)
96 24156 : NULLIFY (rho_atom_set)
97 24156 : NULLIFY (rho0_atom_set)
98 :
99 24156 : my_do_rho0 = .TRUE.
100 24156 : IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0
101 :
102 24156 : output_unit = cp_logger_get_default_io_unit()
103 :
104 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
105 : para_env=para_env, &
106 : qs_charges=qs_charges, &
107 : qs_kind_set=my_kind_set, &
108 : atomic_kind_set=atomic_kind_set, &
109 : rho0_mpole=rho0_mpole, &
110 : rho_atom_set=rho_atom_set, &
111 24156 : rho0_atom_set=rho0_atom_set)
112 :
113 24156 : gapw_control => dft_control%qs_control%gapw_control
114 :
115 : ! If TDDFPT%MGRID is defined, overwrite QS grid info accordingly
116 24156 : IF (PRESENT(local_rho_set)) THEN
117 6676 : rho_atom_set => local_rho_set%rho_atom_set
118 6676 : IF (my_do_rho0) THEN
119 2322 : rho0_mpole => local_rho_set%rho0_mpole
120 2322 : rho0_atom_set => local_rho_set%rho0_atom_set
121 : END IF
122 : END IF
123 :
124 24156 : extern = .FALSE.
125 24156 : IF (PRESENT(kind_set_external)) THEN
126 3204 : CPASSERT(ASSOCIATED(kind_set_external))
127 3204 : my_kind_set => kind_set_external
128 3204 : extern = .TRUE.
129 : END IF
130 :
131 24156 : nspins = dft_control%nspins
132 :
133 24156 : rho0_h_tot = 0.0_dp
134 120780 : ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins))
135 52850 : rho1_h_tot = 0.0_dp
136 52850 : rho1_s_tot = 0.0_dp
137 :
138 73984 : DO ikind = 1, SIZE(atomic_kind_set)
139 49828 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
140 49828 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
141 :
142 : !Calculate rho1_h and rho1_s on the radial grids centered on the atomic position
143 49828 : IF (paw_atom) THEN
144 : CALL calculate_rho_atom(para_env, rho_atom_set, my_kind_set(ikind), &
145 45466 : atom_list, natom, nspins, rho1_h_tot, rho1_s_tot)
146 : END IF
147 :
148 : !Calculate rho0_h and rho0_s on the radial grids centered on the atomic position
149 49828 : IF (my_do_rho0) &
150 : CALL calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, rho0_mpole, &
151 107904 : atom_list, natom, ikind, my_kind_set(ikind), rho0_h_tot)
152 :
153 : END DO
154 :
155 : !Do not mess with charges if using a non-default kind_set
156 24156 : IF (.NOT. extern) THEN
157 70388 : CALL para_env%sum(rho1_h_tot)
158 70388 : CALL para_env%sum(rho1_s_tot)
159 45670 : DO ispin = 1, nspins
160 24718 : qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin)
161 45670 : qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin)
162 : END DO
163 :
164 20952 : IF (my_do_rho0) THEN
165 16842 : rho0_mpole%total_rho0_h = -rho0_h_tot
166 :
167 : ! When MGRID is defined within TDDFPT
168 16842 : IF (PRESENT(pw_env_sub)) THEN
169 : ! Find pool
170 1098 : NULLIFY (my_pools, my_rs_grids, my_rs_descs)
171 : CALL pw_env_get(pw_env=pw_env_sub, rs_grids=my_rs_grids, &
172 1098 : rs_descs=my_rs_descs, pw_pools=my_pools)
173 : ! Put the rho0_soft on the global grid
174 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int, my_pools=my_pools, &
175 1098 : my_rs_grids=my_rs_grids, my_rs_descs=my_rs_descs)
176 : ELSE
177 : ! Put the rho0_soft on the global grid
178 15744 : CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int)
179 : END IF
180 :
181 16842 : IF (ABS(rho0_h_tot) .GE. 1.0E-5_dp) THEN
182 15842 : IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) .GT. 1.0E-3_dp) THEN
183 1778 : IF (output_unit > 0) THEN
184 889 : WRITE (output_unit, '(/,72("*"))')
185 : WRITE (output_unit, '(T2,A,T66,1E20.8)') &
186 889 : "WARNING: rho0 calculated on the local grid is :", -rho0_h_tot, &
187 1778 : " rho0 calculated on the global grid is :", tot_rs_int
188 : WRITE (output_unit, '(T2,A)') &
189 889 : " bad integration"
190 889 : WRITE (output_unit, '(72("*"),/)')
191 : END IF
192 : END IF
193 : END IF
194 16842 : qs_charges%total_rho0_soft_rspace = tot_rs_int
195 16842 : qs_charges%total_rho0_hard_lebedev = rho0_h_tot
196 : ELSE
197 4110 : qs_charges%total_rho0_hard_lebedev = 0.0_dp
198 : END IF
199 : END IF
200 :
201 24156 : DEALLOCATE (rho1_h_tot, rho1_s_tot)
202 :
203 24156 : CALL timestop(handle)
204 :
205 24156 : END SUBROUTINE prepare_gapw_den
206 :
207 : END MODULE qs_gapw_densities
|
