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 : ! **************************************************************************************************
9 : !> \brief Routines to handle an external electrostatic field
10 : !> The external field can be generic and is provided by user input
11 : ! **************************************************************************************************
12 : MODULE qs_external_potential
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE cell_types, ONLY: cell_type,&
16 : pbc
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_log_handling, ONLY: cp_to_string
19 : USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw
20 : USE force_fields_util, ONLY: get_generic_info
21 : USE fparser, ONLY: evalf,&
22 : evalfd,&
23 : finalizef,&
24 : initf,&
25 : parsef
26 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
27 : section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: default_path_length,&
30 : default_string_length,&
31 : dp,&
32 : int_8
33 : USE maxwell_solver_interface, ONLY: maxwell_solver
34 : USE message_passing, ONLY: mp_comm_type
35 : USE particle_types, ONLY: particle_type
36 : USE pw_grid_types, ONLY: PW_MODE_LOCAL
37 : USE pw_methods, ONLY: pw_zero
38 : USE pw_types, ONLY: pw_r3d_rs_type
39 : USE qs_energy_types, ONLY: qs_energy_type
40 : USE qs_environment_types, ONLY: get_qs_env,&
41 : qs_environment_type
42 : USE qs_force_types, ONLY: qs_force_type
43 : USE qs_kind_types, ONLY: get_qs_kind,&
44 : qs_kind_type
45 : USE string_utilities, ONLY: compress
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 :
50 : PRIVATE
51 :
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_external_potential'
53 :
54 : ! *** Public subroutines ***
55 : PUBLIC :: external_e_potential, &
56 : external_c_potential, &
57 : interpolate_external_potential
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Computes the external potential on the grid
63 : !> \param qs_env ...
64 : !> \date 12.2009
65 : !> \author Teodoro Laino [tlaino]
66 : ! **************************************************************************************************
67 20384 : SUBROUTINE external_e_potential(qs_env)
68 :
69 : TYPE(qs_environment_type), POINTER :: qs_env
70 :
71 : CHARACTER(len=*), PARAMETER :: routineN = 'external_e_potential'
72 :
73 : INTEGER :: handle, i, j, k
74 : INTEGER(kind=int_8) :: npoints
75 : INTEGER, DIMENSION(2, 3) :: bo_global, bo_local
76 : REAL(kind=dp) :: dvol, scaling_factor
77 20384 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: efunc, grid_p_i, grid_p_j, grid_p_k
78 20384 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: grid_p
79 : REAL(kind=dp), DIMENSION(3) :: dr
80 : TYPE(dft_control_type), POINTER :: dft_control
81 : TYPE(pw_r3d_rs_type), POINTER :: v_ee
82 : TYPE(section_vals_type), POINTER :: ext_pot_section, input
83 :
84 20384 : CALL timeset(routineN, handle)
85 20384 : CALL get_qs_env(qs_env, dft_control=dft_control)
86 20384 : IF (dft_control%apply_external_potential) THEN
87 88 : IF (dft_control%eval_external_potential) THEN
88 16 : CALL get_qs_env(qs_env, vee=v_ee)
89 16 : IF (dft_control%expot_control%maxwell_solver) THEN
90 0 : scaling_factor = dft_control%expot_control%scaling_factor
91 : CALL maxwell_solver(dft_control%maxwell_control, v_ee, &
92 : qs_env%sim_step, qs_env%sim_time, &
93 0 : scaling_factor)
94 0 : dft_control%eval_external_potential = .FALSE.
95 16 : ELSE IF (dft_control%expot_control%read_from_cube) THEN
96 2 : scaling_factor = dft_control%expot_control%scaling_factor
97 2 : CALL cp_cube_to_pw(v_ee, 'pot.cube', scaling_factor)
98 2 : dft_control%eval_external_potential = .FALSE.
99 : ELSE
100 14 : CALL get_qs_env(qs_env, input=input)
101 14 : ext_pot_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_POTENTIAL")
102 :
103 56 : dr = v_ee%pw_grid%dr
104 14 : dvol = v_ee%pw_grid%dvol
105 14 : CALL pw_zero(v_ee)
106 :
107 140 : bo_local = v_ee%pw_grid%bounds_local
108 140 : bo_global = v_ee%pw_grid%bounds
109 :
110 : npoints = INT(bo_local(2, 1) - bo_local(1, 1) + 1, kind=int_8)* &
111 : INT(bo_local(2, 2) - bo_local(1, 2) + 1, kind=int_8)* &
112 14 : INT(bo_local(2, 3) - bo_local(1, 3) + 1, kind=int_8)
113 42 : ALLOCATE (efunc(npoints))
114 42 : ALLOCATE (grid_p(3, npoints))
115 42 : ALLOCATE (grid_p_i(bo_local(1, 1):bo_local(2, 1)))
116 42 : ALLOCATE (grid_p_j(bo_local(1, 2):bo_local(2, 2)))
117 42 : ALLOCATE (grid_p_k(bo_local(1, 3):bo_local(2, 3)))
118 :
119 378 : DO i = bo_local(1, 1), bo_local(2, 1)
120 378 : grid_p_i(i) = (i - bo_global(1, 1))*dr(1)
121 : END DO
122 742 : DO j = bo_local(1, 2), bo_local(2, 2)
123 742 : grid_p_j(j) = (j - bo_global(1, 2))*dr(2)
124 : END DO
125 742 : DO k = bo_local(1, 3), bo_local(2, 3)
126 742 : grid_p_k(k) = (k - bo_global(1, 3))*dr(3)
127 : END DO
128 :
129 : npoints = 0
130 742 : DO k = bo_local(1, 3), bo_local(2, 3)
131 38934 : DO j = bo_local(1, 2), bo_local(2, 2)
132 1047704 : DO i = bo_local(1, 1), bo_local(2, 1)
133 1008784 : npoints = npoints + 1
134 1008784 : grid_p(1, npoints) = grid_p_i(i)
135 1008784 : grid_p(2, npoints) = grid_p_j(j)
136 1046976 : grid_p(3, npoints) = grid_p_k(k)
137 : END DO
138 : END DO
139 : END DO
140 :
141 14 : DEALLOCATE (grid_p_i, grid_p_j, grid_p_k)
142 :
143 14 : CALL get_external_potential(grid_p, ext_pot_section, func=efunc)
144 :
145 14 : npoints = 0
146 742 : DO k = bo_local(1, 3), bo_local(2, 3)
147 38934 : DO j = bo_local(1, 2), bo_local(2, 2)
148 1047704 : DO i = bo_local(1, 1), bo_local(2, 1)
149 1008784 : npoints = npoints + 1
150 1046976 : v_ee%array(i, j, k) = v_ee%array(i, j, k) + efunc(npoints)
151 : END DO
152 : END DO
153 : END DO
154 :
155 14 : DEALLOCATE (grid_p, efunc)
156 :
157 14 : dft_control%eval_external_potential = .FALSE.
158 : END IF
159 : END IF
160 : END IF
161 20384 : CALL timestop(handle)
162 20384 : END SUBROUTINE external_e_potential
163 :
164 : ! **************************************************************************************************
165 : !> \brief Computes the force and the energy due to the external potential on the cores
166 : !> \param qs_env ...
167 : !> \param calculate_forces ...
168 : !> \date 12.2009
169 : !> \author Teodoro Laino [tlaino]
170 : ! **************************************************************************************************
171 17378 : SUBROUTINE external_c_potential(qs_env, calculate_forces)
172 :
173 : TYPE(qs_environment_type), POINTER :: qs_env
174 : LOGICAL, OPTIONAL :: calculate_forces
175 :
176 : CHARACTER(len=*), PARAMETER :: routineN = 'external_c_potential'
177 :
178 : INTEGER :: atom_a, handle, iatom, ikind, natom, &
179 : nkind, nparticles
180 17378 : INTEGER, DIMENSION(:), POINTER :: list
181 : LOGICAL :: my_force, pot_on_grid
182 : REAL(KIND=dp) :: ee_core_ener, zeff
183 17378 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: efunc
184 17378 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfunc, r
185 17378 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
186 : TYPE(cell_type), POINTER :: cell
187 : TYPE(dft_control_type), POINTER :: dft_control
188 17378 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
189 : TYPE(pw_r3d_rs_type), POINTER :: v_ee
190 : TYPE(qs_energy_type), POINTER :: energy
191 17378 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
192 17378 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
193 : TYPE(section_vals_type), POINTER :: ext_pot_section, input
194 :
195 17378 : CALL timeset(routineN, handle)
196 17378 : NULLIFY (dft_control)
197 :
198 : CALL get_qs_env(qs_env=qs_env, &
199 : atomic_kind_set=atomic_kind_set, &
200 : qs_kind_set=qs_kind_set, &
201 : energy=energy, &
202 : particle_set=particle_set, &
203 : input=input, &
204 : cell=cell, &
205 17378 : dft_control=dft_control)
206 :
207 17378 : IF (dft_control%apply_external_potential .AND. .NOT. qs_env%mimic) THEN
208 : !ensure that external potential is loaded to grid
209 82 : IF (dft_control%eval_external_potential) THEN
210 0 : CALL external_e_potential(qs_env)
211 : END IF
212 82 : my_force = .FALSE.
213 82 : IF (PRESENT(calculate_forces)) my_force = calculate_forces
214 82 : ee_core_ener = 0.0_dp
215 82 : nkind = SIZE(atomic_kind_set)
216 :
217 : !check if external potential on grid has been loaded from a file instead of giving a function
218 82 : IF (dft_control%expot_control%read_from_cube .OR. &
219 : dft_control%expot_control%maxwell_solver) THEN
220 24 : CALL get_qs_env(qs_env, vee=v_ee)
221 24 : pot_on_grid = .TRUE.
222 : ELSE
223 58 : pot_on_grid = .FALSE.
224 58 : ext_pot_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_POTENTIAL")
225 : END IF
226 :
227 82 : nparticles = 0
228 242 : DO ikind = 1, SIZE(atomic_kind_set)
229 160 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
230 242 : nparticles = nparticles + MAX(natom, 0)
231 : END DO
232 :
233 246 : ALLOCATE (efunc(nparticles))
234 328 : ALLOCATE (dfunc(3, nparticles), r(3, nparticles))
235 :
236 82 : nparticles = 0
237 242 : DO ikind = 1, SIZE(atomic_kind_set)
238 160 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom)
239 :
240 480 : DO iatom = 1, natom
241 238 : atom_a = list(iatom)
242 238 : nparticles = nparticles + 1
243 : !pbc returns r(i) in range [-cell%hmat(i,i)/2, cell%hmat(i,i)/2]
244 : !for periodic dimensions (assuming the cell is orthorombic).
245 : !This is not consistent with the potential on grid, where r(i) is
246 : !in range [0, cell%hmat(i,i)]
247 : !Use new pbc function with switch positive_range=.TRUE.
248 398 : r(:, nparticles) = pbc(particle_set(atom_a)%r(:), cell, positive_range=.TRUE.)
249 : END DO
250 : END DO
251 :
252 : !if potential is on grid, interpolate the value at r,
253 : !otherwise evaluate the given function
254 82 : IF (pot_on_grid) THEN
255 96 : DO iatom = 1, nparticles
256 : CALL interpolate_external_potential(r(:, iatom), v_ee, func=efunc(iatom), &
257 96 : dfunc=dfunc(:, iatom), calc_derivatives=my_force)
258 : END DO
259 : ELSE
260 58 : CALL get_external_potential(r, ext_pot_section, func=efunc, dfunc=dfunc, calc_derivatives=my_force)
261 : END IF
262 :
263 82 : IF (my_force) THEN
264 40 : CALL get_qs_env(qs_env=qs_env, force=force)
265 : END IF
266 :
267 82 : nparticles = 0
268 242 : DO ikind = 1, SIZE(atomic_kind_set)
269 160 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
270 160 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
271 :
272 480 : DO iatom = 1, natom
273 238 : nparticles = nparticles + 1
274 :
275 238 : ee_core_ener = ee_core_ener + zeff*efunc(nparticles)
276 398 : IF (my_force) THEN
277 464 : force(ikind)%eev(1:3, iatom) = dfunc(1:3, nparticles)*zeff
278 : END IF
279 : END DO
280 : END DO
281 82 : energy%ee_core = ee_core_ener
282 :
283 82 : DEALLOCATE (dfunc, r)
284 82 : DEALLOCATE (efunc)
285 : END IF
286 17378 : CALL timestop(handle)
287 34756 : END SUBROUTINE external_c_potential
288 :
289 : ! **************************************************************************************************
290 : !> \brief Low level function for computing the potential and the derivatives
291 : !> \param r position in realspace for each grid-point
292 : !> \param ext_pot_section ...
293 : !> \param func external potential at r
294 : !> \param dfunc derivative of the external potential at r
295 : !> \param calc_derivatives Whether to calculate dfunc
296 : !> \date 12.2009
297 : !> \par History
298 : !> 12.2009 created [tlaino]
299 : !> 11.2014 reading external cube file added [Juha Ritala & Matt Watkins]
300 : !> \author Teodoro Laino [tlaino]
301 : ! **************************************************************************************************
302 72 : SUBROUTINE get_external_potential(r, ext_pot_section, func, dfunc, calc_derivatives)
303 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
304 : TYPE(section_vals_type), POINTER :: ext_pot_section
305 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: func
306 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
307 : OPTIONAL :: dfunc
308 : LOGICAL, INTENT(IN), OPTIONAL :: calc_derivatives
309 :
310 : CHARACTER(len=*), PARAMETER :: routineN = 'get_external_potential'
311 :
312 : CHARACTER(LEN=default_path_length) :: coupling_function
313 : CHARACTER(LEN=default_string_length) :: def_error, this_error
314 : CHARACTER(LEN=default_string_length), &
315 72 : DIMENSION(:), POINTER :: my_par
316 : INTEGER :: handle, j
317 : INTEGER(kind=int_8) :: ipoint, npoints
318 : LOGICAL :: check, my_force
319 : REAL(KIND=dp) :: dedf, dx, err, lerr
320 72 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_val
321 :
322 72 : CALL timeset(routineN, handle)
323 72 : NULLIFY (my_par, my_val)
324 72 : my_force = .FALSE.
325 72 : IF (PRESENT(calc_derivatives)) my_force = calc_derivatives
326 72 : check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives)
327 72 : CPASSERT(check)
328 72 : CALL section_vals_val_get(ext_pot_section, "DX", r_val=dx)
329 72 : CALL section_vals_val_get(ext_pot_section, "ERROR_LIMIT", r_val=lerr)
330 : CALL get_generic_info(ext_pot_section, "FUNCTION", coupling_function, my_par, my_val, &
331 288 : input_variables=["X", "Y", "Z"], i_rep_sec=1)
332 72 : CALL initf(1)
333 72 : CALL parsef(1, TRIM(coupling_function), my_par)
334 :
335 72 : npoints = SIZE(r, 2, kind=int_8)
336 :
337 1009022 : DO ipoint = 1, npoints
338 1008950 : my_val(1) = r(1, ipoint)
339 1008950 : my_val(2) = r(2, ipoint)
340 1008950 : my_val(3) = r(3, ipoint)
341 :
342 1008950 : IF (PRESENT(func)) func(ipoint) = evalf(1, my_val)
343 1009022 : IF (my_force) THEN
344 320 : DO j = 1, 3
345 240 : dedf = evalfd(1, j, my_val, dx, err)
346 240 : IF (ABS(err) > lerr) THEN
347 0 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
348 0 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
349 0 : CALL compress(this_error, .TRUE.)
350 0 : CALL compress(def_error, .TRUE.)
351 : CALL cp_warn(__LOCATION__, &
352 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
353 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
354 0 : TRIM(def_error)//' .')
355 : END IF
356 320 : dfunc(j, ipoint) = dedf
357 : END DO
358 : END IF
359 : END DO
360 72 : DEALLOCATE (my_par)
361 72 : DEALLOCATE (my_val)
362 72 : CALL finalizef()
363 72 : CALL timestop(handle)
364 72 : END SUBROUTINE get_external_potential
365 :
366 : ! **************************************************************************************************
367 : !> \brief subroutine that interpolates the value of the external
368 : !> potential at position r based on the values on the realspace grid
369 : !> \param r ...
370 : !> \param grid external potential pw grid, vee
371 : !> \param func value of vee at r
372 : !> \param dfunc derivatives of vee at r
373 : !> \param calc_derivatives calc dfunc
374 : ! **************************************************************************************************
375 76 : SUBROUTINE interpolate_external_potential(r, grid, func, dfunc, calc_derivatives)
376 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
377 : TYPE(pw_r3d_rs_type), POINTER :: grid
378 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: func, dfunc(3)
379 : LOGICAL, INTENT(IN), OPTIONAL :: calc_derivatives
380 :
381 : CHARACTER(len=*), PARAMETER :: routineN = 'interpolate_external_potential'
382 :
383 : INTEGER :: buffer_i, buffer_j, buffer_k, &
384 : data_source, fd_extra_point, handle, &
385 : i, i_pbc, ip, j, j_pbc, k, k_pbc, &
386 : my_rank, num_pe, tag
387 : INTEGER, DIMENSION(3) :: lbounds, lbounds_local, lower_inds, &
388 : ubounds, ubounds_local, upper_inds
389 : LOGICAL :: check, my_force
390 76 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bcast_buffer
391 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: grid_buffer
392 76 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dgrid
393 : REAL(KIND=dp), DIMENSION(3) :: dr, subgrid_origin
394 : TYPE(mp_comm_type) :: gid
395 :
396 76 : CALL timeset(routineN, handle)
397 76 : my_force = .FALSE.
398 76 : IF (PRESENT(calc_derivatives)) my_force = calc_derivatives
399 76 : check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives)
400 76 : CPASSERT(check)
401 :
402 76 : IF (my_force) THEN
403 36 : ALLOCATE (grid_buffer(0:3, 0:3, 0:3))
404 36 : ALLOCATE (bcast_buffer(0:3))
405 36 : ALLOCATE (dgrid(1:2, 1:2, 1:2, 3))
406 36 : fd_extra_point = 1
407 : ELSE
408 40 : ALLOCATE (grid_buffer(1:2, 1:2, 1:2))
409 40 : ALLOCATE (bcast_buffer(1:2))
410 40 : fd_extra_point = 0
411 : END IF
412 :
413 : ! The values of external potential on grid are distributed among the
414 : ! processes, so first we have to gather them up
415 76 : gid = grid%pw_grid%para%group
416 76 : my_rank = grid%pw_grid%para%group%mepos
417 76 : num_pe = grid%pw_grid%para%group%num_pe
418 76 : tag = 1
419 :
420 304 : dr = grid%pw_grid%dr
421 304 : lbounds = grid%pw_grid%bounds(1, :)
422 304 : ubounds = grid%pw_grid%bounds(2, :)
423 304 : lbounds_local = grid%pw_grid%bounds_local(1, :)
424 304 : ubounds_local = grid%pw_grid%bounds_local(2, :)
425 :
426 : ! Determine the indices of grid points that are needed
427 304 : lower_inds = lbounds + FLOOR(r/dr) - fd_extra_point
428 304 : upper_inds = lower_inds + 1 + 2*fd_extra_point
429 :
430 300 : DO i = lower_inds(1), upper_inds(1)
431 : ! If index is out of global bounds, assume periodic boundary conditions
432 224 : i_pbc = pbc_index(i, lbounds(1), ubounds(1))
433 224 : buffer_i = i - lower_inds(1) + 1 - fd_extra_point
434 1036 : DO j = lower_inds(2), upper_inds(2)
435 736 : j_pbc = pbc_index(j, lbounds(2), ubounds(2))
436 736 : buffer_j = j - lower_inds(2) + 1 - fd_extra_point
437 :
438 : ! Find the process that has the data for indices i_pbc and j_pbc
439 : ! and store the data to bcast_buffer. Assuming that each process has full z data
440 960 : IF (grid%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
441 1200 : DO ip = 0, num_pe - 1
442 : IF (grid%pw_grid%para%bo(1, 1, ip, 1) <= i_pbc - lbounds(1) + 1 .AND. &
443 : grid%pw_grid%para%bo(2, 1, ip, 1) >= i_pbc - lbounds(1) + 1 .AND. &
444 1200 : grid%pw_grid%para%bo(1, 2, ip, 1) <= j_pbc - lbounds(2) + 1 .AND. &
445 0 : grid%pw_grid%para%bo(2, 2, ip, 1) >= j_pbc - lbounds(2) + 1) THEN
446 736 : data_source = ip
447 736 : EXIT
448 : END IF
449 : END DO
450 736 : IF (my_rank == data_source) THEN
451 368 : IF (lower_inds(3) >= lbounds(3) .AND. upper_inds(3) <= ubounds(3)) THEN
452 : bcast_buffer(:) = &
453 1680 : grid%array(i_pbc, j_pbc, lower_inds(3):upper_inds(3))
454 : ELSE
455 0 : DO k = lower_inds(3), upper_inds(3)
456 0 : k_pbc = pbc_index(k, lbounds(3), ubounds(3))
457 0 : buffer_k = k - lower_inds(3) + 1 - fd_extra_point
458 : bcast_buffer(buffer_k) = &
459 0 : grid%array(i_pbc, j_pbc, k_pbc)
460 : END DO
461 : END IF
462 : END IF
463 : ! data_source sends data to everyone
464 736 : CALL gid%bcast(bcast_buffer, data_source)
465 3360 : grid_buffer(buffer_i, buffer_j, :) = bcast_buffer
466 : ELSE
467 0 : grid_buffer(buffer_i, buffer_j, :) = grid%array(i_pbc, j_pbc, lower_inds(3):upper_inds(3))
468 : END IF
469 : END DO
470 : END DO
471 :
472 : ! Now that all the processes have local external potential data around r,
473 : ! interpolate the value at r
474 304 : subgrid_origin = (lower_inds - lbounds + fd_extra_point)*dr
475 76 : func = trilinear_interpolation(r, grid_buffer(1:2, 1:2, 1:2), subgrid_origin, dr)
476 :
477 : ! If the derivative of the potential is needed, approximate the derivative at grid
478 : ! points using finite differences, and then interpolate the value at r
479 76 : IF (my_force) THEN
480 36 : CALL d_finite_difference(grid_buffer, dr, dgrid)
481 144 : DO i = 1, 3
482 144 : dfunc(i) = trilinear_interpolation(r, dgrid(:, :, :, i), subgrid_origin, dr)
483 : END DO
484 36 : DEALLOCATE (dgrid)
485 : END IF
486 :
487 76 : DEALLOCATE (grid_buffer)
488 76 : CALL timestop(handle)
489 152 : END SUBROUTINE interpolate_external_potential
490 :
491 : ! **************************************************************************************************
492 : !> \brief subroutine that uses finite differences to approximate the partial
493 : !> derivatives of the potential based on the given values on grid
494 : !> \param grid tiny bit of external potential vee
495 : !> \param dr step size for finite difference
496 : !> \param dgrid derivatives of grid
497 : ! **************************************************************************************************
498 36 : PURE SUBROUTINE d_finite_difference(grid, dr, dgrid)
499 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(IN) :: grid
500 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dr
501 : REAL(KIND=dp), DIMENSION(1:, 1:, 1:, :), &
502 : INTENT(OUT) :: dgrid
503 :
504 : INTEGER :: i, j, k
505 :
506 108 : DO i = 1, SIZE(dgrid, 1)
507 252 : DO j = 1, SIZE(dgrid, 2)
508 504 : DO k = 1, SIZE(dgrid, 3)
509 288 : dgrid(i, j, k, 1) = 0.5_dp*(grid(i + 1, j, k) - grid(i - 1, j, k))/dr(1)
510 288 : dgrid(i, j, k, 2) = 0.5_dp*(grid(i, j + 1, k) - grid(i, j - 1, k))/dr(2)
511 432 : dgrid(i, j, k, 3) = 0.5_dp*(grid(i, j, k + 1) - grid(i, j, k - 1))/dr(3)
512 : END DO
513 : END DO
514 : END DO
515 36 : END SUBROUTINE d_finite_difference
516 :
517 : ! **************************************************************************************************
518 : !> \brief trilinear interpolation function that interpolates value at r based
519 : !> on 2x2x2 grid points around r in subgrid
520 : !> \param r where to interpolate to
521 : !> \param subgrid part of external potential on a grid
522 : !> \param origin center of grid
523 : !> \param dr step size
524 : !> \return interpolated value of external potential
525 : ! **************************************************************************************************
526 184 : PURE FUNCTION trilinear_interpolation(r, subgrid, origin, dr) RESULT(value_at_r)
527 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
528 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: subgrid
529 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: origin, dr
530 : REAL(KIND=dp) :: value_at_r
531 :
532 : REAL(KIND=dp), DIMENSION(3) :: norm_r, norm_r_rev
533 :
534 736 : norm_r = (r - origin)/dr
535 736 : norm_r_rev = 1 - norm_r
536 : value_at_r = subgrid(1, 1, 1)*PRODUCT(norm_r_rev) + &
537 : subgrid(2, 1, 1)*norm_r(1)*norm_r_rev(2)*norm_r_rev(3) + &
538 : subgrid(1, 2, 1)*norm_r_rev(1)*norm_r(2)*norm_r_rev(3) + &
539 : subgrid(1, 1, 2)*norm_r_rev(1)*norm_r_rev(2)*norm_r(3) + &
540 : subgrid(1, 2, 2)*norm_r_rev(1)*norm_r(2)*norm_r(3) + &
541 : subgrid(2, 1, 2)*norm_r(1)*norm_r_rev(2)*norm_r(3) + &
542 : subgrid(2, 2, 1)*norm_r(1)*norm_r(2)*norm_r_rev(3) + &
543 1472 : subgrid(2, 2, 2)*PRODUCT(norm_r)
544 184 : END FUNCTION trilinear_interpolation
545 :
546 : ! **************************************************************************************************
547 : !> \brief get a correct value for possible out of bounds index using periodic
548 : !> boundary conditions
549 : !> \param i ...
550 : !> \param lowbound ...
551 : !> \param upbound ...
552 : !> \return ...
553 : ! **************************************************************************************************
554 960 : ELEMENTAL FUNCTION pbc_index(i, lowbound, upbound)
555 : INTEGER, INTENT(IN) :: i, lowbound, upbound
556 : INTEGER :: pbc_index
557 :
558 960 : IF (i < lowbound) THEN
559 0 : pbc_index = upbound + i - lowbound + 1
560 960 : ELSE IF (i > upbound) THEN
561 0 : pbc_index = lowbound + i - upbound - 1
562 : ELSE
563 : pbc_index = i
564 : END IF
565 960 : END FUNCTION pbc_index
566 :
567 : END MODULE qs_external_potential
|
