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 : !> \brief Routines to handle an external density
9 : !> The external density can be generic and is provided by user input
10 : !> \author D. Varsano
11 : ! **************************************************************************************************
12 : MODULE qs_external_density
13 : USE cp_control_types, ONLY: dft_control_type
14 : USE cp_files, ONLY: close_file,&
15 : open_file
16 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
17 : USE gaussian_gridlevels, ONLY: gridlevel_info_type
18 : USE hfx_types, ONLY: hfx_type
19 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
20 : section_vals_type,&
21 : section_vals_val_get
22 : USE kinds, ONLY: default_string_length,&
23 : dp
24 : USE pw_env_types, ONLY: pw_env_get,&
25 : pw_env_type
26 : USE pw_methods, ONLY: pw_integrate_function
27 : USE pw_types, ONLY: pw_c1d_gs_type,&
28 : pw_r3d_rs_type
29 : USE qs_environment_types, ONLY: get_qs_env,&
30 : qs_environment_type
31 : USE qs_rho_types, ONLY: qs_rho_get,&
32 : qs_rho_set,&
33 : qs_rho_type
34 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
35 : realspace_grid_type,&
36 : rs_grid_create,&
37 : rs_grid_release,&
38 : rs_grid_zero
39 : USE rs_pw_interface, ONLY: density_rs2pw
40 : #include "./base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 :
44 : PRIVATE
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_external_density'
47 :
48 : PUBLIC :: external_read_density, read_cube_density, read_scf_guess_density
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Read a cube density for one-time use as the first SCF input density
54 : !> \param qs_env QS environment
55 : !> \param filename cube filename
56 : ! **************************************************************************************************
57 4 : SUBROUTINE read_scf_guess_density(qs_env, filename)
58 : TYPE(qs_environment_type), POINTER :: qs_env
59 : CHARACTER(LEN=*), INTENT(IN) :: filename
60 :
61 : TYPE(dft_control_type), POINTER :: dft_control
62 4 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
63 : TYPE(qs_rho_type), POINTER :: rho
64 :
65 4 : NULLIFY (dft_control, rho, x_data)
66 :
67 4 : IF (LEN_TRIM(filename) == 0) THEN
68 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY requires a cube filename")
69 : END IF
70 :
71 4 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho, x_data=x_data)
72 4 : IF (.NOT. dft_control%qs_control%gpw) THEN
73 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY currently supports GPW only")
74 : END IF
75 4 : IF (dft_control%nspins /= 1) THEN
76 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY currently requires a spin-restricted calculation")
77 : END IF
78 4 : IF (dft_control%use_kinetic_energy_density) THEN
79 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY does not provide a kinetic-energy density")
80 : END IF
81 4 : IF (dft_control%drho_by_collocation) THEN
82 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY is incompatible with DRHO_BY_COLLOCATION")
83 : END IF
84 4 : IF (dft_control%dft_plus_u) THEN
85 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY currently does not support DFT+U")
86 : END IF
87 4 : IF (ASSOCIATED(x_data)) THEN
88 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY does not support HFX")
89 : END IF
90 4 : IF (dft_control%apply_external_density) THEN
91 0 : CPABORT("SCF_GUESS EXTERNAL_DENSITY cannot be combined with DFT%EXTERNAL_DENSITY")
92 : END IF
93 :
94 : CALL read_cube_density(qs_env, rho, TRIM(filename), total_density_sign=-1, &
95 4 : source_label="SCF GUESS")
96 : CALL qs_rho_set(rho, rho_r_valid=.TRUE., rho_g_valid=.TRUE., &
97 : drho_r_valid=.FALSE., drho_g_valid=.FALSE., &
98 4 : tau_r_valid=.FALSE., tau_g_valid=.FALSE., soft_valid=.FALSE.)
99 :
100 4 : END SUBROUTINE read_scf_guess_density
101 :
102 : ! **************************************************************************************************
103 : !> \brief Computes the external density on the grid
104 : !> \param qs_env ...
105 : !> \date 03.2011
106 : !> \author D. Varsano
107 : ! **************************************************************************************************
108 12880 : SUBROUTINE external_read_density(qs_env)
109 :
110 : TYPE(qs_environment_type), POINTER :: qs_env
111 :
112 : CHARACTER(len=*), PARAMETER :: routineN = 'external_read_density'
113 :
114 : CHARACTER(LEN=default_string_length) :: filename
115 : INTEGER :: handle
116 : TYPE(dft_control_type), POINTER :: dft_control
117 : TYPE(qs_rho_type), POINTER :: rho_external
118 : TYPE(section_vals_type), POINTER :: ext_den_section, input
119 :
120 12880 : CALL timeset(routineN, handle)
121 12880 : NULLIFY (input, ext_den_section, dft_control, rho_external)
122 :
123 : CALL get_qs_env(qs_env, &
124 : rho_external=rho_external, &
125 : input=input, &
126 12880 : dft_control=dft_control)
127 :
128 12880 : IF (dft_control%apply_external_density .AND. dft_control%read_external_density) THEN
129 4 : ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
130 4 : CALL section_vals_val_get(ext_den_section, "FILE_DENSITY", c_val=filename)
131 : CALL read_cube_density(qs_env, rho_external, TRIM(filename), &
132 4 : total_density_sign=1, source_label="ZMP")
133 : END IF
134 :
135 12880 : CALL timestop(handle)
136 :
137 12880 : END SUBROUTINE external_read_density
138 :
139 : ! **************************************************************************************************
140 : !> \brief Read an electron density from a Gaussian cube file into a QS density grid
141 : !> \param qs_env ...
142 : !> \param rho_target target density structure
143 : !> \param filename cube filename
144 : !> \param total_density_sign sign used for rho_target%tot_rho_r
145 : !> \param source_label label used in output
146 : ! **************************************************************************************************
147 28 : SUBROUTINE read_cube_density(qs_env, rho_target, filename, total_density_sign, source_label)
148 : TYPE(qs_environment_type), POINTER :: qs_env
149 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_target
150 : CHARACTER(LEN=*), INTENT(IN) :: filename
151 : INTEGER, INTENT(IN) :: total_density_sign
152 : CHARACTER(LEN=*), INTENT(IN) :: source_label
153 :
154 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cube_density'
155 :
156 : INTEGER :: extunit, handle, i, igrid_level, j, k, &
157 : nat, ndum, output_unit
158 : INTEGER, DIMENSION(3) :: lbounds, lbounds_local, npoints, &
159 : ubounds, ubounds_local
160 : LOGICAL :: grid_mismatch
161 28 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
162 : REAL(kind=dp), DIMENSION(3) :: cube_origin, voxel
163 28 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
164 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
165 28 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
166 : TYPE(pw_env_type), POINTER :: pw_env
167 28 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
168 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
169 28 : POINTER :: rs_descs
170 : TYPE(realspace_grid_type), ALLOCATABLE, &
171 28 : DIMENSION(:) :: rs_rho
172 :
173 28 : CALL timeset(routineN, handle)
174 28 : NULLIFY (pw_env, rho_r, rho_g, tot_rho_r, rs_descs)
175 28 : output_unit = cp_logger_get_default_io_unit()
176 :
177 28 : IF (total_density_sign /= -1 .AND. total_density_sign /= 1) THEN
178 0 : CPABORT("total_density_sign has to be -1 or 1")
179 : END IF
180 :
181 28 : CALL get_qs_env(qs_env, pw_env=pw_env)
182 28 : CALL qs_rho_get(rho_target, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
183 28 : gridlevel_info => pw_env%gridlevel_info
184 28 : CALL pw_env_get(pw_env, rs_descs=rs_descs)
185 :
186 532 : ALLOCATE (rs_rho(gridlevel_info%ngrid_levels))
187 56 : DO igrid_level = 1, gridlevel_info%ngrid_levels
188 28 : CALL rs_grid_create(rs_rho(igrid_level), rs_descs(igrid_level)%rs_desc)
189 56 : CALL rs_grid_zero(rs_rho(igrid_level))
190 : END DO
191 28 : igrid_level = gridlevel_info%ngrid_levels
192 :
193 112 : npoints = rs_descs(igrid_level)%rs_desc%npts
194 112 : lbounds = rs_descs(igrid_level)%rs_desc%lb
195 112 : ubounds = rs_descs(igrid_level)%rs_desc%ub
196 112 : lbounds_local = rho_r(1)%pw_grid%bounds_local(1, :)
197 112 : ubounds_local = rho_r(1)%pw_grid%bounds_local(2, :)
198 84 : ALLOCATE (buffer(lbounds(3):ubounds(3)))
199 :
200 : ASSOCIATE (gid => rho_r(1)%pw_grid%para%group, &
201 : my_rank => rho_r(1)%pw_grid%para%group%mepos)
202 28 : grid_mismatch = .FALSE.
203 28 : IF (my_rank == 0) THEN
204 : WRITE (output_unit, FMT="(/,T3,A,A)") &
205 14 : TRIM(source_label)//"| Reading electron density: ", TRIM(filename)
206 : CALL open_file(file_name=filename, file_status="OLD", file_form="FORMATTED", &
207 14 : file_action="READ", unit_number=extunit)
208 :
209 14 : READ (extunit, *)
210 14 : READ (extunit, *)
211 14 : READ (extunit, *) nat, cube_origin
212 56 : IF (MAXVAL(ABS(cube_origin)) > 1.0E-4_dp) THEN
213 0 : grid_mismatch = .TRUE.
214 : WRITE (output_unit, FMT="(T3,A,3ES16.8)") TRIM(source_label)// &
215 0 : "| Cube origin is not the CP2K grid origin: ", cube_origin
216 : END IF
217 14 : IF (nat < 0) THEN
218 0 : grid_mismatch = .TRUE.
219 : WRITE (output_unit, FMT="(T3,A)") TRIM(source_label)// &
220 0 : "| Multi-orbital cube files are not supported"
221 : END IF
222 56 : DO i = 1, 3
223 42 : READ (extunit, *) ndum, voxel
224 168 : IF (ndum /= npoints(i) .OR. &
225 14 : MAXVAL(ABS(voxel - rs_descs(igrid_level)%rs_desc%dh(:, i))) > 1.0E-4_dp) THEN
226 0 : grid_mismatch = .TRUE.
227 : WRITE (output_unit, FMT="(T3,A,I0)") TRIM(source_label)// &
228 0 : "| Cube grid does not coincide with CP2K grid along axis ", i
229 0 : WRITE (output_unit, FMT="(T3,A,I0,A,I0)") TRIM(source_label)//"| Grid points: ", &
230 0 : ndum, " instead of ", npoints(i)
231 0 : WRITE (output_unit, FMT="(T3,A,3ES16.8)") TRIM(source_label)//"| Cube vector: ", voxel
232 0 : WRITE (output_unit, FMT="(T3,A,3ES16.8)") TRIM(source_label)//"| CP2K vector: ", &
233 0 : rs_descs(igrid_level)%rs_desc%dh(:, i)
234 : END IF
235 : END DO
236 28 : DO i = 1, ABS(nat)
237 28 : READ (extunit, *)
238 : END DO
239 : END IF
240 :
241 28 : CALL gid%bcast(grid_mismatch, 0)
242 28 : IF (grid_mismatch) THEN
243 0 : IF (my_rank == 0) CALL close_file(unit_number=extunit)
244 0 : CPABORT("Cube density grid is incompatible with the CP2K real-space grid")
245 : END IF
246 :
247 532 : DO i = lbounds(1), ubounds(1)
248 9604 : DO j = lbounds(2), ubounds(2)
249 9072 : IF (my_rank == 0) THEN
250 4536 : READ (extunit, *) (buffer(k), k=lbounds(3), ubounds(3))
251 : END IF
252 9072 : CALL gid%bcast(buffer, 0)
253 : IF ((lbounds_local(1) <= i) .AND. (i <= ubounds_local(1)) .AND. &
254 9576 : (lbounds_local(2) <= j) .AND. (j <= ubounds_local(2))) THEN
255 86184 : rs_rho(igrid_level)%r(i, j, lbounds(3):ubounds(3)) = buffer
256 : END IF
257 : END DO
258 : END DO
259 28 : IF (my_rank == 0) CALL close_file(unit_number=extunit)
260 :
261 28 : CALL density_rs2pw(pw_env, rs_rho, rho=rho_r(1), rho_gspace=rho_g(1))
262 28 : tot_rho_r(1) = pw_integrate_function(rho_r(1), isign=total_density_sign)
263 28 : IF (my_rank == 0) THEN
264 : WRITE (output_unit, FMT="(T3,A,T61,F20.10)") TRIM(source_label)// &
265 14 : "| Integrated electron density:", REAL(total_density_sign, dp)*tot_rho_r(1)
266 : END IF
267 56 : CALL gid%sync()
268 : END ASSOCIATE
269 :
270 56 : DO igrid_level = 1, SIZE(rs_rho)
271 56 : CALL rs_grid_release(rs_rho(igrid_level))
272 : END DO
273 56 : DEALLOCATE (buffer, rs_rho)
274 :
275 28 : CALL timestop(handle)
276 :
277 56 : END SUBROUTINE read_cube_density
278 :
279 : END MODULE qs_external_density
|
