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 : MODULE fde_methods
8 : USE cp_control_types, ONLY: dft_control_type,&
9 : fde_control_type
10 : USE cp_files, ONLY: close_file,&
11 : open_file
12 : USE cp_log_handling, ONLY: cp_to_string
13 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
14 : section_vals_type
15 : USE kinds, ONLY: dp
16 : USE message_passing, ONLY: mp_comm_type,&
17 : mp_para_env_type
18 : USE pw_env_types, ONLY: pw_env_get,&
19 : pw_env_type
20 : USE pw_grid_types, ONLY: PW_MODE_LOCAL
21 : USE pw_methods, ONLY: pw_axpy,&
22 : pw_copy,&
23 : pw_transfer,&
24 : pw_zero
25 : USE pw_poisson_methods, ONLY: pw_poisson_solve
26 : USE pw_poisson_types, ONLY: pw_poisson_type
27 : USE pw_pool_types, ONLY: pw_pool_type
28 : USE pw_types, ONLY: pw_c1d_gs_type,&
29 : pw_r3d_rs_type
30 : USE qs_environment_types, ONLY: get_qs_env,&
31 : qs_environment_type
32 : USE qs_ks_methods, ONLY: calc_rho_tot_gspace
33 : USE qs_ks_types, ONLY: get_ks_env,&
34 : qs_ks_env_type
35 : USE qs_rho_types, ONLY: qs_rho_get,&
36 : qs_rho_type
37 : USE xc, ONLY: xc_vxc_pw_create
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fde_methods'
44 :
45 : PUBLIC :: &
46 : fde_build_v_emb, &
47 : fde_get_env_density, &
48 : fde_spawn, &
49 : fde_v_es, &
50 : fde_write_pw, &
51 : fde_read_pw, &
52 : fde_read_results
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Build the electrostatic potential of a QS calculation
58 : !> \param qs_env ...
59 : !> \param v_es ...
60 : ! **************************************************************************************************
61 16 : SUBROUTINE fde_v_es(qs_env, v_es)
62 : TYPE(qs_environment_type), POINTER :: qs_env
63 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: v_es
64 :
65 : REAL(KIND=dp) :: e_es
66 : TYPE(pw_c1d_gs_type) :: rho_tot_g, v_es_g
67 : TYPE(pw_env_type), POINTER :: pw_env
68 : TYPE(pw_poisson_type), POINTER :: poisson_env
69 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
70 : TYPE(qs_rho_type), POINTER :: rho
71 :
72 8 : NULLIFY (pw_env, poisson_env, auxbas_pw_pool, rho)
73 8 : CALL get_qs_env(qs_env, pw_env=pw_env, rho=rho)
74 8 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
75 :
76 8 : CALL auxbas_pw_pool%create_pw(rho_tot_g)
77 8 : CALL pw_zero(rho_tot_g)
78 8 : CALL calc_rho_tot_gspace(rho_tot_g, qs_env, rho)
79 :
80 8 : CALL auxbas_pw_pool%create_pw(v_es_g)
81 8 : CALL pw_poisson_solve(poisson_env, rho_tot_g, e_es, v_es_g)
82 :
83 8 : CALL auxbas_pw_pool%create_pw(v_es)
84 8 : CALL pw_transfer(v_es_g, v_es)
85 :
86 8 : CALL auxbas_pw_pool%give_back_pw(rho_tot_g)
87 8 : CALL auxbas_pw_pool%give_back_pw(v_es_g)
88 :
89 8 : END SUBROUTINE fde_v_es
90 :
91 : ! **************************************************************************************************
92 : !> \brief Copy the electron density from a QS calculation into a new grid
93 : !> \param qs_env ...
94 : !> \param rho ...
95 : ! **************************************************************************************************
96 4 : SUBROUTINE fde_get_env_density(qs_env, rho)
97 : TYPE(qs_environment_type), POINTER :: qs_env
98 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: rho
99 :
100 : INTEGER :: ispin, nspins
101 : TYPE(dft_control_type), POINTER :: dft_control
102 : TYPE(pw_env_type), POINTER :: pw_env
103 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
104 4 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
105 : TYPE(qs_rho_type), POINTER :: rho_qs
106 :
107 4 : NULLIFY (dft_control, pw_env, auxbas_pw_pool, rho_r, rho_qs)
108 4 : CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, rho=rho_qs)
109 4 : nspins = dft_control%nspins
110 :
111 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
112 4 : CALL auxbas_pw_pool%create_pw(rho)
113 4 : CALL pw_zero(rho)
114 :
115 4 : CALL qs_rho_get(rho_qs, rho_r=rho_r)
116 8 : DO ispin = 1, nspins
117 8 : CALL pw_axpy(rho_r(ispin), rho, 1.0_dp)
118 : END DO
119 :
120 4 : END SUBROUTINE fde_get_env_density
121 :
122 : ! **************************************************************************************************
123 : !> \brief Evaluate an XC+KEDF potential and energy for a density
124 : !> \param ks_env ...
125 : !> \param pw_env ...
126 : !> \param xc_section ...
127 : !> \param rho ...
128 : !> \param v_xc ...
129 : !> \param e_xc ...
130 : ! **************************************************************************************************
131 24 : SUBROUTINE fde_vxc_on_grid(ks_env, pw_env, xc_section, rho, v_xc, e_xc)
132 : TYPE(qs_ks_env_type), POINTER :: ks_env
133 : TYPE(pw_env_type), POINTER :: pw_env
134 : TYPE(section_vals_type), POINTER :: xc_section
135 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
136 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: v_xc
137 : REAL(KIND=dp), INTENT(OUT) :: e_xc
138 :
139 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc
140 : TYPE(pw_c1d_gs_type) :: rho_g_single
141 24 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
142 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
143 24 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc_rho, vxc_tau
144 : TYPE(pw_r3d_rs_type), POINTER :: weights
145 :
146 24 : NULLIFY (auxbas_pw_pool, xc_pw_pool, rho_r, rho_g, tau, vxc_rho, vxc_tau, weights)
147 24 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
148 24 : CALL get_ks_env(ks_env, xcint_weights=weights)
149 :
150 72 : ALLOCATE (rho_r(1), rho_g(1))
151 24 : rho_r(1) = rho
152 24 : CALL auxbas_pw_pool%create_pw(rho_g_single)
153 24 : CALL pw_transfer(rho, rho_g_single)
154 24 : rho_g(1) = rho_g_single
155 :
156 : CALL xc_vxc_pw_create(vxc_rho=vxc_rho, &
157 : vxc_tau=vxc_tau, &
158 : exc=e_xc, &
159 : rho_r=rho_r, &
160 : rho_g=rho_g, &
161 : tau=tau, &
162 : xc_section=xc_section, &
163 : weights=weights, &
164 : pw_pool=xc_pw_pool, &
165 : compute_virial=.FALSE., &
166 24 : virial_xc=virial_xc)
167 :
168 24 : CALL auxbas_pw_pool%create_pw(v_xc)
169 24 : CALL pw_copy(vxc_rho(1), v_xc)
170 :
171 24 : CALL xc_pw_pool%give_back_pw(vxc_rho(1))
172 24 : DEALLOCATE (vxc_rho)
173 24 : IF (ASSOCIATED(vxc_tau)) THEN
174 0 : CALL xc_pw_pool%give_back_pw(vxc_tau(1))
175 0 : DEALLOCATE (vxc_tau)
176 : END IF
177 24 : CALL auxbas_pw_pool%give_back_pw(rho_g_single)
178 24 : DEALLOCATE (rho_r, rho_g)
179 :
180 24 : END SUBROUTINE fde_vxc_on_grid
181 :
182 : ! **************************************************************************************************
183 : !> \brief Build the embedding potential and evaluate nonadditive terms of a QS calculation.
184 : !> \param qs_env ...
185 : !> \param rho ...
186 : !> \param v_emb ...
187 : !> \param e_xc_na ...
188 : !> \param e_kin_na ...
189 : ! **************************************************************************************************
190 36 : SUBROUTINE fde_build_v_emb(qs_env, rho, v_emb, e_xc_na, e_kin_na)
191 : TYPE(qs_environment_type), POINTER :: qs_env
192 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
193 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: v_emb
194 : REAL(KIND=dp), INTENT(OUT) :: e_xc_na, e_kin_na
195 :
196 : REAL(KIND=dp) :: exc_a, exc_ab, exc_b
197 : TYPE(dft_control_type), POINTER :: dft_control
198 : TYPE(pw_env_type), POINTER :: pw_env
199 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
200 : TYPE(pw_r3d_rs_type) :: rho_ab, rho_b_r, v_a, v_ab, v_tmp
201 : TYPE(qs_ks_env_type), POINTER :: ks_env
202 : TYPE(section_vals_type), POINTER :: input, kg_xc_section, xc_section
203 :
204 4 : NULLIFY (dft_control, pw_env, auxbas_pw_pool, ks_env, input, xc_section, kg_xc_section)
205 4 : CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, ks_env=ks_env, input=input)
206 4 : CPASSERT(dft_control%nspins == 1)
207 4 : CPASSERT(.NOT. dft_control%qs_control%gapw)
208 :
209 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
210 :
211 4 : CALL fde_get_env_density(qs_env, rho_b_r)
212 4 : CALL auxbas_pw_pool%create_pw(rho_ab)
213 4 : CALL pw_copy(rho_b_r, rho_ab)
214 4 : CALL pw_axpy(rho, rho_ab, 1.0_dp)
215 :
216 4 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
217 4 : kg_xc_section => section_vals_get_subs_vals(input, "DFT%KG_METHOD%XC")
218 :
219 4 : CALL fde_v_es(qs_env, v_emb)
220 :
221 4 : CALL fde_vxc_on_grid(ks_env, pw_env, xc_section, rho_ab, v_ab, exc_ab)
222 4 : CALL fde_vxc_on_grid(ks_env, pw_env, xc_section, rho, v_a, exc_a)
223 4 : CALL fde_vxc_on_grid(ks_env, pw_env, xc_section, rho_b_r, v_tmp, exc_b) ! energy only
224 4 : CALL auxbas_pw_pool%give_back_pw(v_tmp)
225 4 : CALL pw_axpy(v_ab, v_emb, 1.0_dp)
226 4 : CALL pw_axpy(v_a, v_emb, -1.0_dp)
227 4 : CALL auxbas_pw_pool%give_back_pw(v_ab)
228 4 : CALL auxbas_pw_pool%give_back_pw(v_a)
229 4 : e_xc_na = exc_ab - exc_a - exc_b
230 :
231 4 : CALL fde_vxc_on_grid(ks_env, pw_env, kg_xc_section, rho_ab, v_ab, exc_ab)
232 4 : CALL fde_vxc_on_grid(ks_env, pw_env, kg_xc_section, rho, v_a, exc_a)
233 4 : CALL fde_vxc_on_grid(ks_env, pw_env, kg_xc_section, rho_b_r, v_tmp, exc_b) ! energy only
234 4 : CALL auxbas_pw_pool%give_back_pw(v_tmp)
235 4 : CALL pw_axpy(v_ab, v_emb, 1.0_dp)
236 4 : CALL pw_axpy(v_a, v_emb, -1.0_dp)
237 4 : CALL auxbas_pw_pool%give_back_pw(v_ab)
238 4 : CALL auxbas_pw_pool%give_back_pw(v_a)
239 4 : e_kin_na = exc_ab - exc_a - exc_b
240 :
241 4 : CALL auxbas_pw_pool%give_back_pw(rho_ab)
242 4 : CALL auxbas_pw_pool%give_back_pw(rho_b_r)
243 :
244 4 : END SUBROUTINE fde_build_v_emb
245 :
246 : ! **************************************************************************************************
247 : !> \brief Spawns the subsystem solver. Only the source rank launches the process; the other ranks
248 : !> wait for it on the broadcasts below
249 : !> \param fde_control ...
250 : !> \param input_file ...
251 : !> \param para_env ...
252 : ! **************************************************************************************************
253 4 : SUBROUTINE fde_spawn(fde_control, input_file, para_env)
254 : TYPE(fde_control_type), POINTER :: fde_control
255 : CHARACTER(len=*), INTENT(IN) :: input_file
256 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
257 :
258 4 : CHARACTER(len=:), ALLOCATABLE :: command
259 : INTEGER :: cmdstat, exitstat
260 :
261 4 : IF (LEN_TRIM(fde_control%workdir) > 0) THEN
262 : command = "cd "//TRIM(fde_control%workdir)//" && "// &
263 0 : TRIM(fde_control%command)//" "//TRIM(input_file)
264 : ELSE
265 4 : command = TRIM(fde_control%command)//" "//TRIM(input_file)
266 : END IF
267 :
268 4 : exitstat = 0
269 4 : cmdstat = 0
270 4 : IF (para_env%is_source()) THEN
271 2 : CALL EXECUTE_COMMAND_LINE(command, wait=.TRUE., exitstat=exitstat, cmdstat=cmdstat)
272 : END IF
273 :
274 4 : CALL para_env%bcast(cmdstat)
275 4 : CALL para_env%bcast(exitstat)
276 :
277 4 : IF (cmdstat /= 0) THEN
278 0 : CPABORT("FDE: failed to launch command: "//command)
279 : END IF
280 4 : IF (exitstat /= 0) THEN
281 0 : CPABORT("FDE: command returned a non-zero exit status: "//TRIM(cp_to_string(exitstat)))
282 : END IF
283 :
284 4 : END SUBROUTINE fde_spawn
285 :
286 : ! **************************************************************************************************
287 : !> \brief Write grid points to disk in the Molcas format, optionally with weights and values.
288 : !> write_value=.TRUE.: x y z weight value
289 : !> write_value=.FALSE.: x y z
290 : !> \param pw ...
291 : !> \param filename ...
292 : !> \param write_value ...
293 : ! **************************************************************************************************
294 8 : SUBROUTINE fde_write_pw(pw, filename, write_value)
295 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
296 : CHARACTER(len=*), INTENT(IN) :: filename
297 : LOGICAL, INTENT(IN) :: write_value
298 :
299 : CHARACTER(len=*), PARAMETER :: routineN = 'fde_write_pw'
300 :
301 : INTEGER :: dest, handle, I1, I2, I3, ip, L1, L2, &
302 : L3, mepos, num_pe, source, tag, U1, &
303 : U2, U3, unit_nr
304 : REAL(KIND=dp) :: dvol, x, y, z
305 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buf
306 : TYPE(mp_comm_type) :: group
307 :
308 8 : CALL timeset(routineN, handle)
309 :
310 8 : L1 = pw%pw_grid%bounds(1, 1)
311 8 : U1 = pw%pw_grid%bounds(2, 1)
312 8 : L2 = pw%pw_grid%bounds(1, 2)
313 8 : U2 = pw%pw_grid%bounds(2, 2)
314 8 : L3 = pw%pw_grid%bounds(1, 3)
315 8 : U3 = pw%pw_grid%bounds(2, 3)
316 8 : dvol = pw%pw_grid%dvol
317 :
318 8 : group = pw%pw_grid%para%group
319 8 : mepos = pw%pw_grid%para%group%mepos
320 8 : num_pe = pw%pw_grid%para%group%num_pe
321 8 : tag = 1
322 8 : dest = 0
323 :
324 8 : unit_nr = -1
325 8 : IF (mepos == dest) THEN
326 : CALL open_file(file_name=TRIM(filename), file_status="REPLACE", file_form="FORMATTED", &
327 4 : file_action="WRITE", unit_number=unit_nr)
328 4 : WRITE (unit_nr, '(I10)') pw%pw_grid%ngpts
329 : END IF
330 :
331 24 : ALLOCATE (buf(L3:U3))
332 :
333 152 : DO I1 = L1, U1
334 2744 : DO I2 = L2, U2
335 2592 : IF (pw%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
336 2592 : source = dest
337 7776 : DO ip = 0, num_pe - 1
338 : IF (pw%pw_grid%para%bo(1, 1, ip, 1) <= I1 - L1 + 1 &
339 : .AND. pw%pw_grid%para%bo(2, 1, ip, 1) >= I1 - L1 + 1 &
340 : .AND. pw%pw_grid%para%bo(1, 2, ip, 1) <= I2 - L2 + 1 &
341 7776 : .AND. pw%pw_grid%para%bo(2, 2, ip, 1) >= I2 - L2 + 1) THEN
342 2592 : source = ip
343 : END IF
344 : END DO
345 : ELSE
346 0 : source = dest
347 : END IF
348 :
349 2592 : IF (source == dest) THEN
350 12960 : IF (mepos == source) buf(:) = pw%array(I1, I2, :)
351 : ELSE
352 1296 : IF (mepos == source) THEN
353 12312 : buf(:) = pw%array(I1, I2, :)
354 648 : CALL group%send(buf, dest, tag)
355 : END IF
356 1296 : IF (mepos == dest) CALL group%recv(buf, source, tag)
357 : END IF
358 :
359 2592 : IF (mepos == dest) THEN
360 24624 : DO I3 = L3, U3
361 : x = pw%pw_grid%dh(1, 1)*(I1 - L1) &
362 : + pw%pw_grid%dh(2, 1)*(I2 - L2) &
363 23328 : + pw%pw_grid%dh(3, 1)*(I3 - L3)
364 : y = pw%pw_grid%dh(1, 2)*(I1 - L1) &
365 : + pw%pw_grid%dh(2, 2)*(I2 - L2) &
366 23328 : + pw%pw_grid%dh(3, 2)*(I3 - L3)
367 : z = pw%pw_grid%dh(1, 3)*(I1 - L1) &
368 : + pw%pw_grid%dh(2, 3)*(I2 - L2) &
369 23328 : + pw%pw_grid%dh(3, 3)*(I3 - L3)
370 24624 : IF (write_value) THEN
371 11664 : WRITE (unit_nr, '(3F16.8,2ES20.10E3)') x, y, z, dvol, buf(I3)
372 : ELSE
373 11664 : WRITE (unit_nr, '(3F16.8)') x, y, z
374 : END IF
375 : END DO
376 : END IF
377 :
378 2736 : CALL group%sync()
379 : END DO
380 : END DO
381 :
382 8 : DEALLOCATE (buf)
383 :
384 8 : IF (mepos == dest) CALL close_file(unit_number=unit_nr)
385 :
386 8 : CALL timestop(handle)
387 :
388 8 : END SUBROUTINE fde_write_pw
389 :
390 : ! **************************************************************************************************
391 : !> \brief Read the density in x y z weight value format
392 : !> \param pw ...
393 : !> \param filename ...
394 : ! **************************************************************************************************
395 4 : SUBROUTINE fde_read_pw(pw, filename)
396 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
397 : CHARACTER(len=*), INTENT(IN) :: filename
398 :
399 : CHARACTER(len=*), PARAMETER :: routineN = 'fde_read_pw'
400 :
401 : INTEGER :: dest, handle, I1, I2, I3, ip, L1, L2, &
402 : L3, mepos, npoints, num_pe, source, &
403 : tag, U1, U2, U3, unit_nr
404 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buf
405 : TYPE(mp_comm_type) :: group
406 :
407 4 : CALL timeset(routineN, handle)
408 :
409 4 : L1 = pw%pw_grid%bounds(1, 1)
410 4 : U1 = pw%pw_grid%bounds(2, 1)
411 4 : L2 = pw%pw_grid%bounds(1, 2)
412 4 : U2 = pw%pw_grid%bounds(2, 2)
413 4 : L3 = pw%pw_grid%bounds(1, 3)
414 4 : U3 = pw%pw_grid%bounds(2, 3)
415 :
416 4 : group = pw%pw_grid%para%group
417 4 : mepos = pw%pw_grid%para%group%mepos
418 4 : num_pe = pw%pw_grid%para%group%num_pe
419 4 : tag = 1
420 4 : dest = 0
421 :
422 4 : unit_nr = -1
423 4 : IF (mepos == dest) THEN
424 : CALL open_file(file_name=TRIM(filename), file_status="OLD", file_form="FORMATTED", &
425 2 : file_action="READ", unit_number=unit_nr)
426 2 : READ (unit_nr, *) npoints
427 2 : IF (npoints /= pw%pw_grid%ngpts) THEN
428 0 : CPABORT("Number of grid points in file does not match the pw grid")
429 : END IF
430 : END IF
431 :
432 12 : ALLOCATE (buf(L3:U3))
433 :
434 76 : DO I1 = L1, U1
435 1372 : DO I2 = L2, U2
436 1296 : IF (pw%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
437 1296 : source = dest
438 3888 : DO ip = 0, num_pe - 1
439 : IF (pw%pw_grid%para%bo(1, 1, ip, 1) <= I1 - L1 + 1 &
440 : .AND. pw%pw_grid%para%bo(2, 1, ip, 1) >= I1 - L1 + 1 &
441 : .AND. pw%pw_grid%para%bo(1, 2, ip, 1) <= I2 - L2 + 1 &
442 3888 : .AND. pw%pw_grid%para%bo(2, 2, ip, 1) >= I2 - L2 + 1) THEN
443 1296 : source = ip
444 : END IF
445 : END DO
446 : ELSE
447 0 : source = dest
448 : END IF
449 :
450 1296 : IF (mepos == dest) THEN
451 12312 : DO I3 = L3, U3
452 12312 : READ (unit_nr, *) buf(I3)
453 : END DO
454 : END IF
455 :
456 1296 : IF (source == dest) THEN
457 6480 : IF (mepos == source) pw%array(I1, I2, :) = buf(:)
458 : ELSE
459 648 : IF (mepos == dest) CALL group%send(buf, source, tag)
460 648 : IF (mepos == source) THEN
461 324 : CALL group%recv(buf, dest, tag)
462 6156 : pw%array(I1, I2, :) = buf(:)
463 : END IF
464 : END IF
465 :
466 1368 : CALL group%sync()
467 : END DO
468 : END DO
469 :
470 4 : DEALLOCATE (buf)
471 :
472 4 : IF (mepos == dest) CALL close_file(unit_number=unit_nr)
473 :
474 4 : CALL timestop(handle)
475 :
476 4 : END SUBROUTINE fde_read_pw
477 :
478 : ! **************************************************************************************************
479 : !> \brief Read scalar results returned by the subsystem solver, one value per line, in this order:
480 : !> E_sub : embedded subsystem total energy
481 : !> E_emb : embedding energy
482 : !> E_nuc : nuclear repulsion energy
483 : !> n_e : number of active-subsystem electrons
484 : !> n_a : number of nuclei of subsystem A
485 : !> Followed by n_a lines of
486 : !> Z_I x_I y_I z_I (charge and coordinates in the CP2K frame)
487 : !> \param filename ...
488 : !> \param para_env ...
489 : !> \param e_sub ...
490 : !> \param e_emb ...
491 : !> \param e_nuc ...
492 : !> \param n_e ...
493 : !> \param n_a ...
494 : !> \param z_a ...
495 : !> \param r_a ...
496 : ! **************************************************************************************************
497 4 : SUBROUTINE fde_read_results(filename, para_env, e_sub, e_emb, e_nuc, n_e, n_a, z_a, r_a)
498 : CHARACTER(len=*), INTENT(IN) :: filename
499 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
500 : REAL(KIND=dp), INTENT(OUT) :: e_sub, e_emb, e_nuc, n_e
501 : INTEGER, INTENT(OUT) :: n_a
502 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
503 : INTENT(OUT) :: z_a
504 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
505 : INTENT(OUT) :: r_a
506 :
507 : INTEGER :: iatom, unit_nr
508 : REAL(KIND=dp), DIMENSION(4) :: vals
509 :
510 4 : vals = 0.0_dp
511 4 : n_a = 0
512 4 : IF (para_env%is_source()) THEN
513 : CALL open_file(file_name=TRIM(filename), file_status="OLD", file_form="FORMATTED", &
514 2 : file_action="READ", unit_number=unit_nr)
515 2 : READ (unit_nr, *) vals(1)
516 2 : READ (unit_nr, *) vals(2)
517 2 : READ (unit_nr, *) vals(3)
518 2 : READ (unit_nr, *) vals(4)
519 2 : READ (unit_nr, *) n_a
520 : END IF
521 4 : CALL para_env%bcast(vals)
522 4 : CALL para_env%bcast(n_a)
523 :
524 20 : ALLOCATE (z_a(n_a), r_a(3, n_a))
525 4 : z_a = 0.0_dp
526 4 : r_a = 0.0_dp
527 4 : IF (para_env%is_source()) THEN
528 4 : DO iatom = 1, n_a
529 4 : READ (unit_nr, *) z_a(iatom), r_a(1, iatom), r_a(2, iatom), r_a(3, iatom)
530 : END DO
531 2 : CALL close_file(unit_number=unit_nr)
532 : END IF
533 4 : CALL para_env%bcast(z_a)
534 4 : CALL para_env%bcast(r_a)
535 :
536 4 : e_sub = vals(1)
537 4 : e_emb = vals(2)
538 4 : e_nuc = vals(3)
539 4 : n_e = vals(4)
540 :
541 4 : END SUBROUTINE fde_read_results
542 :
543 : END MODULE fde_methods
|
