Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 :
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 16427 : 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 16427 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: efunc, grid_p_i, grid_p_j, grid_p_k
77 16427 : 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 16427 : CALL timeset(routineN, handle)
84 16427 : CALL get_qs_env(qs_env, dft_control=dft_control)
85 16427 : 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 : ELSEIF (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 16427 : CALL timestop(handle)
161 16427 : 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 14647 : 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 14647 : INTEGER, DIMENSION(:), POINTER :: list
180 : LOGICAL :: my_force, pot_on_grid
181 : REAL(KIND=dp) :: ee_core_ener, zeff
182 14647 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: efunc
183 14647 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfunc, r
184 14647 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
185 : TYPE(cell_type), POINTER :: cell
186 : TYPE(dft_control_type), POINTER :: dft_control
187 14647 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
188 : TYPE(pw_r3d_rs_type), POINTER :: v_ee
189 : TYPE(qs_energy_type), POINTER :: energy
190 14647 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
191 14647 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
192 : TYPE(section_vals_type), POINTER :: ext_pot_section, input
193 :
194 14647 : CALL timeset(routineN, handle)
195 14647 : 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 14647 : dft_control=dft_control)
205 :
206 14647 : IF (dft_control%apply_external_potential) 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 410 : 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 14647 : CALL timestop(handle)
286 29294 : 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 : DIMENSION(:), POINTER :: 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), DIMENSION(:), POINTER :: my_val
320 :
321 72 : CALL timeset(routineN, handle)
322 72 : NULLIFY (my_par, my_val)
323 72 : my_force = .FALSE.
324 72 : IF (PRESENT(calc_derivatives)) my_force = calc_derivatives
325 72 : check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives)
326 72 : CPASSERT(check)
327 72 : CALL section_vals_val_get(ext_pot_section, "DX", r_val=dx)
328 72 : CALL section_vals_val_get(ext_pot_section, "ERROR_LIMIT", r_val=lerr)
329 : CALL get_generic_info(ext_pot_section, "FUNCTION", coupling_function, my_par, my_val, &
330 288 : input_variables=(/"X", "Y", "Z"/), i_rep_sec=1)
331 72 : CALL initf(1)
332 72 : CALL parsef(1, TRIM(coupling_function), my_par)
333 :
334 72 : npoints = SIZE(r, 2, kind=int_8)
335 :
336 1009022 : DO ipoint = 1, npoints
337 1008950 : my_val(1) = r(1, ipoint)
338 1008950 : my_val(2) = r(2, ipoint)
339 1008950 : my_val(3) = r(3, ipoint)
340 :
341 1008950 : IF (PRESENT(func)) func(ipoint) = evalf(1, my_val)
342 1009022 : IF (my_force) THEN
343 320 : DO j = 1, 3
344 240 : dedf = evalfd(1, j, my_val, dx, err)
345 240 : IF (ABS(err) > lerr) THEN
346 0 : WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
347 0 : WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
348 0 : CALL compress(this_error, .TRUE.)
349 0 : CALL compress(def_error, .TRUE.)
350 : CALL cp_warn(__LOCATION__, &
351 : 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// &
352 : ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// &
353 0 : TRIM(def_error)//' .')
354 : END IF
355 320 : dfunc(j, ipoint) = dedf
356 : END DO
357 : END IF
358 : END DO
359 72 : DEALLOCATE (my_par)
360 72 : DEALLOCATE (my_val)
361 72 : CALL finalizef()
362 72 : CALL timestop(handle)
363 72 : END SUBROUTINE get_external_potential
364 :
365 : ! **************************************************************************************************
366 : !> \brief subroutine that interpolates the value of the external
367 : !> potential at position r based on the values on the realspace grid
368 : !> \param r ...
369 : !> \param grid external potential pw grid, vee
370 : !> \param func value of vee at r
371 : !> \param dfunc derivatives of vee at r
372 : !> \param calc_derivatives calc dfunc
373 : ! **************************************************************************************************
374 72 : SUBROUTINE interpolate_external_potential(r, grid, func, dfunc, calc_derivatives)
375 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
376 : TYPE(pw_r3d_rs_type), POINTER :: grid
377 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: func, dfunc(3)
378 : LOGICAL, INTENT(IN), OPTIONAL :: calc_derivatives
379 :
380 : CHARACTER(len=*), PARAMETER :: routineN = 'interpolate_external_potential'
381 :
382 : INTEGER :: buffer_i, buffer_j, buffer_k, &
383 : data_source, fd_extra_point, handle, &
384 : i, i_pbc, ip, j, j_pbc, k, k_pbc, &
385 : my_rank, num_pe, tag
386 : INTEGER, DIMENSION(3) :: lbounds, lbounds_local, lower_inds, &
387 : ubounds, ubounds_local, upper_inds
388 : LOGICAL :: check, my_force
389 72 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bcast_buffer
390 72 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: grid_buffer
391 72 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dgrid
392 : REAL(KIND=dp), DIMENSION(3) :: dr, subgrid_origin
393 : TYPE(mp_comm_type) :: gid
394 :
395 72 : CALL timeset(routineN, handle)
396 72 : my_force = .FALSE.
397 72 : IF (PRESENT(calc_derivatives)) my_force = calc_derivatives
398 72 : check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives)
399 72 : CPASSERT(check)
400 :
401 72 : IF (my_force) THEN
402 36 : ALLOCATE (grid_buffer(0:3, 0:3, 0:3))
403 36 : ALLOCATE (bcast_buffer(0:3))
404 36 : ALLOCATE (dgrid(1:2, 1:2, 1:2, 3))
405 36 : fd_extra_point = 1
406 : ELSE
407 36 : ALLOCATE (grid_buffer(1:2, 1:2, 1:2))
408 36 : ALLOCATE (bcast_buffer(1:2))
409 36 : fd_extra_point = 0
410 : END IF
411 :
412 : ! The values of external potential on grid are distributed among the
413 : ! processes, so first we have to gather them up
414 72 : gid = grid%pw_grid%para%group
415 72 : my_rank = grid%pw_grid%para%group%mepos
416 72 : num_pe = grid%pw_grid%para%group%num_pe
417 72 : tag = 1
418 :
419 288 : dr = grid%pw_grid%dr
420 288 : lbounds = grid%pw_grid%bounds(1, :)
421 288 : ubounds = grid%pw_grid%bounds(2, :)
422 288 : lbounds_local = grid%pw_grid%bounds_local(1, :)
423 288 : ubounds_local = grid%pw_grid%bounds_local(2, :)
424 :
425 : ! Determine the indices of grid points that are needed
426 288 : lower_inds = lbounds + FLOOR(r/dr) - fd_extra_point
427 288 : upper_inds = lower_inds + 1 + 2*fd_extra_point
428 :
429 288 : DO i = lower_inds(1), upper_inds(1)
430 : ! If index is out of global bounds, assume periodic boundary conditions
431 216 : i_pbc = pbc_index(i, lbounds(1), ubounds(1))
432 216 : buffer_i = i - lower_inds(1) + 1 - fd_extra_point
433 1008 : DO j = lower_inds(2), upper_inds(2)
434 720 : j_pbc = pbc_index(j, lbounds(2), ubounds(2))
435 720 : buffer_j = j - lower_inds(2) + 1 - fd_extra_point
436 :
437 : ! Find the process that has the data for indices i_pbc and j_pbc
438 : ! and store the data to bcast_buffer. Assuming that each process has full z data
439 936 : IF (grid%pw_grid%para%mode .NE. PW_MODE_LOCAL) THEN
440 1176 : DO ip = 0, num_pe - 1
441 : IF (grid%pw_grid%para%bo(1, 1, ip, 1) <= i_pbc - lbounds(1) + 1 .AND. &
442 : grid%pw_grid%para%bo(2, 1, ip, 1) >= i_pbc - lbounds(1) + 1 .AND. &
443 1176 : grid%pw_grid%para%bo(1, 2, ip, 1) <= j_pbc - lbounds(2) + 1 .AND. &
444 0 : grid%pw_grid%para%bo(2, 2, ip, 1) >= j_pbc - lbounds(2) + 1) THEN
445 720 : data_source = ip
446 720 : EXIT
447 : END IF
448 : END DO
449 720 : IF (my_rank == data_source) THEN
450 360 : IF (lower_inds(3) >= lbounds(3) .AND. upper_inds(3) <= ubounds(3)) THEN
451 : bcast_buffer(:) = &
452 1656 : grid%array(i_pbc, j_pbc, lower_inds(3):upper_inds(3))
453 : ELSE
454 0 : DO k = lower_inds(3), upper_inds(3)
455 0 : k_pbc = pbc_index(k, lbounds(3), ubounds(3))
456 0 : buffer_k = k - lower_inds(3) + 1 - fd_extra_point
457 : bcast_buffer(buffer_k) = &
458 0 : grid%array(i_pbc, j_pbc, k_pbc)
459 : END DO
460 : END IF
461 : END IF
462 : ! data_source sends data to everyone
463 720 : CALL gid%bcast(bcast_buffer, data_source)
464 3312 : grid_buffer(buffer_i, buffer_j, :) = bcast_buffer
465 : ELSE
466 0 : grid_buffer(buffer_i, buffer_j, :) = grid%array(i_pbc, j_pbc, lower_inds(3):upper_inds(3))
467 : END IF
468 : END DO
469 : END DO
470 :
471 : ! Now that all the processes have local external potential data around r,
472 : ! interpolate the value at r
473 288 : subgrid_origin = (lower_inds - lbounds + fd_extra_point)*dr
474 72 : func = trilinear_interpolation(r, grid_buffer(1:2, 1:2, 1:2), subgrid_origin, dr)
475 :
476 : ! If the derivative of the potential is needed, approximate the derivative at grid
477 : ! points using finite differences, and then interpolate the value at r
478 72 : IF (my_force) THEN
479 36 : CALL d_finite_difference(grid_buffer, dr, dgrid)
480 144 : DO i = 1, 3
481 144 : dfunc(i) = trilinear_interpolation(r, dgrid(:, :, :, i), subgrid_origin, dr)
482 : END DO
483 36 : DEALLOCATE (dgrid)
484 : END IF
485 :
486 72 : DEALLOCATE (grid_buffer)
487 72 : CALL timestop(handle)
488 144 : END SUBROUTINE interpolate_external_potential
489 :
490 : ! **************************************************************************************************
491 : !> \brief subroutine that uses finite differences to approximate the partial
492 : !> derivatives of the potential based on the given values on grid
493 : !> \param grid tiny bit of external potential vee
494 : !> \param dr step size for finite difference
495 : !> \param dgrid derivatives of grid
496 : ! **************************************************************************************************
497 36 : PURE SUBROUTINE d_finite_difference(grid, dr, dgrid)
498 : REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(IN) :: grid
499 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dr
500 : REAL(KIND=dp), DIMENSION(1:, 1:, 1:, :), &
501 : INTENT(OUT) :: dgrid
502 :
503 : INTEGER :: i, j, k
504 :
505 108 : DO i = 1, SIZE(dgrid, 1)
506 252 : DO j = 1, SIZE(dgrid, 2)
507 504 : DO k = 1, SIZE(dgrid, 3)
508 288 : dgrid(i, j, k, 1) = 0.5*(grid(i + 1, j, k) - grid(i - 1, j, k))/dr(1)
509 288 : dgrid(i, j, k, 2) = 0.5*(grid(i, j + 1, k) - grid(i, j - 1, k))/dr(2)
510 432 : dgrid(i, j, k, 3) = 0.5*(grid(i, j, k + 1) - grid(i, j, k - 1))/dr(3)
511 : END DO
512 : END DO
513 : END DO
514 36 : END SUBROUTINE d_finite_difference
515 :
516 : ! **************************************************************************************************
517 : !> \brief trilinear interpolation function that interpolates value at r based
518 : !> on 2x2x2 grid points around r in subgrid
519 : !> \param r where to interpolate to
520 : !> \param subgrid part of external potential on a grid
521 : !> \param origin center of grid
522 : !> \param dr step size
523 : !> \return interpolated value of external potential
524 : ! **************************************************************************************************
525 180 : PURE FUNCTION trilinear_interpolation(r, subgrid, origin, dr) RESULT(value_at_r)
526 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
527 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: subgrid
528 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: origin, dr
529 : REAL(KIND=dp) :: value_at_r
530 :
531 : REAL(KIND=dp), DIMENSION(3) :: norm_r, norm_r_rev
532 :
533 720 : norm_r = (r - origin)/dr
534 720 : norm_r_rev = 1 - norm_r
535 : value_at_r = subgrid(1, 1, 1)*PRODUCT(norm_r_rev) + &
536 : subgrid(2, 1, 1)*norm_r(1)*norm_r_rev(2)*norm_r_rev(3) + &
537 : subgrid(1, 2, 1)*norm_r_rev(1)*norm_r(2)*norm_r_rev(3) + &
538 : subgrid(1, 1, 2)*norm_r_rev(1)*norm_r_rev(2)*norm_r(3) + &
539 : subgrid(1, 2, 2)*norm_r_rev(1)*norm_r(2)*norm_r(3) + &
540 : subgrid(2, 1, 2)*norm_r(1)*norm_r_rev(2)*norm_r(3) + &
541 : subgrid(2, 2, 1)*norm_r(1)*norm_r(2)*norm_r_rev(3) + &
542 1440 : subgrid(2, 2, 2)*PRODUCT(norm_r)
543 180 : END FUNCTION trilinear_interpolation
544 :
545 : ! **************************************************************************************************
546 : !> \brief get a correct value for possible out of bounds index using periodic
547 : !> boundary conditions
548 : !> \param i ...
549 : !> \param lowbound ...
550 : !> \param upbound ...
551 : !> \return ...
552 : ! **************************************************************************************************
553 936 : ELEMENTAL FUNCTION pbc_index(i, lowbound, upbound)
554 : INTEGER, INTENT(IN) :: i, lowbound, upbound
555 : INTEGER :: pbc_index
556 :
557 936 : IF (i < lowbound) THEN
558 0 : pbc_index = upbound + i - lowbound + 1
559 936 : ELSE IF (i > upbound) THEN
560 0 : pbc_index = lowbound + i - upbound - 1
561 : ELSE
562 : pbc_index = i
563 : END IF
564 936 : END FUNCTION pbc_index
565 :
566 : END MODULE qs_external_potential
|
