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