(git:0de0cc2)
qs_external_potential.F
Go to the documentation of this file.
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 ! **************************************************************************************************
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
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
21  USE fparser, ONLY: evalf,&
22  evalfd,&
23  finalizef,&
24  initf,&
25  parsef
27  section_vals_type,&
29  USE kinds, ONLY: default_path_length,&
31  dp,&
32  int_8
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, &
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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: efunc, grid_p_i, grid_p_j, grid_p_k
77  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  CALL timeset(routinen, handle)
84  CALL get_qs_env(qs_env, dft_control=dft_control)
85  IF (dft_control%apply_external_potential) THEN
86  IF (dft_control%eval_external_potential) THEN
87  CALL get_qs_env(qs_env, vee=v_ee)
88  IF (dft_control%expot_control%maxwell_solver) THEN
89  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  scaling_factor)
93  dft_control%eval_external_potential = .false.
94  ELSEIF (dft_control%expot_control%read_from_cube) THEN
95  scaling_factor = dft_control%expot_control%scaling_factor
96  CALL cp_cube_to_pw(v_ee, 'pot.cube', scaling_factor)
97  dft_control%eval_external_potential = .false.
98  ELSE
99  CALL get_qs_env(qs_env, input=input)
100  ext_pot_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_POTENTIAL")
101 
102  dr = v_ee%pw_grid%dr
103  dvol = v_ee%pw_grid%dvol
104  CALL pw_zero(v_ee)
105 
106  bo_local = v_ee%pw_grid%bounds_local
107  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  int(bo_local(2, 3) - bo_local(1, 3) + 1, kind=int_8)
112  ALLOCATE (efunc(npoints))
113  ALLOCATE (grid_p(3, npoints))
114  ALLOCATE (grid_p_i(bo_local(1, 1):bo_local(2, 1)))
115  ALLOCATE (grid_p_j(bo_local(1, 2):bo_local(2, 2)))
116  ALLOCATE (grid_p_k(bo_local(1, 3):bo_local(2, 3)))
117 
118  DO i = bo_local(1, 1), bo_local(2, 1)
119  grid_p_i(i) = (i - bo_global(1, 1))*dr(1)
120  END DO
121  DO j = bo_local(1, 2), bo_local(2, 2)
122  grid_p_j(j) = (j - bo_global(1, 2))*dr(2)
123  END DO
124  DO k = bo_local(1, 3), bo_local(2, 3)
125  grid_p_k(k) = (k - bo_global(1, 3))*dr(3)
126  END DO
127 
128  npoints = 0
129  DO k = bo_local(1, 3), bo_local(2, 3)
130  DO j = bo_local(1, 2), bo_local(2, 2)
131  DO i = bo_local(1, 1), bo_local(2, 1)
132  npoints = npoints + 1
133  grid_p(1, npoints) = grid_p_i(i)
134  grid_p(2, npoints) = grid_p_j(j)
135  grid_p(3, npoints) = grid_p_k(k)
136  END DO
137  END DO
138  END DO
139 
140  DEALLOCATE (grid_p_i, grid_p_j, grid_p_k)
141 
142  CALL get_external_potential(grid_p, ext_pot_section, func=efunc)
143 
144  npoints = 0
145  DO k = bo_local(1, 3), bo_local(2, 3)
146  DO j = bo_local(1, 2), bo_local(2, 2)
147  DO i = bo_local(1, 1), bo_local(2, 1)
148  npoints = npoints + 1
149  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  DEALLOCATE (grid_p, efunc)
155 
156  dft_control%eval_external_potential = .false.
157  END IF
158  END IF
159  END IF
160  CALL timestop(handle)
161  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  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  INTEGER, DIMENSION(:), POINTER :: list
180  LOGICAL :: my_force, pot_on_grid
181  REAL(kind=dp) :: ee_core_ener, zeff
182  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: efunc
183  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dfunc, r
184  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
185  TYPE(cell_type), POINTER :: cell
186  TYPE(dft_control_type), POINTER :: dft_control
187  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
188  TYPE(pw_r3d_rs_type), POINTER :: v_ee
189  TYPE(qs_energy_type), POINTER :: energy
190  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
191  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
192  TYPE(section_vals_type), POINTER :: ext_pot_section, input
193 
194  CALL timeset(routinen, handle)
195  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  dft_control=dft_control)
205 
206  IF (dft_control%apply_external_potential) THEN
207  !ensure that external potential is loaded to grid
208  IF (dft_control%eval_external_potential) THEN
209  CALL external_e_potential(qs_env)
210  END IF
211  my_force = .false.
212  IF (PRESENT(calculate_forces)) my_force = calculate_forces
213  ee_core_ener = 0.0_dp
214  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  IF (dft_control%expot_control%read_from_cube .OR. &
218  dft_control%expot_control%maxwell_solver) THEN
219  CALL get_qs_env(qs_env, vee=v_ee)
220  pot_on_grid = .true.
221  ELSE
222  pot_on_grid = .false.
223  ext_pot_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_POTENTIAL")
224  END IF
225 
226  nparticles = 0
227  DO ikind = 1, SIZE(atomic_kind_set)
228  CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
229  nparticles = nparticles + max(natom, 0)
230  END DO
231 
232  ALLOCATE (efunc(nparticles))
233  ALLOCATE (dfunc(3, nparticles), r(3, nparticles))
234 
235  nparticles = 0
236  DO ikind = 1, SIZE(atomic_kind_set)
237  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom)
238 
239  DO iatom = 1, natom
240  atom_a = list(iatom)
241  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  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  IF (pot_on_grid) THEN
254  DO iatom = 1, nparticles
255  CALL interpolate_external_potential(r(:, iatom), v_ee, func=efunc(iatom), &
256  dfunc=dfunc(:, iatom), calc_derivatives=my_force)
257  END DO
258  ELSE
259  CALL get_external_potential(r, ext_pot_section, func=efunc, dfunc=dfunc, calc_derivatives=my_force)
260  END IF
261 
262  IF (my_force) THEN
263  CALL get_qs_env(qs_env=qs_env, force=force)
264  END IF
265 
266  nparticles = 0
267  DO ikind = 1, SIZE(atomic_kind_set)
268  CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
269  CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
270 
271  DO iatom = 1, natom
272  nparticles = nparticles + 1
273 
274  ee_core_ener = ee_core_ener + zeff*efunc(nparticles)
275  IF (my_force) THEN
276  force(ikind)%eev(1:3, iatom) = dfunc(1:3, nparticles)*zeff
277  END IF
278  END DO
279  END DO
280  energy%ee_core = ee_core_ener
281 
282  DEALLOCATE (dfunc, r)
283  DEALLOCATE (efunc)
284  END IF
285  CALL timestop(handle)
286  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  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: my_val
320 
321  CALL timeset(routinen, handle)
322  NULLIFY (my_par, my_val)
323  my_force = .false.
324  IF (PRESENT(calc_derivatives)) my_force = calc_derivatives
325  check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives)
326  cpassert(check)
327  CALL section_vals_val_get(ext_pot_section, "DX", r_val=dx)
328  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  input_variables=(/"X", "Y", "Z"/), i_rep_sec=1)
331  CALL initf(1)
332  CALL parsef(1, trim(coupling_function), my_par)
333 
334  npoints = SIZE(r, 2, kind=int_8)
335 
336  DO ipoint = 1, npoints
337  my_val(1) = r(1, ipoint)
338  my_val(2) = r(2, ipoint)
339  my_val(3) = r(3, ipoint)
340 
341  IF (PRESENT(func)) func(ipoint) = evalf(1, my_val)
342  IF (my_force) THEN
343  DO j = 1, 3
344  dedf = evalfd(1, j, my_val, dx, err)
345  IF (abs(err) > lerr) THEN
346  WRITE (this_error, "(A,G12.6,A)") "(", err, ")"
347  WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")"
348  CALL compress(this_error, .true.)
349  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  trim(def_error)//' .')
354  END IF
355  dfunc(j, ipoint) = dedf
356  END DO
357  END IF
358  END DO
359  DEALLOCATE (my_par)
360  DEALLOCATE (my_val)
361  CALL finalizef()
362  CALL timestop(handle)
363  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: bcast_buffer
390  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: grid_buffer
391  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dgrid
392  REAL(kind=dp), DIMENSION(3) :: dr, subgrid_origin
393  TYPE(mp_comm_type) :: gid
394 
395  CALL timeset(routinen, handle)
396  my_force = .false.
397  IF (PRESENT(calc_derivatives)) my_force = calc_derivatives
398  check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives)
399  cpassert(check)
400 
401  IF (my_force) THEN
402  ALLOCATE (grid_buffer(0:3, 0:3, 0:3))
403  ALLOCATE (bcast_buffer(0:3))
404  ALLOCATE (dgrid(1:2, 1:2, 1:2, 3))
405  fd_extra_point = 1
406  ELSE
407  ALLOCATE (grid_buffer(1:2, 1:2, 1:2))
408  ALLOCATE (bcast_buffer(1:2))
409  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  gid = grid%pw_grid%para%group
415  my_rank = grid%pw_grid%para%my_pos
416  num_pe = grid%pw_grid%para%group_size
417  tag = 1
418 
419  dr = grid%pw_grid%dr
420  lbounds = grid%pw_grid%bounds(1, :)
421  ubounds = grid%pw_grid%bounds(2, :)
422  lbounds_local = grid%pw_grid%bounds_local(1, :)
423  ubounds_local = grid%pw_grid%bounds_local(2, :)
424 
425  ! Determine the indices of grid points that are needed
426  lower_inds = lbounds + floor(r/dr) - fd_extra_point
427  upper_inds = lower_inds + 1 + 2*fd_extra_point
428 
429  DO i = lower_inds(1), upper_inds(1)
430  ! If index is out of global bounds, assume periodic boundary conditions
431  i_pbc = pbc_index(i, lbounds(1), ubounds(1))
432  buffer_i = i - lower_inds(1) + 1 - fd_extra_point
433  DO j = lower_inds(2), upper_inds(2)
434  j_pbc = pbc_index(j, lbounds(2), ubounds(2))
435  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  IF (grid%pw_grid%para%mode .NE. pw_mode_local) THEN
440  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  grid%pw_grid%para%bo(1, 2, ip, 1) <= j_pbc - lbounds(2) + 1 .AND. &
444  grid%pw_grid%para%bo(2, 2, ip, 1) >= j_pbc - lbounds(2) + 1) THEN
445  data_source = ip
446  EXIT
447  END IF
448  END DO
449  IF (my_rank == data_source) THEN
450  IF (lower_inds(3) >= lbounds(3) .AND. upper_inds(3) <= ubounds(3)) THEN
451  bcast_buffer(:) = &
452  grid%array(i_pbc, j_pbc, lower_inds(3):upper_inds(3))
453  ELSE
454  DO k = lower_inds(3), upper_inds(3)
455  k_pbc = pbc_index(k, lbounds(3), ubounds(3))
456  buffer_k = k - lower_inds(3) + 1 - fd_extra_point
457  bcast_buffer(buffer_k) = &
458  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  CALL gid%bcast(bcast_buffer, data_source)
464  grid_buffer(buffer_i, buffer_j, :) = bcast_buffer
465  ELSE
466  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  subgrid_origin = (lower_inds - lbounds + fd_extra_point)*dr
474  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  IF (my_force) THEN
479  CALL d_finite_difference(grid_buffer, dr, dgrid)
480  DO i = 1, 3
481  dfunc(i) = trilinear_interpolation(r, dgrid(:, :, :, i), subgrid_origin, dr)
482  END DO
483  DEALLOCATE (dgrid)
484  END IF
485 
486  DEALLOCATE (grid_buffer)
487  CALL timestop(handle)
488  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  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  DO i = 1, SIZE(dgrid, 1)
506  DO j = 1, SIZE(dgrid, 2)
507  DO k = 1, SIZE(dgrid, 3)
508  dgrid(i, j, k, 1) = 0.5*(grid(i + 1, j, k) - grid(i - 1, j, k))/dr(1)
509  dgrid(i, j, k, 2) = 0.5*(grid(i, j + 1, k) - grid(i, j - 1, k))/dr(2)
510  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  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  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  norm_r = (r - origin)/dr
534  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  subgrid(2, 2, 2)*product(norm_r)
543  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  ELEMENTAL FUNCTION pbc_index(i, lowbound, upbound)
554  INTEGER, INTENT(IN) :: i, lowbound, upbound
555  INTEGER :: pbc_index
556 
557  IF (i < lowbound) THEN
558  pbc_index = upbound + i - lowbound + 1
559  ELSE IF (i > upbound) THEN
560  pbc_index = lowbound + i - upbound - 1
561  ELSE
562  pbc_index = i
563  END IF
564  END FUNCTION pbc_index
565 
566 END MODULE qs_external_potential
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_cube_to_pw(grid, filename, scaling, silent)
Thin wrapper around routine cube_to_pw.
subroutine, public get_generic_info(gen_section, func_name, xfunction, parameters, values, var_values, size_variables, i_rep_sec, input_variables)
Reads from the input structure all information for generic functions.
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition: fparser.F:17
real(rn) function, public evalf(i, Val)
...
Definition: fparser.F:180
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
Definition: fparser.F:976
subroutine, public finalizef()
...
Definition: fparser.F:101
subroutine, public initf(n)
...
Definition: fparser.F:130
subroutine, public parsef(i, FuncStr, Var)
Parse ith function string FuncStr and compile it into bytecode.
Definition: fparser.F:148
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public int_8
Definition: kinds.F:54
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
integer, parameter, public default_path_length
Definition: kinds.F:58
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
Interface to Maxwell equation solver.
subroutine, public maxwell_solver(maxwell_control, v_ee, sim_step, sim_time, scaling_factor)
Computes the external potential on the grid.
Interface to the message passing library MPI.
Define the data structure for the particle information.
integer, parameter, public pw_mode_local
Definition: pw_grid_types.F:29
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Routines to handle an external electrostatic field The external field can be generic and is provided ...
subroutine, public external_c_potential(qs_env, calculate_forces)
Computes the force and the energy due to the external potential on the cores.
subroutine, public external_e_potential(qs_env)
Computes the external potential on the grid.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.