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 gaussian_gridlevels, ONLY: gridlevel_info_type
17 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
18 : section_vals_type,&
19 : section_vals_val_get
20 : USE kinds, ONLY: default_string_length,&
21 : dp
22 : USE pw_env_types, ONLY: pw_env_get,&
23 : pw_env_type
24 : USE pw_methods, ONLY: pw_integrate_function
25 : USE pw_types, ONLY: pw_c1d_gs_type,&
26 : pw_r3d_rs_type
27 : USE qs_environment_types, ONLY: get_qs_env,&
28 : qs_environment_type
29 : USE qs_rho_types, ONLY: qs_rho_get,&
30 : qs_rho_type
31 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
32 : realspace_grid_type,&
33 : rs_grid_create,&
34 : rs_grid_release,&
35 : rs_grid_zero
36 : USE rs_pw_interface, ONLY: density_rs2pw
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_external_density'
44 :
45 : PUBLIC :: external_read_density, read_cube_density
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Computes the external density on the grid
51 : !> \param qs_env ...
52 : !> \date 03.2011
53 : !> \author D. Varsano
54 : ! **************************************************************************************************
55 12836 : SUBROUTINE external_read_density(qs_env)
56 :
57 : TYPE(qs_environment_type), POINTER :: qs_env
58 :
59 : CHARACTER(len=*), PARAMETER :: routineN = 'external_read_density'
60 :
61 : CHARACTER(LEN=default_string_length) :: filename
62 : INTEGER :: handle
63 : TYPE(dft_control_type), POINTER :: dft_control
64 : TYPE(qs_rho_type), POINTER :: rho_external
65 : TYPE(section_vals_type), POINTER :: ext_den_section, input
66 :
67 12836 : CALL timeset(routineN, handle)
68 12836 : NULLIFY (input, ext_den_section, dft_control, rho_external)
69 :
70 : CALL get_qs_env(qs_env, &
71 : rho_external=rho_external, &
72 : input=input, &
73 12836 : dft_control=dft_control)
74 :
75 12836 : IF (dft_control%apply_external_density .AND. dft_control%read_external_density) THEN
76 4 : ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
77 4 : CALL section_vals_val_get(ext_den_section, "FILE_DENSITY", c_val=filename)
78 : CALL read_cube_density(qs_env, rho_external, TRIM(filename), &
79 4 : total_density_sign=1, source_label="ZMP")
80 : END IF
81 :
82 12836 : CALL timestop(handle)
83 :
84 12836 : END SUBROUTINE external_read_density
85 :
86 : ! **************************************************************************************************
87 : !> \brief Read an electron density from a Gaussian cube file into a QS density grid
88 : !> \param qs_env ...
89 : !> \param rho_target target density structure
90 : !> \param filename cube filename
91 : !> \param total_density_sign sign used for rho_target%tot_rho_r
92 : !> \param source_label label used in output
93 : ! **************************************************************************************************
94 24 : SUBROUTINE read_cube_density(qs_env, rho_target, filename, total_density_sign, source_label)
95 : TYPE(qs_environment_type), POINTER :: qs_env
96 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_target
97 : CHARACTER(LEN=*), INTENT(IN) :: filename
98 : INTEGER, INTENT(IN) :: total_density_sign
99 : CHARACTER(LEN=*), INTENT(IN) :: source_label
100 :
101 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cube_density'
102 :
103 : INTEGER :: extunit, handle, i, igrid_level, j, k, &
104 : nat, ndum
105 : INTEGER, DIMENSION(3) :: lbounds, lbounds_local, npoints, &
106 : ubounds, ubounds_local
107 : LOGICAL :: grid_mismatch
108 24 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
109 : REAL(kind=dp), DIMENSION(3) :: cube_origin, voxel
110 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
111 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
112 24 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
113 : TYPE(pw_env_type), POINTER :: pw_env
114 24 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
115 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
116 24 : POINTER :: rs_descs
117 : TYPE(realspace_grid_type), ALLOCATABLE, &
118 24 : DIMENSION(:) :: rs_rho
119 :
120 24 : CALL timeset(routineN, handle)
121 24 : NULLIFY (pw_env, rho_r, rho_g, tot_rho_r, rs_descs)
122 :
123 24 : IF (total_density_sign /= -1 .AND. total_density_sign /= 1) THEN
124 0 : CPABORT("total_density_sign has to be -1 or 1")
125 : END IF
126 :
127 24 : CALL get_qs_env(qs_env, pw_env=pw_env)
128 24 : CALL qs_rho_get(rho_target, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
129 24 : gridlevel_info => pw_env%gridlevel_info
130 24 : CALL pw_env_get(pw_env, rs_descs=rs_descs)
131 :
132 456 : ALLOCATE (rs_rho(gridlevel_info%ngrid_levels))
133 48 : DO igrid_level = 1, gridlevel_info%ngrid_levels
134 24 : CALL rs_grid_create(rs_rho(igrid_level), rs_descs(igrid_level)%rs_desc)
135 48 : CALL rs_grid_zero(rs_rho(igrid_level))
136 : END DO
137 24 : igrid_level = gridlevel_info%ngrid_levels
138 :
139 96 : npoints = rs_descs(igrid_level)%rs_desc%npts
140 96 : lbounds = rs_descs(igrid_level)%rs_desc%lb
141 96 : ubounds = rs_descs(igrid_level)%rs_desc%ub
142 96 : lbounds_local = rho_r(1)%pw_grid%bounds_local(1, :)
143 96 : ubounds_local = rho_r(1)%pw_grid%bounds_local(2, :)
144 72 : ALLOCATE (buffer(lbounds(3):ubounds(3)))
145 :
146 : ASSOCIATE (gid => rho_r(1)%pw_grid%para%group, &
147 : my_rank => rho_r(1)%pw_grid%para%group%mepos)
148 24 : grid_mismatch = .FALSE.
149 24 : IF (my_rank == 0) THEN
150 12 : WRITE (*, FMT="(/,T3,A,A)") TRIM(source_label)//"| Reading electron density: ", TRIM(filename)
151 : CALL open_file(file_name=filename, file_status="OLD", file_form="FORMATTED", &
152 12 : file_action="READ", unit_number=extunit)
153 :
154 12 : READ (extunit, *)
155 12 : READ (extunit, *)
156 12 : READ (extunit, *) nat, cube_origin
157 48 : IF (MAXVAL(ABS(cube_origin)) > 1.0E-4_dp) THEN
158 0 : grid_mismatch = .TRUE.
159 : WRITE (*, FMT="(T3,A,3ES16.8)") TRIM(source_label)// &
160 0 : "| Cube origin is not the CP2K grid origin: ", cube_origin
161 : END IF
162 12 : IF (nat < 0) THEN
163 0 : grid_mismatch = .TRUE.
164 : WRITE (*, FMT="(T3,A)") TRIM(source_label)// &
165 0 : "| Multi-orbital cube files are not supported"
166 : END IF
167 48 : DO i = 1, 3
168 36 : READ (extunit, *) ndum, voxel
169 144 : IF (ndum /= npoints(i) .OR. &
170 12 : MAXVAL(ABS(voxel - rs_descs(igrid_level)%rs_desc%dh(:, i))) > 1.0E-4_dp) THEN
171 0 : grid_mismatch = .TRUE.
172 : WRITE (*, FMT="(T3,A,I0)") TRIM(source_label)// &
173 0 : "| Cube grid does not coincide with CP2K grid along axis ", i
174 0 : WRITE (*, FMT="(T3,A,I0,A,I0)") TRIM(source_label)//"| Grid points: ", &
175 0 : ndum, " instead of ", npoints(i)
176 0 : WRITE (*, FMT="(T3,A,3ES16.8)") TRIM(source_label)//"| Cube vector: ", voxel
177 0 : WRITE (*, FMT="(T3,A,3ES16.8)") TRIM(source_label)//"| CP2K vector: ", &
178 0 : rs_descs(igrid_level)%rs_desc%dh(:, i)
179 : END IF
180 : END DO
181 24 : DO i = 1, ABS(nat)
182 24 : READ (extunit, *)
183 : END DO
184 : END IF
185 :
186 24 : CALL gid%bcast(grid_mismatch, 0)
187 24 : IF (grid_mismatch) THEN
188 0 : IF (my_rank == 0) CALL close_file(unit_number=extunit)
189 0 : CPABORT("Cube density grid is incompatible with the CP2K real-space grid")
190 : END IF
191 :
192 456 : DO i = lbounds(1), ubounds(1)
193 8232 : DO j = lbounds(2), ubounds(2)
194 7776 : IF (my_rank == 0) THEN
195 3888 : READ (extunit, *) (buffer(k), k=lbounds(3), ubounds(3))
196 : END IF
197 7776 : CALL gid%bcast(buffer, 0)
198 : IF ((lbounds_local(1) <= i) .AND. (i <= ubounds_local(1)) .AND. &
199 8208 : (lbounds_local(2) <= j) .AND. (j <= ubounds_local(2))) THEN
200 73872 : rs_rho(igrid_level)%r(i, j, lbounds(3):ubounds(3)) = buffer
201 : END IF
202 : END DO
203 : END DO
204 24 : IF (my_rank == 0) CALL close_file(unit_number=extunit)
205 :
206 24 : CALL density_rs2pw(pw_env, rs_rho, rho=rho_r(1), rho_gspace=rho_g(1))
207 24 : tot_rho_r(1) = pw_integrate_function(rho_r(1), isign=total_density_sign)
208 24 : IF (my_rank == 0) THEN
209 : WRITE (*, FMT="(T3,A,T61,F20.10)") TRIM(source_label)// &
210 12 : "| Integrated electron density:", REAL(total_density_sign, dp)*tot_rho_r(1)
211 : END IF
212 48 : CALL gid%sync()
213 : END ASSOCIATE
214 :
215 48 : DO igrid_level = 1, SIZE(rs_rho)
216 48 : CALL rs_grid_release(rs_rho(igrid_level))
217 : END DO
218 48 : DEALLOCATE (buffer, rs_rho)
219 :
220 24 : CALL timestop(handle)
221 :
222 48 : END SUBROUTINE read_cube_density
223 :
224 : END MODULE qs_external_density
|
