(git:34ef472)
qmmm_gpw_forces.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 compute energy and forces in a QM/MM calculation
10 !> \par History
11 !> 05.2004 created [tlaino]
12 !> \author Teodoro Laino
13 ! **************************************************************************************************
14 MODULE qmmm_gpw_forces
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_get_default_logger,&
19  cp_logger_type
20  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
21  cp_print_key_unit_nr
22  USE cp_spline_utils, ONLY: pw_restrict_s3,&
23  spline3_nopbc_interp,&
24  spline3_pbc_interp
25  USE cube_utils, ONLY: cube_info_type
26  USE input_constants, ONLY: do_par_atom,&
27  do_qmmm_coulomb,&
28  do_qmmm_gauss,&
29  do_qmmm_none,&
30  do_qmmm_pcharge,&
31  do_qmmm_swave
32  USE input_section_types, ONLY: section_vals_get_subs_vals,&
33  section_vals_type,&
34  section_vals_val_get
35  USE kinds, ONLY: dp
36  USE message_passing, ONLY: mp_comm_type,&
37  mp_para_env_type,&
38  mp_request_type
39  USE mm_collocate_potential, ONLY: collocate_gf_rspace_nopbc,&
40  integrate_gf_rspace_nopbc
41  USE particle_types, ONLY: particle_type
42  USE pw_env_types, ONLY: pw_env_get,&
43  pw_env_type
44  USE pw_methods, ONLY: pw_axpy,&
45  pw_integral_ab,&
46  pw_transfer,&
47  pw_zero
48  USE pw_pool_types, ONLY: pw_pool_p_type,&
49  pw_pool_type,&
50  pw_pools_create_pws,&
51  pw_pools_give_back_pws
52  USE pw_types, ONLY: pw_c1d_gs_type,&
53  pw_r3d_rs_type
54  USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,&
55  qmmm_gaussian_type
56  USE qmmm_gpw_energy, ONLY: qmmm_elec_with_gaussian,&
57  qmmm_elec_with_gaussian_lg,&
58  qmmm_elec_with_gaussian_lr
59  USE qmmm_se_forces, ONLY: deriv_se_qmmm_matrix
60  USE qmmm_tb_methods, ONLY: deriv_tb_qmmm_matrix,&
61  deriv_tb_qmmm_matrix_pc
62  USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
63  qmmm_per_pot_p_type,&
64  qmmm_per_pot_type,&
65  qmmm_pot_p_type,&
66  qmmm_pot_type
67  USE qmmm_util, ONLY: spherical_cutoff_factor
68  USE qs_energy_types, ONLY: qs_energy_type
69  USE qs_environment_types, ONLY: get_qs_env,&
70  qs_environment_type
71  USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
72  USE qs_rho_types, ONLY: qs_rho_get,&
73  qs_rho_type
74 #include "./base/base_uses.f90"
75 
76  IMPLICIT NONE
77 
78  PRIVATE
79  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
80  REAL(KIND=dp), PARAMETER, PRIVATE :: dx = 0.01_dp ! Debug Variables
81  REAL(KIND=dp), PARAMETER, PRIVATE :: maxerr = 10.0_dp ! Debug Variables
82  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gpw_forces'
83  PUBLIC :: qmmm_forces
84 
85 CONTAINS
86 
87 ! **************************************************************************************************
88 !> \brief General driver to Compute the contribution
89 !> to the forces due to the QM/MM potential
90 !> \param qs_env ...
91 !> \param qmmm_env ...
92 !> \param mm_particles ...
93 !> \param calc_force ...
94 !> \param mm_cell ...
95 !> \par History
96 !> 06.2004 created [tlaino]
97 !> \author Teodoro Laino
98 ! **************************************************************************************************
99  SUBROUTINE qmmm_forces(qs_env, qmmm_env, mm_particles, calc_force, mm_cell)
100  TYPE(qs_environment_type), POINTER :: qs_env
101  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
102  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
103  LOGICAL, INTENT(in), OPTIONAL :: calc_force
104  TYPE(cell_type), POINTER :: mm_cell
105 
106  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces'
107 
108  INTEGER :: handle, iatom, image_indmm, imm, indmm, &
109  ispin, iw
110  LOGICAL :: gapw, need_f, periodic
111  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges, &
112  forces_added_shells
113  TYPE(cp_logger_type), POINTER :: logger
114  TYPE(dft_control_type), POINTER :: dft_control
115  TYPE(mp_para_env_type), POINTER :: para_env
116  TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
117  TYPE(pw_env_type), POINTER :: pw_env
118  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
119  TYPE(pw_pool_type), POINTER :: auxbas_pool
120  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
121  TYPE(pw_r3d_rs_type), POINTER :: rho_tot_r, rho_tot_r2
122  TYPE(qs_energy_type), POINTER :: energy
123  TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
124  TYPE(qs_rho_type), POINTER :: rho
125  TYPE(section_vals_type), POINTER :: input_section, interp_section, &
126  print_section
127 
128  CALL timeset(routinen, handle)
129  need_f = .true.
130  periodic = qmmm_env%periodic
131  IF (PRESENT(calc_force)) need_f = calc_force
132  NULLIFY (dft_control, ks_qmmm_env_loc, rho, pw_env, rho_tot_r, energy, forces, &
133  forces_added_charges, input_section, rho0_s_gs, rho_r)
134  CALL get_qs_env(qs_env=qs_env, &
135  rho=rho, &
136  rho_core=rho_core, &
137  pw_env=pw_env, &
138  energy=energy, &
139  para_env=para_env, &
140  input=input_section, &
141  rho0_s_gs=rho0_s_gs, &
142  dft_control=dft_control)
143 
144  CALL qs_rho_get(rho, rho_r=rho_r)
145 
146  logger => cp_get_default_logger()
147  ks_qmmm_env_loc => qs_env%ks_qmmm_env
148  interp_section => section_vals_get_subs_vals(input_section, "QMMM%INTERPOLATOR")
149  print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
150  iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
151  extension=".qmmmLog")
152  gapw = dft_control%qs_control%gapw
153  ! If forces are required allocate these temporary arrays
154  IF (need_f) THEN
155  ALLOCATE (forces(3, qmmm_env%num_mm_atoms))
156  ALLOCATE (forces_added_charges(3, qmmm_env%added_charges%num_mm_atoms))
157  ALLOCATE (forces_added_shells(3, qmmm_env%added_shells%num_mm_atoms))
158  forces(:, :) = 0.0_dp
159  forces_added_charges(:, :) = 0.0_dp
160  forces_added_shells(:, :) = 0.0_dp
161  END IF
162  IF (dft_control%qs_control%semi_empirical) THEN
163  ! SEMIEMPIRICAL
164  SELECT CASE (qmmm_env%qmmm_coupl_type)
165  CASE (do_qmmm_coulomb)
166  CALL deriv_se_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
167  need_f, forces, forces_added_charges)
168  CASE (do_qmmm_pcharge)
169  cpabort("Point Charge QM/MM electrostatic coupling not yet implemented for SE.")
170  CASE (do_qmmm_gauss, do_qmmm_swave)
171  cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.")
172  CASE (do_qmmm_none)
173  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
174  "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
175  CASE DEFAULT
176  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
177  cpabort("")
178  END SELECT
179  ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
180  ! DFTB
181  SELECT CASE (qmmm_env%qmmm_coupl_type)
182  CASE (do_qmmm_none)
183  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
184  "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
185  CASE (do_qmmm_coulomb)
186  CALL deriv_tb_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
187  need_f, forces, forces_added_charges)
188  CASE (do_qmmm_pcharge)
189  CALL deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
190  need_f, forces, forces_added_charges)
191  CASE (do_qmmm_gauss, do_qmmm_swave)
192  cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for DFTB.")
193  CASE DEFAULT
194  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
195  cpabort("")
196  END SELECT
197  IF (need_f) THEN
198  forces(:, :) = forces(:, :)/real(para_env%num_pe, kind=dp)
199  forces_added_charges(:, :) = forces_added_charges(:, :)/real(para_env%num_pe, kind=dp)
200  END IF
201  ELSE
202  ! GPW/GAPW
203  CALL pw_env_get(pw_env=pw_env, &
204  pw_pools=pw_pools, &
205  auxbas_pw_pool=auxbas_pool)
206  NULLIFY (rho_tot_r)
207  ALLOCATE (rho_tot_r)
208  CALL auxbas_pool%create_pw(rho_tot_r)
209  ! IF GAPW the core charge is replaced by the compensation charge
210  IF (gapw) THEN
211  IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
212  CALL pw_transfer(rho_core, rho_tot_r)
213  energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
214  NULLIFY (rho_tot_r2)
215  ALLOCATE (rho_tot_r2)
216  CALL auxbas_pool%create_pw(rho_tot_r2)
217  CALL pw_transfer(rho0_s_gs, rho_tot_r2)
218  CALL pw_axpy(rho_tot_r2, rho_tot_r)
219  CALL auxbas_pool%give_back_pw(rho_tot_r2)
220  DEALLOCATE (rho_tot_r2)
221  ELSE
222  CALL pw_transfer(rho0_s_gs, rho_tot_r)
223  !
224  ! QM/MM Nuclear Electrostatic Potential already included through rho0
225  !
226  energy%qmmm_nu = 0.0_dp
227  END IF
228  ELSE
229  CALL pw_transfer(rho_core, rho_tot_r)
230  !
231  ! Computes the QM/MM Nuclear Electrostatic Potential
232  !
233  energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
234  END IF
235  IF (need_f) THEN
236  !
237  DO ispin = 1, SIZE(rho_r)
238  CALL pw_axpy(rho_r(ispin), rho_tot_r)
239  END DO
240  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Evaluating forces on MM atoms due to the:"
241  ! Electrostatic Interaction type...
242  SELECT CASE (qmmm_env%qmmm_coupl_type)
243  CASE (do_qmmm_coulomb)
244  cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
245  CASE (do_qmmm_pcharge)
246  cpabort("Point Charge QM/MM electrostatic coupling not yet implemented for GPW/GAPW.")
247  CASE (do_qmmm_gauss, do_qmmm_swave)
248  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
249  "- QM/MM Coupling computed collocating the Gaussian Potential Functions."
250  CALL qmmm_forces_with_gaussian(rho=rho_tot_r, &
251  qmmm_env=qmmm_env, &
252  mm_particles=mm_particles, &
253  aug_pools=qmmm_env%aug_pools, &
254  auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
255  coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
256  para_env=para_env, &
257  pw_pools=pw_pools, &
258  eps_mm_rspace=qmmm_env%eps_mm_rspace, &
259  cube_info=ks_qmmm_env_loc%cube_info, &
260  forces=forces, &
261  forces_added_charges=forces_added_charges, &
262  forces_added_shells=forces_added_shells, &
263  interp_section=interp_section, &
264  iw=iw, &
265  mm_cell=mm_cell)
266  CASE (do_qmmm_none)
267  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
268  "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
269  CASE DEFAULT
270  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "- Unknown Coupling..."
271  cpabort("")
272  END SELECT
273  END IF
274  END IF
275  ! Correct Total Energy adding the contribution of the QM/MM nuclear interaction
276  energy%total = energy%total + energy%qmmm_nu
277  ! Proceed if gradients are requested..
278  IF (need_f) THEN
279  !ikuo Temporary change to alleviate compiler problems on Intel with
280  !array dimension of 0
281  IF (qmmm_env%num_mm_atoms /= 0) CALL para_env%sum(forces)
282  IF (qmmm_env%added_charges%num_mm_atoms /= 0) CALL para_env%sum(forces_added_charges)
283  IF (qmmm_env%added_shells%num_mm_atoms /= 0) CALL para_env%sum(forces_added_shells)
284  ! Debug Forces
285  IF (debug_this_module) THEN
286  IF (dft_control%qs_control%semi_empirical .OR. &
287  dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
288  WRITE (iw, *) "NO DEBUG AVAILABLE in module"//trim(routinen)
289  ELSE
290  ! Print Out Forces
291  IF (iw > 0) THEN
292  DO imm = 1, SIZE(qmmm_env%mm_atom_index)
293  WRITE (iw, *) "ANALYTICAL FORCES:"
294  indmm = qmmm_env%mm_atom_index(imm)
295  WRITE (iw, '(I6,3F15.9)') indmm, forces(:, imm)
296  END DO
297  END IF
298  CALL qmmm_debug_forces(rho=rho_tot_r, &
299  qs_env=qs_env, &
300  qmmm_env=qmmm_env, &
301  analytical_forces=forces, &
302  mm_particles=mm_particles, &
303  mm_atom_index=qmmm_env%mm_atom_index, &
304  num_mm_atoms=qmmm_env%num_mm_atoms, &
305  interp_section=interp_section, &
306  mm_cell=mm_cell)
307  END IF
308  END IF
309  END IF
310  ! Give back rho_tot_t to auxbas_pool only for GPW/GAPW
311  IF ((.NOT. dft_control%qs_control%semi_empirical) .AND. &
312  (.NOT. dft_control%qs_control%dftb) .AND. (.NOT. dft_control%qs_control%xtb)) THEN
313  CALL auxbas_pool%give_back_pw(rho_tot_r)
314  DEALLOCATE (rho_tot_r)
315  END IF
316  IF (iw > 0) THEN
317  IF (.NOT. gapw) WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
318  "QM/MM Nuclear Electrostatic Potential :", energy%qmmm_nu
319  WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
320  "QMMM Total Energy (QM + QMMM electronic + QMMM nuclear):", energy%total
321  WRITE (iw, '(T2,"QMMM|",1X,A)') "MM energy NOT included in the above term!"// &
322  " Check for: FORCE_EVAL ( QMMM )"
323  WRITE (iw, '(T2,"QMMM|",1X,A)') "that includes both QM, QMMM and MM energy terms!"
324  END IF
325  IF (need_f) THEN
326  ! Transfer Forces
327  DO imm = 1, qmmm_env%num_mm_atoms
328  indmm = qmmm_env%mm_atom_index(imm)
329 
330  !add image forces to Forces
331  IF (qmmm_env%image_charge) THEN
332  DO iatom = 1, qmmm_env%num_image_mm_atoms
333  image_indmm = qmmm_env%image_charge_pot%image_mm_list(iatom)
334  IF (image_indmm .EQ. indmm) THEN
335  forces(:, imm) = forces(:, imm) &
336  + qmmm_env%image_charge_pot%image_forcesMM(:, iatom)
337  END IF
338  END DO
339  END IF
340 
341  ! Hack: In Forces there the gradients indeed...
342  ! Minux sign to take care of this misunderstanding...
343  mm_particles(indmm)%f(:) = -forces(:, imm) + mm_particles(indmm)%f(:)
344  END DO
345  DEALLOCATE (forces)
346  IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
347  DO imm = 1, qmmm_env%added_charges%num_mm_atoms
348  indmm = qmmm_env%added_charges%mm_atom_index(imm)
349  ! Hack: In Forces there the gradients indeed...
350  ! Minux sign to take care of this misunderstanding...
351  qmmm_env%added_charges%added_particles(indmm)%f(:) = -forces_added_charges(:, imm)
352  END DO
353  END IF
354  DEALLOCATE (forces_added_charges)
355  IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN
356  DO imm = 1, qmmm_env%added_shells%num_mm_atoms
357  indmm = qmmm_env%added_shells%mm_core_index(imm)
358  ! Hack: In Forces there the gradients indeed...
359  ! Minux sign to take care of this misunderstanding...
360  qmmm_env%added_shells%added_particles(imm)%f(:) = qmmm_env%added_shells%added_particles(imm)%f(:) - &
361  forces_added_shells(:, imm)
362 
363  END DO
364  END IF
365  DEALLOCATE (forces_added_shells)
366  END IF
367  CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
368  CALL timestop(handle)
369 
370  END SUBROUTINE qmmm_forces
371 
372 ! **************************************************************************************************
373 !> \brief Evaluates the contribution to the forces due to the
374 !> QM/MM potential computed collocating the Electrostatic
375 !> Gaussian Potential.
376 !> \param rho ...
377 !> \param qmmm_env ...
378 !> \param mm_particles ...
379 !> \param aug_pools ...
380 !> \param auxbas_grid ...
381 !> \param coarser_grid ...
382 !> \param cube_info ...
383 !> \param para_env ...
384 !> \param eps_mm_rspace ...
385 !> \param pw_pools ...
386 !> \param Forces ...
387 !> \param Forces_added_charges ...
388 !> \param Forces_added_shells ...
389 !> \param interp_section ...
390 !> \param iw ...
391 !> \param mm_cell ...
392 !> \par History
393 !> 06.2004 created [tlaino]
394 !> \author Teodoro Laino
395 ! **************************************************************************************************
396  SUBROUTINE qmmm_forces_with_gaussian(rho, qmmm_env, mm_particles, &
397  aug_pools, auxbas_grid, coarser_grid, cube_info, para_env, &
398  eps_mm_rspace, pw_pools, Forces, Forces_added_charges, Forces_added_shells, &
399  interp_section, iw, mm_cell)
400  TYPE(pw_r3d_rs_type), POINTER :: rho
401  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
402  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
403  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
404  INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
405  TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
406  TYPE(mp_para_env_type), POINTER :: para_env
407  REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
408  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
409  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges, &
410  forces_added_shells
411  TYPE(section_vals_type), POINTER :: interp_section
412  INTEGER, INTENT(IN) :: iw
413  TYPE(cell_type), POINTER :: mm_cell
414 
415  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian'
416 
417  INTEGER :: handle, i, igrid, j, k, kind_interp, me, &
418  ngrids, stat
419  INTEGER, DIMENSION(3) :: glb, gub, lb, ub
420  INTEGER, DIMENSION(:), POINTER :: pos_of_x
421  LOGICAL :: shells
422  REAL(kind=dp), DIMENSION(:, :), POINTER :: tmp
423  TYPE(mp_comm_type) :: group
424  TYPE(mp_request_type) :: request
425  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
426 
427 ! Statements
428 
429  CALL timeset(routinen, handle)
430  NULLIFY (tmp)
431  cpassert(ASSOCIATED(mm_particles))
432  cpassert(ASSOCIATED(qmmm_env%mm_atom_chrg))
433  cpassert(ASSOCIATED(qmmm_env%mm_atom_index))
434  cpassert(ASSOCIATED(forces))
435  !Statements
436  ngrids = SIZE(pw_pools)
437  CALL pw_pools_create_pws(aug_pools, grids)
438  DO igrid = 1, ngrids
439  CALL pw_zero(grids(igrid))
440  END DO
441  ! Collocate Density on multigrids
442  lb = rho%pw_grid%bounds_local(1, :)
443  ub = rho%pw_grid%bounds_local(2, :)
444  grids(auxbas_grid)%array(lb(1):ub(1), &
445  lb(2):ub(2), &
446  lb(3):ub(3)) = rho%array
447  ! copy the boundaries
448  DO i = lb(1), ub(1)
449  grids(auxbas_grid)%array(i, ub(2) + 1, ub(3) + 1) = rho%array(i, lb(2), lb(3))
450  END DO
451  DO k = lb(3), ub(3)
452  DO i = lb(1), ub(1)
453  grids(auxbas_grid)%array(i, ub(2) + 1, k) = rho%array(i, lb(2), k)
454  END DO
455  END DO
456  DO j = lb(2), ub(2)
457  DO i = lb(1), ub(1)
458  grids(auxbas_grid)%array(i, j, ub(3) + 1) = rho%array(i, j, lb(3))
459  END DO
460  END DO
461  pos_of_x => grids(auxbas_grid)%pw_grid%para%pos_of_x
462  group = grids(auxbas_grid)%pw_grid%para%group
463  me = grids(auxbas_grid)%pw_grid%para%my_pos
464  glb = rho%pw_grid%bounds(1, :)
465  gub = rho%pw_grid%bounds(2, :)
466  IF ((pos_of_x(glb(1)) .EQ. me) .AND. (pos_of_x(gub(1)) .EQ. me)) THEN
467  DO k = lb(3), ub(3)
468  DO j = lb(2), ub(2)
469  grids(auxbas_grid)%array(ub(1) + 1, j, k) = rho%array(lb(1), j, k)
470  END DO
471  grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, k) = rho%array(lb(1), lb(2), k)
472  END DO
473  DO j = lb(2), ub(2)
474  grids(auxbas_grid)%array(ub(1) + 1, j, ub(3) + 1) = rho%array(lb(1), j, lb(3))
475  END DO
476  grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, ub(3) + 1) = rho%array(lb(1), lb(2), lb(3))
477  ELSE IF (pos_of_x(glb(1)) .EQ. me) THEN
478  ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
479  rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)), &
480  stat=stat)
481  cpassert(stat == 0)
482  tmp = rho%array(lb(1), :, :)
483  CALL group%isend(msgin=tmp, dest=pos_of_x(rho%pw_grid%bounds(2, 1)), &
484  request=request, tag=112)
485  CALL request%wait()
486  ELSE IF (pos_of_x(gub(1)) .EQ. me) THEN
487  ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
488  rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)), &
489  stat=stat)
490  cpassert(stat == 0)
491  CALL group%irecv(msgout=tmp, source=pos_of_x(rho%pw_grid%bounds(1, 1)), &
492  request=request, tag=112)
493  CALL request%wait()
494 
495  DO k = lb(3), ub(3)
496  DO j = lb(2), ub(2)
497  grids(auxbas_grid)%array(ub(1) + 1, j, k) = tmp(j, k)
498  END DO
499  grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, k) = tmp(lb(2), k)
500  END DO
501  DO j = lb(2), ub(2)
502  grids(auxbas_grid)%array(ub(1) + 1, j, ub(3) + 1) = tmp(j, lb(3))
503  END DO
504  grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, ub(3) + 1) = tmp(lb(2), lb(3))
505  END IF
506  IF (ASSOCIATED(tmp)) THEN
507  DEALLOCATE (tmp)
508  END IF
509  ! Further setup of parallelization scheme
510  IF (qmmm_env%par_scheme == do_par_atom) THEN
511  CALL para_env%sum(grids(auxbas_grid)%array)
512  END IF
513  ! RealSpace Interpolation
514  CALL section_vals_val_get(interp_section, "kind", i_val=kind_interp)
515  SELECT CASE (kind_interp)
516  CASE (spline3_nopbc_interp, spline3_pbc_interp)
517  ! Spline Interpolator
518  DO igrid = auxbas_grid, SIZE(grids) - 1
519  CALL pw_restrict_s3(grids(igrid), &
520  grids(igrid + 1), &
521  aug_pools(igrid + 1)%pool, &
522  param_section=interp_section)
523  END DO
524  CASE DEFAULT
525  cpabort("")
526  END SELECT
527 
528  shells = .false.
529  CALL qmmm_force_with_gaussian_low(grids, mm_particles, &
530  qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
531  qmmm_env%num_mm_atoms, cube_info, para_env, eps_mm_rspace, auxbas_grid, &
532  coarser_grid, qmmm_env%pgfs, qmmm_env%potentials, forces, aug_pools, &
533  mm_cell, qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%per_potentials, &
534  iw, qmmm_env%par_scheme, qmmm_env%spherical_cutoff, shells)
535 
536  IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
537  CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_charges%added_particles, &
538  qmmm_env%added_charges%mm_atom_chrg, &
539  qmmm_env%added_charges%mm_atom_index, qmmm_env%added_charges%num_mm_atoms, &
540  cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_charges%pgfs, &
541  qmmm_env%added_charges%potentials, forces_added_charges, aug_pools, mm_cell, &
542  qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%added_charges%per_potentials, iw, qmmm_env%par_scheme, &
543  qmmm_env%spherical_cutoff, shells)
544  END IF
545 
546  IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN
547  shells = .true.
548  CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_shells%added_particles, &
549  qmmm_env%added_shells%mm_core_chrg, &
550  qmmm_env%added_shells%mm_core_index, qmmm_env%added_shells%num_mm_atoms, &
551  cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_shells%pgfs, &
552  qmmm_env%added_shells%potentials, forces_added_shells, aug_pools, mm_cell, &
553  qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%added_shells%per_potentials, iw, qmmm_env%par_scheme, &
554  qmmm_env%spherical_cutoff, shells)
555  END IF
556 
557  CALL pw_pools_give_back_pws(aug_pools, grids)
558  CALL timestop(handle)
559 
560  END SUBROUTINE qmmm_forces_with_gaussian
561 
562 ! **************************************************************************************************
563 !> \brief Evaluates the contribution to the forces due to the
564 !> QM/MM potential computed collocating the Electrostatic
565 !> Gaussian Potential. Low Level
566 !> \param grids ...
567 !> \param mm_particles ...
568 !> \param mm_charges ...
569 !> \param mm_atom_index ...
570 !> \param num_mm_atoms ...
571 !> \param cube_info ...
572 !> \param para_env ...
573 !> \param eps_mm_rspace ...
574 !> \param auxbas_grid ...
575 !> \param coarser_grid ...
576 !> \param pgfs ...
577 !> \param potentials ...
578 !> \param Forces ...
579 !> \param aug_pools ...
580 !> \param mm_cell ...
581 !> \param dOmmOqm ...
582 !> \param periodic ...
583 !> \param per_potentials ...
584 !> \param iw ...
585 !> \param par_scheme ...
586 !> \param qmmm_spherical_cutoff ...
587 !> \param shells ...
588 !> \par History
589 !> 06.2004 created [tlaino]
590 !> \author Teodoro Laino
591 ! **************************************************************************************************
592  SUBROUTINE qmmm_force_with_gaussian_low(grids, mm_particles, mm_charges, &
593  mm_atom_index, num_mm_atoms, cube_info, para_env, &
594  eps_mm_rspace, auxbas_grid, coarser_grid, pgfs, potentials, Forces, &
595  aug_pools, mm_cell, dOmmOqm, periodic, per_potentials, iw, par_scheme, &
596  qmmm_spherical_cutoff, shells)
597  TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: grids
598  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
599  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
600  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
601  INTEGER, INTENT(IN) :: num_mm_atoms
602  TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
603  TYPE(mp_para_env_type), POINTER :: para_env
604  REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
605  INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
606  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
607  TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
608  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
609  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
610  TYPE(cell_type), POINTER :: mm_cell
611  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
612  LOGICAL, INTENT(in) :: periodic
613  TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
614  INTEGER, INTENT(IN) :: iw, par_scheme
615  REAL(kind=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
616  LOGICAL, INTENT(in) :: shells
617 
618  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_force_with_gaussian_low', &
619  routinenb = 'qmmm_forces_gaussian_low'
620 
621  INTEGER :: handle, handle2, igauss, ilevel, imm, &
622  indmm, iradtyp, lindmm, myind, &
623  n_rep_real(3)
624  INTEGER, DIMENSION(2, 3) :: bo
625  REAL(kind=dp) :: alpha, dvol, height, sph_chrg_factor, w
626  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xdat, ydat, zdat
627  REAL(kind=dp), DIMENSION(3) :: force, ra
628  TYPE(qmmm_gaussian_type), POINTER :: pgf
629  TYPE(qmmm_per_pot_type), POINTER :: per_pot
630  TYPE(qmmm_pot_type), POINTER :: pot
631 
632  CALL timeset(routinen, handle)
633  CALL timeset(routinenb//"_G", handle2)
634  NULLIFY (pgf, pot, per_pot)
635  IF (par_scheme == do_par_atom) myind = 0
636  radius: DO iradtyp = 1, SIZE(pgfs)
637  pgf => pgfs(iradtyp)%pgf
638  pot => potentials(iradtyp)%pot
639  n_rep_real = 0
640  IF (periodic) THEN
641  per_pot => per_potentials(iradtyp)%pot
642  n_rep_real = per_pot%n_rep_real
643  END IF
644  gaussian: DO igauss = 1, pgf%Number_of_Gaussians
645  alpha = 1.0_dp/pgf%Gk(igauss)
646  alpha = alpha*alpha
647  height = pgf%Ak(igauss)
648  ilevel = pgf%grid_level(igauss)
649  dvol = grids(ilevel)%pw_grid%dvol
650  bo = grids(ilevel)%pw_grid%bounds_local
651  ALLOCATE (xdat(2, bo(1, 1):bo(2, 1)))
652  ALLOCATE (ydat(2, bo(1, 2):bo(2, 2)))
653  ALLOCATE (zdat(2, bo(1, 3):bo(2, 3)))
654  !$OMP PARALLEL DO DEFAULT(NONE) &
655  !$OMP SHARED(pot, par_scheme, dvol, alpha, para_env, mm_atom_index, shells) &
656  !$OMP SHARED(mm_particles, dOmmOqm, mm_cell, height, mm_charges, qmmm_spherical_cutoff) &
657  !$OMP SHARED(grids, cube_info, bo, n_rep_real, eps_mm_rspace, Forces, ilevel) &
658  !$OMP SHARED(IGauss, pgf, IRadTyp, iw, aug_pools, auxbas_grid) &
659  !$OMP PRIVATE(xdat, ydat, zdat) &
660  !$OMP PRIVATE(Imm, LIndMM, IndMM, ra, W, force, sph_chrg_factor, myind)
661  atoms: DO imm = 1, SIZE(pot%mm_atom_index)
662  IF (par_scheme == do_par_atom) THEN
663  myind = imm + (igauss - 1)*SIZE(pot%mm_atom_index) + (iradtyp - 1)*pgf%Number_of_Gaussians
664  IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
665  END IF
666  lindmm = pot%mm_atom_index(imm)
667  indmm = mm_atom_index(lindmm)
668  IF (shells) THEN
669  ra(:) = pbc(mm_particles(imm)%r - dommoqm, mm_cell) + dommoqm
670  ELSE
671  ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
672  END IF
673  w = mm_charges(lindmm)*height
674  force = 0.0_dp
675  ! Possible Spherical Cutoff
676  IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
677  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
678  w = w*sph_chrg_factor
679  END IF
680  IF (abs(w) <= epsilon(0.0_dp)) cycle atoms
681  CALL integrate_gf_rspace_nopbc(zetp=alpha, &
682  rp=ra, &
683  scale=-1.0_dp, &
684  w=w, &
685  pwgrid=grids(ilevel), &
686  cube_info=cube_info(ilevel), &
687  eps_mm_rspace=eps_mm_rspace, &
688  xdat=xdat, &
689  ydat=ydat, &
690  zdat=zdat, &
691  bo=bo, &
692  force=force, &
693  n_rep_real=n_rep_real, &
694  mm_cell=mm_cell)
695  force = force*dvol
696  forces(:, lindmm) = forces(:, lindmm) + force(:)
697  !
698  ! Debug Statement
699  !
700  IF (debug_this_module) THEN
701  CALL debug_integrate_gf_rspace_nopbc(ilevel=ilevel, &
702  zetp=alpha, &
703  rp=ra, &
704  w=w, &
705  pwgrid=grids(ilevel), &
706  cube_info=cube_info(ilevel), &
707  eps_mm_rspace=eps_mm_rspace, &
708  aug_pools=aug_pools, &
709  debug_force=force, &
710  mm_cell=mm_cell, &
711  auxbas_grid=auxbas_grid, &
712  n_rep_real=n_rep_real, &
713  iw=iw)
714  END IF
715  END DO atoms
716  !$OMP END PARALLEL DO
717  DEALLOCATE (xdat)
718  DEALLOCATE (ydat)
719  DEALLOCATE (zdat)
720  END DO gaussian
721  END DO radius
722  CALL timestop(handle2)
723  CALL timeset(routinenb//"_R", handle2)
724  IF (periodic) THEN
725  CALL qmmm_forces_with_gaussian_lg(pgfs=pgfs, &
726  cgrid=grids(coarser_grid), &
727  num_mm_atoms=num_mm_atoms, &
728  mm_charges=mm_charges, &
729  mm_atom_index=mm_atom_index, &
730  mm_particles=mm_particles, &
731  para_env=para_env, &
732  coarser_grid_level=coarser_grid, &
733  forces=forces, &
734  per_potentials=per_potentials, &
735  aug_pools=aug_pools, &
736  mm_cell=mm_cell, &
737  dommoqm=dommoqm, &
738  iw=iw, &
739  par_scheme=par_scheme, &
740  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
741  shells=shells)
742  ELSE
743  CALL qmmm_forces_with_gaussian_lr(pgfs=pgfs, &
744  cgrid=grids(coarser_grid), &
745  num_mm_atoms=num_mm_atoms, &
746  mm_charges=mm_charges, &
747  mm_atom_index=mm_atom_index, &
748  mm_particles=mm_particles, &
749  para_env=para_env, &
750  coarser_grid_level=coarser_grid, &
751  forces=forces, &
752  potentials=potentials, &
753  aug_pools=aug_pools, &
754  mm_cell=mm_cell, &
755  dommoqm=dommoqm, &
756  iw=iw, &
757  par_scheme=par_scheme, &
758  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
759  shells=shells)
760  END IF
761  CALL timestop(handle2)
762  CALL timestop(handle)
763  END SUBROUTINE qmmm_force_with_gaussian_low
764 
765 ! **************************************************************************************************
766 !> \brief Evaluates the contribution to the forces due to the Long Range
767 !> part of the QM/MM potential computed collocating the Electrostatic
768 !> Gaussian Potential.
769 !> \param pgfs ...
770 !> \param cgrid ...
771 !> \param num_mm_atoms ...
772 !> \param mm_charges ...
773 !> \param mm_atom_index ...
774 !> \param mm_particles ...
775 !> \param para_env ...
776 !> \param coarser_grid_level ...
777 !> \param Forces ...
778 !> \param per_potentials ...
779 !> \param aug_pools ...
780 !> \param mm_cell ...
781 !> \param dOmmOqm ...
782 !> \param iw ...
783 !> \param par_scheme ...
784 !> \param qmmm_spherical_cutoff ...
785 !> \param shells ...
786 !> \par History
787 !> 08.2004 created [tlaino]
788 !> \author Teodoro Laino
789 ! **************************************************************************************************
790  SUBROUTINE qmmm_forces_with_gaussian_lg(pgfs, cgrid, num_mm_atoms, mm_charges, mm_atom_index, &
791  mm_particles, para_env, coarser_grid_level, Forces, per_potentials, &
792  aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
793  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
794  TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
795  INTEGER, INTENT(IN) :: num_mm_atoms
796  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
797  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
798  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
799  TYPE(mp_para_env_type), POINTER :: para_env
800  INTEGER, INTENT(IN) :: coarser_grid_level
801  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
802  TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
803  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
804  TYPE(cell_type), POINTER :: mm_cell
805  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
806  INTEGER, INTENT(IN) :: iw, par_scheme
807  REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
808  LOGICAL :: shells
809 
810  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian_LG'
811 
812  INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, ij3, ij4, ik1, ik2, ik3, ik4, imm, &
813  indmm, iradtyp, ivec(3), j, k, lindmm, my_i, my_j, my_k, myind, npts(3)
814  INTEGER, DIMENSION(2, 3) :: bo, gbo
815  REAL(kind=dp) :: a1, a2, a3, abc_x(4, 4), abc_x_y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, &
816  dr1, dr1c, dr1i, dr2, dr2c, dr2i, dr3, dr3c, dr3i, dvol, e1, e2, e3, f1, f2, f3, fac, &
817  ft1, ft2, ft3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, &
818  rt3, rv1, rv2, rv3, s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, &
819  sph_chrg_factor, t1, t1d, t1o, t2, t2d, t2o, t3, t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, &
820  v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, v4o, xd1, xd2, xd3, xs1, xs2, xs3
821  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: lforces
822  REAL(kind=dp), DIMENSION(3) :: ra, val, vec
823  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid, grid2
824  TYPE(pw_r3d_rs_type), POINTER :: pw
825  TYPE(qmmm_per_pot_type), POINTER :: per_pot
826 
827  CALL timeset(routinen, handle)
828  NULLIFY (grid)
829  ALLOCATE (lforces(3, num_mm_atoms))
830  lforces = 0.0_dp
831  dr1c = cgrid%pw_grid%dr(1)
832  dr2c = cgrid%pw_grid%dr(2)
833  dr3c = cgrid%pw_grid%dr(3)
834  dvol = cgrid%pw_grid%dvol
835  gbo = cgrid%pw_grid%bounds
836  bo = cgrid%pw_grid%bounds_local
837  grid => cgrid%array
838  IF (par_scheme == do_par_atom) myind = 0
839  radius: DO iradtyp = 1, SIZE(pgfs)
840  per_pot => per_potentials(iradtyp)%pot
841  pw => per_pot%TabLR
842  grid2 => pw%array(:, :, :)
843  npts = pw%pw_grid%npts
844  dr1 = pw%pw_grid%dr(1)
845  dr2 = pw%pw_grid%dr(2)
846  dr3 = pw%pw_grid%dr(3)
847  dr1i = 1.0_dp/dr1
848  dr2i = 1.0_dp/dr2
849  dr3i = 1.0_dp/dr3
850 
851  !$OMP PARALLEL DO DEFAULT(NONE) &
852  !$OMP SHARED(bo, grid, grid2, pw, npts, gbo, per_pot, mm_atom_index) &
853  !$OMP SHARED(dr1, dr2, dr3, dr1i, dr2i, dr3i, dr1c, dr2c, dr3c, par_scheme, mm_charges) &
854  !$OMP SHARED(mm_cell, dOmmOqm, dvol, shells, para_env, IRadTyp) &
855  !$OMP SHARED(qmmm_spherical_cutoff, mm_particles, Forces, LForces) &
856  !$OMP PRIVATE(qt, Imm, LIndMM, IndMM, sph_chrg_factor, ra, myind) &
857  !$OMP PRIVATE(rt1, rt2, rt3, ft1, ft2, ft3, my_k, my_j, my_i, xs3, xs2, xs1) &
858  !$OMP PRIVATE(rv3, rv2, rv1, vec, ivec, ik1, ik2, ik3, ik4, xd3, xd2, xd1) &
859  !$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, u1, u2, u3, v1o, v2o, v3o, v4o) &
860  !$OMP PRIVATE(v1d, v2d, v3d, v4d, ij1, ij2, ij3, ij4, e1, e2, e3, f1, f2, f3) &
861  !$OMP PRIVATE(g1, g2, g3, h1, h2, h3, s1o, s2o, s3o, s4o, s1d, s2d, s3d, s4d) &
862  !$OMP PRIVATE(ii1, ii2, ii3, ii4, a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3) &
863  !$OMP PRIVATE(t1o, t2o, t3o, t4o, t1d, t2d, t3d, t4d, t1, t2, t3, t4, s1, s2, s3, s4) &
864  !$OMP PRIVATE(v1, v2, v3, v4, abc_x, abc_x_y, val, fac)
865  atoms: DO imm = 1, SIZE(per_pot%mm_atom_index)
866  IF (par_scheme == do_par_atom) THEN
867  myind = imm + (iradtyp - 1)*SIZE(per_pot%mm_atom_index)
868  IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle
869  END IF
870  lindmm = per_pot%mm_atom_index(imm)
871  indmm = mm_atom_index(lindmm)
872  IF (shells) THEN
873  ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
874  ELSE
875  ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
876  END IF
877  qt = mm_charges(lindmm)
878  ! Possible Spherical Cutoff
879  IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
880  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
881  qt = qt*sph_chrg_factor
882  END IF
883  IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
884  rt1 = ra(1)
885  rt2 = ra(2)
886  rt3 = ra(3)
887  ft1 = 0.0_dp
888  ft2 = 0.0_dp
889  ft3 = 0.0_dp
890  loopongrid: DO k = bo(1, 3), bo(2, 3)
891  my_k = k - gbo(1, 3)
892  xs3 = real(my_k, dp)*dr3c
893  my_j = bo(1, 2) - gbo(1, 2)
894  xs2 = real(my_j, dp)*dr2c
895  rv3 = rt3 - xs3
896  vec(3) = rv3
897  ivec(3) = floor(vec(3)/pw%pw_grid%dr(3))
898  ik1 = modulo(ivec(3) - 1, npts(3)) + 1
899  ik2 = modulo(ivec(3), npts(3)) + 1
900  ik3 = modulo(ivec(3) + 1, npts(3)) + 1
901  ik4 = modulo(ivec(3) + 2, npts(3)) + 1
902  xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
903  p1 = 3.0_dp + xd3
904  p2 = p1*p1
905  p3 = p2*p1
906  q1 = 2.0_dp + xd3
907  q2 = q1*q1
908  q3 = q2*q1
909  r1 = 1.0_dp + xd3
910  r2 = r1*r1
911  r3 = r2*r1
912  u1 = xd3
913  u2 = u1*u1
914  u3 = u2*u1
915  v1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
916  v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
917  v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
918  v4o = 1.0_dp/6.0_dp*u3
919  v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2
920  v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2
921  v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2
922  v4d = 0.5_dp*u2
923  DO j = bo(1, 2), bo(2, 2)
924  my_i = bo(1, 1) - gbo(1, 1)
925  xs1 = real(my_i, dp)*dr1c
926  rv2 = rt2 - xs2
927  vec(2) = rv2
928  ivec(2) = floor(vec(2)/pw%pw_grid%dr(2))
929  ij1 = modulo(ivec(2) - 1, npts(2)) + 1
930  ij2 = modulo(ivec(2), npts(2)) + 1
931  ij3 = modulo(ivec(2) + 1, npts(2)) + 1
932  ij4 = modulo(ivec(2) + 2, npts(2)) + 1
933  xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
934  e1 = 3.0_dp + xd2
935  e2 = e1*e1
936  e3 = e2*e1
937  f1 = 2.0_dp + xd2
938  f2 = f1*f1
939  f3 = f2*f1
940  g1 = 1.0_dp + xd2
941  g2 = g1*g1
942  g3 = g2*g1
943  h1 = xd2
944  h2 = h1*h1
945  h3 = h2*h1
946  s1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
947  s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
948  s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
949  s4o = 1.0_dp/6.0_dp*h3
950  s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2
951  s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2
952  s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2
953  s4d = 0.5_dp*h2
954  DO i = bo(1, 1), bo(2, 1)
955  rv1 = rt1 - xs1
956  vec(1) = rv1
957  ivec(1) = floor(vec(1)/pw%pw_grid%dr(1))
958  ii1 = modulo(ivec(1) - 1, npts(1)) + 1
959  ii2 = modulo(ivec(1), npts(1)) + 1
960  ii3 = modulo(ivec(1) + 1, npts(1)) + 1
961  ii4 = modulo(ivec(1) + 2, npts(1)) + 1
962  xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
963  a1 = 3.0_dp + xd1
964  a2 = a1*a1
965  a3 = a2*a1
966  b1 = 2.0_dp + xd1
967  b2 = b1*b1
968  b3 = b2*b1
969  c1 = 1.0_dp + xd1
970  c2 = c1*c1
971  c3 = c2*c1
972  d1 = xd1
973  d2 = d1*d1
974  d3 = d2*d1
975  t1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
976  t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
977  t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
978  t4o = 1.0_dp/6.0_dp*d3
979  t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2
980  t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2
981  t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2
982  t4d = 0.5_dp*d2
983 
984  t1 = t1d*dr1i
985  t2 = t2d*dr1i
986  t3 = t3d*dr1i
987  t4 = t4d*dr1i
988  s1 = s1o
989  s2 = s2o
990  s3 = s3o
991  s4 = s4o
992  v1 = v1o
993  v2 = v2o
994  v3 = v3o
995  v4 = v4o
996 
997  abc_x(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
998  abc_x(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
999  abc_x(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
1000  abc_x(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
1001  abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1002 
1003  abc_x(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
1004  abc_x(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
1005  abc_x(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
1006  abc_x(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
1007  abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1008 
1009  abc_x(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
1010  abc_x(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
1011  abc_x(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
1012  abc_x(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
1013  abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1014 
1015  abc_x(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
1016  abc_x(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
1017  abc_x(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
1018  abc_x(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
1019  abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1020 
1021  val(1) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1022 
1023  t1 = t1o
1024  t2 = t2o
1025  t3 = t3o
1026  t4 = t4o
1027  s1 = s1d*dr2i
1028  s2 = s2d*dr2i
1029  s3 = s3d*dr2i
1030  s4 = s4d*dr2i
1031 
1032  abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1033  abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1034  abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1035  abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1036 
1037  val(2) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1038 
1039  t1 = t1o
1040  t2 = t2o
1041  t3 = t3o
1042  t4 = t4o
1043  s1 = s1o
1044  s2 = s2o
1045  s3 = s3o
1046  s4 = s4o
1047  v1 = v1d*dr3i
1048  v2 = v2d*dr3i
1049  v3 = v3d*dr3i
1050  v4 = v4d*dr3i
1051 
1052  abc_x(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
1053  abc_x(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
1054  abc_x(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
1055  abc_x(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
1056  abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1057  abc_x(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
1058  abc_x(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
1059  abc_x(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
1060  abc_x(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
1061  abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1062  abc_x(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
1063  abc_x(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
1064  abc_x(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
1065  abc_x(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
1066  abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1067  abc_x(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
1068  abc_x(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
1069  abc_x(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
1070  abc_x(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
1071  abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1072 
1073  val(3) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1074 
1075  fac = grid(i, j, k)
1076  ft1 = ft1 + val(1)*fac
1077  ft2 = ft2 + val(2)*fac
1078  ft3 = ft3 + val(3)*fac
1079  xs1 = xs1 + dr1c
1080  END DO
1081  xs2 = xs2 + dr2c
1082  END DO
1083  END DO loopongrid
1084  qt = -qt*dvol
1085  lforces(1, lindmm) = ft1*qt
1086  lforces(2, lindmm) = ft2*qt
1087  lforces(3, lindmm) = ft3*qt
1088 
1089  forces(1, lindmm) = forces(1, lindmm) + lforces(1, lindmm)
1090  forces(2, lindmm) = forces(2, lindmm) + lforces(2, lindmm)
1091  forces(3, lindmm) = forces(3, lindmm) + lforces(3, lindmm)
1092  END DO atoms
1093  !$OMP END PARALLEL DO
1094  END DO radius
1095  !
1096  ! Debug Statement
1097  !
1098  IF (debug_this_module) THEN
1099  CALL debug_qmmm_forces_with_gauss_lg(pgfs=pgfs, &
1100  aug_pools=aug_pools, &
1101  rho=cgrid, &
1102  num_mm_atoms=num_mm_atoms, &
1103  mm_charges=mm_charges, &
1104  mm_atom_index=mm_atom_index, &
1105  mm_particles=mm_particles, &
1106  coarser_grid_level=coarser_grid_level, &
1107  debug_force=lforces, &
1108  per_potentials=per_potentials, &
1109  para_env=para_env, &
1110  mm_cell=mm_cell, &
1111  dommoqm=dommoqm, &
1112  iw=iw, &
1113  par_scheme=par_scheme, &
1114  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1115  shells=shells)
1116  END IF
1117  DEALLOCATE (lforces)
1118  CALL timestop(handle)
1119  END SUBROUTINE qmmm_forces_with_gaussian_lg
1120 
1121 ! **************************************************************************************************
1122 !> \brief Evaluates the contribution to the forces due to the Long Range
1123 !> part of the QM/MM potential computed collocating the Electrostatic
1124 !> Gaussian Potential.
1125 !> \param pgfs ...
1126 !> \param cgrid ...
1127 !> \param num_mm_atoms ...
1128 !> \param mm_charges ...
1129 !> \param mm_atom_index ...
1130 !> \param mm_particles ...
1131 !> \param para_env ...
1132 !> \param coarser_grid_level ...
1133 !> \param Forces ...
1134 !> \param potentials ...
1135 !> \param aug_pools ...
1136 !> \param mm_cell ...
1137 !> \param dOmmOqm ...
1138 !> \param iw ...
1139 !> \param par_scheme ...
1140 !> \param qmmm_spherical_cutoff ...
1141 !> \param shells ...
1142 !> \par History
1143 !> 08.2004 created [tlaino]
1144 !> \author Teodoro Laino
1145 ! **************************************************************************************************
1146  SUBROUTINE qmmm_forces_with_gaussian_lr(pgfs, cgrid, num_mm_atoms, mm_charges, mm_atom_index, &
1147  mm_particles, para_env, coarser_grid_level, Forces, potentials, &
1148  aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1149  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1150  TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
1151  INTEGER, INTENT(IN) :: num_mm_atoms
1152  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1153  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1154  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1155  TYPE(mp_para_env_type), POINTER :: para_env
1156  INTEGER, INTENT(IN) :: coarser_grid_level
1157  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
1158  TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1159  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1160  TYPE(cell_type), POINTER :: mm_cell
1161  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1162  INTEGER, INTENT(IN) :: iw, par_scheme
1163  REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1164  LOGICAL :: shells
1165 
1166  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian_LR'
1167 
1168  INTEGER :: handle, i, imm, indmm, iradtyp, ix, j, &
1169  k, lindmm, my_i, my_j, my_k, myind, &
1170  n1, n2, n3
1171  INTEGER, DIMENSION(2, 3) :: bo, gbo
1172  REAL(kind=dp) :: dr1, dr2, dr3, dvol, dx, fac, ft1, ft2, &
1173  ft3, qt, r, r2, rd1, rd2, rd3, rt1, &
1174  rt2, rt3, rv1, rv2, rv3, rx, rx2, &
1175  sph_chrg_factor, term, xs1, xs2, xs3
1176  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: lforces
1177  REAL(kind=dp), DIMENSION(3) :: ra
1178  REAL(kind=dp), DIMENSION(:, :), POINTER :: pot0_2
1179  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
1180  TYPE(qmmm_pot_type), POINTER :: pot
1181 
1182  CALL timeset(routinen, handle)
1183  ALLOCATE (lforces(3, num_mm_atoms))
1184  lforces = 0.0_dp
1185  n1 = cgrid%pw_grid%npts(1)
1186  n2 = cgrid%pw_grid%npts(2)
1187  n3 = cgrid%pw_grid%npts(3)
1188  dr1 = cgrid%pw_grid%dr(1)
1189  dr2 = cgrid%pw_grid%dr(2)
1190  dr3 = cgrid%pw_grid%dr(3)
1191  dvol = cgrid%pw_grid%dvol
1192  gbo = cgrid%pw_grid%bounds
1193  bo = cgrid%pw_grid%bounds_local
1194  grid => cgrid%array
1195  IF (par_scheme == do_par_atom) myind = 0
1196  radius: DO iradtyp = 1, SIZE(pgfs)
1197  pot => potentials(iradtyp)%pot
1198  dx = pot%dx
1199  pot0_2 => pot%pot0_2
1200  !$OMP PARALLEL DO DEFAULT(NONE) &
1201  !$OMP SHARED(pot, par_scheme, para_env, dvol, mm_atom_index, mm_particles, dOmmOqm) &
1202  !$OMP SHARED(mm_cell, mm_charges, dx, LForces, Forces, qmmm_spherical_cutoff, shells, dr1, dr2, dr3, gbo, bo) &
1203  !$OMP SHARED(IRadTyp, pot0_2, grid) &
1204  !$OMP PRIVATE(Imm, myind, ra, LIndMM, IndMM, qt, rt1, rt2, rt3, ft1, ft2, ft3, i, j, k, sph_chrg_factor) &
1205  !$OMP PRIVATE(my_k, my_j, my_i, xs3, xs2, xs1, rv1, rv2, rv3, r, ix, rx, rx2, r2, Term, fac) &
1206  !$OMP PRIVATE(rd1, rd2, rd3)
1207  atoms: DO imm = 1, SIZE(pot%mm_atom_index)
1208  IF (par_scheme == do_par_atom) THEN
1209  myind = imm + (iradtyp - 1)*SIZE(pot%mm_atom_index)
1210  IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle
1211  END IF
1212  lindmm = pot%mm_atom_index(imm)
1213  indmm = mm_atom_index(lindmm)
1214  ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
1215  IF (shells) &
1216  ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
1217  qt = mm_charges(lindmm)
1218  ! Possible Spherical Cutoff
1219  IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1220  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
1221  qt = qt*sph_chrg_factor
1222  END IF
1223  IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
1224  rt1 = ra(1)
1225  rt2 = ra(2)
1226  rt3 = ra(3)
1227  ft1 = 0.0_dp
1228  ft2 = 0.0_dp
1229  ft3 = 0.0_dp
1230  loopongrid: DO k = bo(1, 3), bo(2, 3)
1231  my_k = k - gbo(1, 3)
1232  xs3 = real(my_k, dp)*dr3
1233  my_j = bo(1, 2) - gbo(1, 2)
1234  xs2 = real(my_j, dp)*dr2
1235  rv3 = rt3 - xs3
1236  DO j = bo(1, 2), bo(2, 2)
1237  my_i = bo(1, 1) - gbo(1, 1)
1238  xs1 = real(my_i, dp)*dr1
1239  rv2 = rt2 - xs2
1240  DO i = bo(1, 1), bo(2, 1)
1241  rv1 = rt1 - xs1
1242  r2 = rv1*rv1 + rv2*rv2 + rv3*rv3
1243  r = sqrt(r2)
1244  ix = floor(r/dx) + 1
1245  rx = (r - real(ix - 1, dp)*dx)/dx
1246  rx2 = rx*rx
1247  term = pot0_2(1, ix)*(-6.0_dp*(rx - rx2)) &
1248  + pot0_2(2, ix)*(1.0_dp - 4.0_dp*rx + 3.0_dp*rx2) &
1249  + pot0_2(1, ix + 1)*(6.0_dp*(rx - rx2)) &
1250  + pot0_2(2, ix + 1)*(-2.0_dp*rx + 3.0_dp*rx2)
1251  fac = grid(i, j, k)*term
1252  IF (r == 0.0_dp) THEN
1253  rd1 = 1.0_dp
1254  rd2 = 1.0_dp
1255  rd3 = 1.0_dp
1256  ELSE
1257  rd1 = rv1/r
1258  rd2 = rv2/r
1259  rd3 = rv3/r
1260  END IF
1261  ft1 = ft1 + fac*rd1
1262  ft2 = ft2 + fac*rd2
1263  ft3 = ft3 + fac*rd3
1264  xs1 = xs1 + dr1
1265  END DO
1266  xs2 = xs2 + dr2
1267  END DO
1268  END DO loopongrid
1269  qt = -qt*dvol/dx
1270  lforces(1, lindmm) = ft1*qt
1271  lforces(2, lindmm) = ft2*qt
1272  lforces(3, lindmm) = ft3*qt
1273 
1274  forces(1, lindmm) = forces(1, lindmm) + lforces(1, lindmm)
1275  forces(2, lindmm) = forces(2, lindmm) + lforces(2, lindmm)
1276  forces(3, lindmm) = forces(3, lindmm) + lforces(3, lindmm)
1277  END DO atoms
1278  !$OMP END PARALLEL DO
1279  END DO radius
1280  !
1281  ! Debug Statement
1282  !
1283  IF (debug_this_module) THEN
1284  CALL debug_qmmm_forces_with_gauss_lr(pgfs=pgfs, &
1285  aug_pools=aug_pools, &
1286  rho=cgrid, &
1287  num_mm_atoms=num_mm_atoms, &
1288  mm_charges=mm_charges, &
1289  mm_atom_index=mm_atom_index, &
1290  mm_particles=mm_particles, &
1291  coarser_grid_level=coarser_grid_level, &
1292  debug_force=lforces, &
1293  potentials=potentials, &
1294  para_env=para_env, &
1295  mm_cell=mm_cell, &
1296  dommoqm=dommoqm, &
1297  iw=iw, &
1298  par_scheme=par_scheme, &
1299  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1300  shells=shells)
1301  END IF
1302 
1303  DEALLOCATE (lforces)
1304  CALL timestop(handle)
1305  END SUBROUTINE qmmm_forces_with_gaussian_lr
1306 
1307 ! **************************************************************************************************
1308 !> \brief Evaluates numerically QM/MM forces and compares them with
1309 !> the analytically computed ones.
1310 !> It is evaluated only when debug_this_module is set to .TRUE.
1311 !> \param rho ...
1312 !> \param qs_env ...
1313 !> \param qmmm_env ...
1314 !> \param Analytical_Forces ...
1315 !> \param mm_particles ...
1316 !> \param mm_atom_index ...
1317 !> \param num_mm_atoms ...
1318 !> \param interp_section ...
1319 !> \param mm_cell ...
1320 !> \par History
1321 !> 08.2004 created [tlaino]
1322 !> \author Teodoro Laino
1323 ! **************************************************************************************************
1324  SUBROUTINE qmmm_debug_forces(rho, qs_env, qmmm_env, Analytical_Forces, &
1325  mm_particles, mm_atom_index, num_mm_atoms, &
1326  interp_section, mm_cell)
1327  TYPE(pw_r3d_rs_type), POINTER :: rho
1328  TYPE(qs_environment_type), POINTER :: qs_env
1329  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1330  REAL(kind=dp), DIMENSION(:, :), POINTER :: analytical_forces
1331  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1332  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1333  INTEGER, INTENT(IN) :: num_mm_atoms
1334  TYPE(section_vals_type), POINTER :: interp_section
1335  TYPE(cell_type), POINTER :: mm_cell
1336 
1337  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_debug_forces'
1338 
1339  INTEGER :: handle, i, indmm, iw, j, k
1340  REAL(kind=dp) :: coord_save
1341  REAL(kind=dp), DIMENSION(2) :: energy
1342  REAL(kind=dp), DIMENSION(3) :: err
1343  REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1344  TYPE(cp_logger_type), POINTER :: logger
1345  TYPE(mp_para_env_type), POINTER :: para_env
1346  TYPE(pw_env_type), POINTER :: pw_env
1347  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1348  TYPE(pw_r3d_rs_type) :: v_qmmm_rspace
1349  TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
1350  TYPE(section_vals_type), POINTER :: input_section, print_section
1351 
1352  CALL timeset(routinen, handle)
1353  NULLIFY (num_forces)
1354  CALL get_qs_env(qs_env=qs_env, &
1355  pw_env=pw_env, &
1356  input=input_section, &
1357  para_env=para_env)
1358 
1359  print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
1360  logger => cp_get_default_logger()
1361  iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".qmmmLog")
1362  CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools)
1363  CALL pw_pools(1)%pool%create_pw(v_qmmm_rspace)
1364  ALLOCATE (num_forces(3, num_mm_atoms))
1365  ks_qmmm_env_loc => qs_env%ks_qmmm_env
1366  IF (iw > 0) WRITE (iw, '(/A)') "DEBUG SECTION:"
1367  atoms: DO i = 1, num_mm_atoms
1368  indmm = mm_atom_index(i)
1369  coords: DO j = 1, 3
1370  coord_save = mm_particles(indmm)%r(j)
1371  energy = 0.0_dp
1372  diff: DO k = 1, 2
1373  mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1374  CALL pw_zero(v_qmmm_rspace)
1375  SELECT CASE (qmmm_env%qmmm_coupl_type)
1376  CASE (do_qmmm_coulomb)
1377  cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1378  CASE (do_qmmm_pcharge)
1379  cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1380  CASE (do_qmmm_gauss, do_qmmm_swave)
1381  CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env, &
1382  v_qmmm=v_qmmm_rspace, &
1383  mm_particles=mm_particles, &
1384  aug_pools=qmmm_env%aug_pools, &
1385  para_env=para_env, &
1386  eps_mm_rspace=qmmm_env%eps_mm_rspace, &
1387  cube_info=ks_qmmm_env_loc%cube_info, &
1388  pw_pools=pw_pools, &
1389  auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
1390  coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
1391  interp_section=interp_section, &
1392  mm_cell=mm_cell)
1393  CASE (do_qmmm_none)
1394  cycle diff
1395  CASE DEFAULT
1396  IF (iw > 0) WRITE (iw, '(T3,A)') "Unknown Coupling..."
1397  cpabort("")
1398  END SELECT
1399  energy(k) = pw_integral_ab(rho, v_qmmm_rspace)
1400  END DO diff
1401  IF (iw > 0) THEN
1402  WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1403  "DEBUG :: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1404  END IF
1405  num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1406  mm_particles(indmm)%r(j) = coord_save
1407  END DO coords
1408  END DO atoms
1409 
1410  SELECT CASE (qmmm_env%qmmm_coupl_type)
1411  CASE (do_qmmm_coulomb)
1412  cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1413  CASE (do_qmmm_pcharge)
1414  cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1415  CASE (do_qmmm_gauss, do_qmmm_swave)
1416  IF (iw > 0) WRITE (iw, '(/A/)') "CHECKING NUMERICAL Vs ANALYTICAL FORCES (Err%):"
1417  DO i = 1, num_mm_atoms
1418  indmm = mm_atom_index(i)
1419  err = 0.0_dp
1420  DO k = 1, 3
1421  IF (abs(num_forces(k, i)) >= 5.0e-5_dp) THEN
1422  err(k) = (analytical_forces(k, i) - num_forces(k, i))/num_forces(k, i)*100.0_dp
1423  END IF
1424  END DO
1425  IF (iw > 0) &
1426  WRITE (iw, 100) indmm, analytical_forces(1, i), num_forces(1, i), err(1), &
1427  analytical_forces(2, i), num_forces(2, i), err(2), &
1428  analytical_forces(3, i), num_forces(3, i), err(3)
1429  cpassert(abs(err(1)) <= maxerr)
1430  cpassert(abs(err(2)) <= maxerr)
1431  cpassert(abs(err(3)) <= maxerr)
1432  END DO
1433  CASE (do_qmmm_none)
1434  IF (iw > 0) WRITE (iw, '(T3,A)') "No QM/MM Derivatives to debug. Just Mechanical Coupling!"
1435  CASE DEFAULT
1436  IF (iw > 0) WRITE (iw, '(T3,A)') "Unknown Coupling..."
1437  cpabort("")
1438  END SELECT
1439  CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
1440 
1441  CALL pw_pools(1)%pool%give_back_pw(v_qmmm_rspace)
1442  DEALLOCATE (num_forces)
1443  CALL timestop(handle)
1444 100 FORMAT(i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1445  END SUBROUTINE qmmm_debug_forces
1446 
1447 ! **************************************************************************************************
1448 !> \brief Debugs the integrate_gf_rspace_NoPBC.. It may helps ;-P
1449 !> \param ilevel ...
1450 !> \param zetp ...
1451 !> \param rp ...
1452 !> \param W ...
1453 !> \param pwgrid ...
1454 !> \param cube_info ...
1455 !> \param eps_mm_rspace ...
1456 !> \param aug_pools ...
1457 !> \param debug_force ...
1458 !> \param mm_cell ...
1459 !> \param auxbas_grid ...
1460 !> \param n_rep_real ...
1461 !> \param iw ...
1462 !> \par History
1463 !> 08.2004 created [tlaino]
1464 !> \author Teodoro Laino
1465 ! **************************************************************************************************
1466  SUBROUTINE debug_integrate_gf_rspace_nopbc(ilevel, zetp, rp, W, pwgrid, cube_info, &
1467  eps_mm_rspace, aug_pools, debug_force, &
1468  mm_cell, auxbas_grid, n_rep_real, iw)
1469  INTEGER, INTENT(IN) :: ilevel
1470  REAL(kind=dp), INTENT(IN) :: zetp
1471  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rp
1472  REAL(kind=dp), INTENT(IN) :: w
1473  TYPE(pw_r3d_rs_type), INTENT(IN) :: pwgrid
1474  TYPE(cube_info_type), INTENT(IN) :: cube_info
1475  REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
1476  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1477  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: debug_force
1478  TYPE(cell_type), POINTER :: mm_cell
1479  INTEGER, INTENT(IN) :: auxbas_grid
1480  INTEGER, DIMENSION(3), INTENT(IN) :: n_rep_real
1481  INTEGER, INTENT(IN) :: iw
1482 
1483  CHARACTER(len=*), PARAMETER :: routinen = 'debug_integrate_gf_rspace_NoPBC'
1484 
1485  INTEGER :: handle, i, igrid, k, ngrids
1486  INTEGER, DIMENSION(2, 3) :: bo2
1487  INTEGER, SAVE :: icount
1488  REAL(kind=dp), DIMENSION(2) :: energy
1489  REAL(kind=dp), DIMENSION(3) :: err, force, myrp
1490  REAL(kind=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat
1491  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1492 
1493  DATA icount/0/
1494  ! Statements
1495  CALL timeset(routinen, handle)
1496  !Statements
1497  ngrids = SIZE(aug_pools)
1498  CALL pw_pools_create_pws(aug_pools, grids)
1499  DO igrid = 1, ngrids
1500  CALL pw_zero(grids(igrid))
1501  END DO
1502  bo2 = grids(auxbas_grid)%pw_grid%bounds
1503  ALLOCATE (xdat(bo2(1, 1):bo2(2, 1)))
1504  ALLOCATE (ydat(bo2(1, 2):bo2(2, 2)))
1505  ALLOCATE (zdat(bo2(1, 3):bo2(2, 3)))
1506 
1507  icount = icount + 1
1508  DO i = 1, 3
1509  DO k = 1, 2
1510  myrp = rp
1511  myrp(i) = myrp(i) + (-1.0_dp)**k*dx
1512  CALL pw_zero(grids(ilevel))
1513  CALL collocate_gf_rspace_nopbc(zetp=zetp, &
1514  rp=myrp, &
1515  scale=-1.0_dp, &
1516  w=w, &
1517  pwgrid=grids(ilevel), &
1518  cube_info=cube_info, &
1519  eps_mm_rspace=eps_mm_rspace, &
1520  xdat=xdat, &
1521  ydat=ydat, &
1522  zdat=zdat, &
1523  bo2=bo2, &
1524  n_rep_real=n_rep_real, &
1525  mm_cell=mm_cell)
1526 
1527  energy(k) = pw_integral_ab(pwgrid, grids(ilevel))
1528  END DO
1529  force(i) = (energy(2) - energy(1))/(2.0_dp*dx)
1530  END DO
1531  err = 0.0_dp
1532  IF (all(force /= 0.0_dp)) THEN
1533  err(1) = (debug_force(1) - force(1))/force(1)*100.0_dp
1534  err(2) = (debug_force(2) - force(2))/force(2)*100.0_dp
1535  err(3) = (debug_force(3) - force(3))/force(3)*100.0_dp
1536  END IF
1537  IF (iw > 0) &
1538  WRITE (iw, 100) icount, debug_force(1), force(1), err(1), &
1539  debug_force(2), force(2), err(2), &
1540  debug_force(3), force(3), err(3)
1541  cpassert(abs(err(1)) <= maxerr)
1542  cpassert(abs(err(2)) <= maxerr)
1543  cpassert(abs(err(3)) <= maxerr)
1544 
1545  IF (ASSOCIATED(xdat)) THEN
1546  DEALLOCATE (xdat)
1547  END IF
1548  IF (ASSOCIATED(ydat)) THEN
1549  DEALLOCATE (ydat)
1550  END IF
1551  IF (ASSOCIATED(zdat)) THEN
1552  DEALLOCATE (zdat)
1553  END IF
1554 
1555  CALL pw_pools_give_back_pws(aug_pools, grids)
1556  CALL timestop(handle)
1557 100 FORMAT("Collocation : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1558  END SUBROUTINE debug_integrate_gf_rspace_nopbc
1559 
1560 ! **************************************************************************************************
1561 !> \brief Debugs qmmm_forces_with_gaussian_LG.. It may helps too ... ;-]
1562 !> \param pgfs ...
1563 !> \param aug_pools ...
1564 !> \param rho ...
1565 !> \param mm_charges ...
1566 !> \param mm_atom_index ...
1567 !> \param mm_particles ...
1568 !> \param num_mm_atoms ...
1569 !> \param coarser_grid_level ...
1570 !> \param per_potentials ...
1571 !> \param debug_force ...
1572 !> \param para_env ...
1573 !> \param mm_cell ...
1574 !> \param dOmmOqm ...
1575 !> \param iw ...
1576 !> \param par_scheme ...
1577 !> \param qmmm_spherical_cutoff ...
1578 !> \param shells ...
1579 !> \par History
1580 !> 08.2004 created [tlaino]
1581 !> \author Teodoro Laino
1582 ! **************************************************************************************************
1583  SUBROUTINE debug_qmmm_forces_with_gauss_lg(pgfs, aug_pools, rho, mm_charges, mm_atom_index, &
1584  mm_particles, num_mm_atoms, coarser_grid_level, per_potentials, &
1585  debug_force, para_env, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1586 
1587  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1588  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1589  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1590  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1591  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1592  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1593  INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level
1594  TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
1595  REAL(kind=dp), DIMENSION(:, :) :: debug_force
1596  TYPE(mp_para_env_type), POINTER :: para_env
1597  TYPE(cell_type), POINTER :: mm_cell
1598  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1599  INTEGER, INTENT(IN) :: iw, par_scheme
1600  REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1601  LOGICAL :: shells
1602 
1603  CHARACTER(len=*), PARAMETER :: routinen = 'debug_qmmm_forces_with_gauss_LG'
1604 
1605  INTEGER :: handle, i, igrid, indmm, j, k, ngrids
1606  REAL(kind=dp) :: coord_save
1607  REAL(kind=dp), DIMENSION(2) :: energy
1608  REAL(kind=dp), DIMENSION(3) :: err
1609  REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1610  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1611 
1612  ALLOCATE (num_forces(3, num_mm_atoms))
1613  CALL timeset(routinen, handle)
1614  ngrids = SIZE(aug_pools)
1615  CALL pw_pools_create_pws(aug_pools, grids)
1616  DO igrid = 1, ngrids
1617  CALL pw_zero(grids(igrid))
1618  END DO
1619  atoms: DO i = 1, num_mm_atoms
1620  indmm = mm_atom_index(i)
1621  coords: DO j = 1, 3
1622  coord_save = mm_particles(indmm)%r(j)
1623  energy = 0.0_dp
1624  diff: DO k = 1, 2
1625  mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1626  CALL pw_zero(grids(coarser_grid_level))
1627 
1628  CALL qmmm_elec_with_gaussian_lg(pgfs=pgfs, &
1629  cgrid=grids(coarser_grid_level), &
1630  mm_charges=mm_charges, &
1631  mm_atom_index=mm_atom_index, &
1632  mm_particles=mm_particles, &
1633  para_env=para_env, &
1634  per_potentials=per_potentials, &
1635  mm_cell=mm_cell, &
1636  dommoqm=dommoqm, &
1637  par_scheme=par_scheme, &
1638  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1639  shells=shells)
1640 
1641  energy(k) = pw_integral_ab(rho, grids(coarser_grid_level))
1642  END DO diff
1643  IF (iw > 0) &
1644  WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1645  "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1646  num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1647  mm_particles(indmm)%r(j) = coord_save
1648  END DO coords
1649  END DO atoms
1650 
1651  DO i = 1, num_mm_atoms
1652  indmm = mm_atom_index(i)
1653  err = 0.0_dp
1654  IF (all(num_forces /= 0.0_dp)) THEN
1655  err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1656  err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1657  err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1658  END IF
1659  IF (iw > 0) &
1660  WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1661  debug_force(2, i), num_forces(2, i), err(2), &
1662  debug_force(3, i), num_forces(3, i), err(3)
1663  cpassert(abs(err(1)) <= maxerr)
1664  cpassert(abs(err(2)) <= maxerr)
1665  cpassert(abs(err(3)) <= maxerr)
1666  END DO
1667 
1668  DEALLOCATE (num_forces)
1669  CALL pw_pools_give_back_pws(aug_pools, grids)
1670  CALL timestop(handle)
1671 100 FORMAT("MM Atom LR : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1672  END SUBROUTINE debug_qmmm_forces_with_gauss_lg
1673 
1674 ! **************************************************************************************************
1675 !> \brief Debugs qmmm_forces_with_gaussian_LR.. It may helps too ... ;-]
1676 !> \param pgfs ...
1677 !> \param aug_pools ...
1678 !> \param rho ...
1679 !> \param mm_charges ...
1680 !> \param mm_atom_index ...
1681 !> \param mm_particles ...
1682 !> \param num_mm_atoms ...
1683 !> \param coarser_grid_level ...
1684 !> \param potentials ...
1685 !> \param debug_force ...
1686 !> \param para_env ...
1687 !> \param mm_cell ...
1688 !> \param dOmmOqm ...
1689 !> \param iw ...
1690 !> \param par_scheme ...
1691 !> \param qmmm_spherical_cutoff ...
1692 !> \param shells ...
1693 !> \par History
1694 !> 08.2004 created [tlaino]
1695 !> \author Teodoro Laino
1696 ! **************************************************************************************************
1697  SUBROUTINE debug_qmmm_forces_with_gauss_lr(pgfs, aug_pools, rho, mm_charges, mm_atom_index, &
1698  mm_particles, num_mm_atoms, coarser_grid_level, potentials, &
1699  debug_force, para_env, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1700 
1701  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1702  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1703  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1704  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1705  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1706  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1707  INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level
1708  TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1709  REAL(kind=dp), DIMENSION(:, :) :: debug_force
1710  TYPE(mp_para_env_type), POINTER :: para_env
1711  TYPE(cell_type), POINTER :: mm_cell
1712  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1713  INTEGER, INTENT(IN) :: iw, par_scheme
1714  REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1715  LOGICAL :: shells
1716 
1717  CHARACTER(len=*), PARAMETER :: routinen = 'debug_qmmm_forces_with_gauss_LR'
1718 
1719  INTEGER :: handle, i, igrid, indmm, j, k, ngrids
1720  REAL(kind=dp) :: coord_save
1721  REAL(kind=dp), DIMENSION(2) :: energy
1722  REAL(kind=dp), DIMENSION(3) :: err
1723  REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1724  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1725 
1726  ALLOCATE (num_forces(3, num_mm_atoms))
1727  CALL timeset(routinen, handle)
1728  ngrids = SIZE(aug_pools)
1729  CALL pw_pools_create_pws(aug_pools, grids)
1730  DO igrid = 1, ngrids
1731  CALL pw_zero(grids(igrid))
1732  END DO
1733  atoms: DO i = 1, num_mm_atoms
1734  indmm = mm_atom_index(i)
1735  coords: DO j = 1, 3
1736  coord_save = mm_particles(indmm)%r(j)
1737  energy = 0.0_dp
1738  diff: DO k = 1, 2
1739  mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1740  CALL pw_zero(grids(coarser_grid_level))
1741 
1742  CALL qmmm_elec_with_gaussian_lr(pgfs=pgfs, &
1743  grid=grids(coarser_grid_level), &
1744  mm_charges=mm_charges, &
1745  mm_atom_index=mm_atom_index, &
1746  mm_particles=mm_particles, &
1747  para_env=para_env, &
1748  potentials=potentials, &
1749  mm_cell=mm_cell, &
1750  dommoqm=dommoqm, &
1751  par_scheme=par_scheme, &
1752  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1753  shells=shells)
1754 
1755  energy(k) = pw_integral_ab(rho, grids(coarser_grid_level))
1756  END DO diff
1757  IF (iw > 0) &
1758  WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1759  "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1760  num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1761  mm_particles(indmm)%r(j) = coord_save
1762  END DO coords
1763  END DO atoms
1764 
1765  DO i = 1, num_mm_atoms
1766  indmm = mm_atom_index(i)
1767  err = 0.0_dp
1768  IF (all(num_forces(:, i) /= 0.0_dp)) THEN
1769  err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1770  err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1771  err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1772  END IF
1773  IF (iw > 0) &
1774  WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1775  debug_force(2, i), num_forces(2, i), err(2), &
1776  debug_force(3, i), num_forces(3, i), err(3)
1777  cpassert(abs(err(1)) <= maxerr)
1778  cpassert(abs(err(2)) <= maxerr)
1779  cpassert(abs(err(3)) <= maxerr)
1780  END DO
1781 
1782  DEALLOCATE (num_forces)
1783  CALL pw_pools_give_back_pws(aug_pools, grids)
1784  CALL timestop(handle)
1785 100 FORMAT("MM Atom LR : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1786  END SUBROUTINE debug_qmmm_forces_with_gauss_lr
1787 
1788 END MODULE qmmm_gpw_forces
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_spline_utils
utils to manipulate splines on the regular grid of a pw
Definition: cp_spline_utils.F:14
cp_spline_utils::spline3_nopbc_interp
integer, parameter, public spline3_nopbc_interp
Definition: cp_spline_utils.F:38
cp_spline_utils::spline3_pbc_interp
integer, parameter, public spline3_pbc_interp
Definition: cp_spline_utils.F:38
cp_spline_utils::pw_restrict_s3
subroutine, public pw_restrict_s3(pw_fine_in, pw_coarse_out, coarse_pool, param_section)
restricts the function from a fine grid to a coarse one
Definition: cp_spline_utils.F:55
cube_utils
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
Definition: cube_utils.F:18
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_par_atom
integer, parameter, public do_par_atom
Definition: input_constants.F:286
input_constants::do_qmmm_none
integer, parameter, public do_qmmm_none
Definition: input_constants.F:384
input_constants::do_qmmm_pcharge
integer, parameter, public do_qmmm_pcharge
Definition: input_constants.F:384
input_constants::do_qmmm_coulomb
integer, parameter, public do_qmmm_coulomb
Definition: input_constants.F:384
input_constants::do_qmmm_swave
integer, parameter, public do_qmmm_swave
Definition: input_constants.F:384
input_constants::do_qmmm_gauss
integer, parameter, public do_qmmm_gauss
Definition: input_constants.F:384
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
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
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
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
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mm_collocate_potential
Calculate the MM potential by collocating the primitive Gaussian functions (pgf)
Definition: mm_collocate_potential.F:15
mm_collocate_potential::collocate_gf_rspace_nopbc
subroutine, public collocate_gf_rspace_nopbc(zetp, rp, scale, W, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo2, n_rep_real, mm_cell)
Main driver to collocate gaussian functions on grid without using periodic boundary conditions (NoPBC...
Definition: mm_collocate_potential.F:173
mm_collocate_potential::integrate_gf_rspace_nopbc
subroutine, public integrate_gf_rspace_nopbc(zetp, rp, scale, W, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo, force, n_rep_real, mm_cell)
Main driver to integrate gaussian functions on a grid function without using periodic boundary condit...
Definition: mm_collocate_potential.F:296
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qmmm_gaussian_types
Sets the typo for the gaussian treatment of the qm/mm interaction.
Definition: qmmm_gaussian_types.F:14
qmmm_gpw_energy
A collection of methods to treat the QM/MM electrostatic coupling.
Definition: qmmm_gpw_energy.F:14
qmmm_gpw_energy::qmmm_elec_with_gaussian
subroutine, public qmmm_elec_with_gaussian(qmmm_env, v_qmmm, mm_particles, aug_pools, cube_info, para_env, eps_mm_rspace, pw_pools, auxbas_grid, coarser_grid, interp_section, mm_cell)
Compute the QM/MM electrostatic Interaction collocating the gaussian Electrostatic Potential.
Definition: qmmm_gpw_energy.F:243
qmmm_gpw_energy::qmmm_elec_with_gaussian_lr
subroutine, public qmmm_elec_with_gaussian_lr(pgfs, grid, mm_charges, mm_atom_index, mm_particles, para_env, potentials, mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
Compute the QM/MM electrostatic Interaction collocating (1/R - Sum_NG Gaussians) on the coarser grid ...
Definition: qmmm_gpw_energy.F:757
qmmm_gpw_energy::qmmm_elec_with_gaussian_lg
subroutine, public qmmm_elec_with_gaussian_lg(pgfs, cgrid, mm_charges, mm_atom_index, mm_particles, para_env, per_potentials, mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
Compute the QM/MM electrostatic Interaction collocating (1/R - Sum_NG Gaussians) on the coarser grid ...
Definition: qmmm_gpw_energy.F:532
qmmm_gpw_forces
Routines to compute energy and forces in a QM/MM calculation.
Definition: qmmm_gpw_forces.F:14
qmmm_gpw_forces::qmmm_forces
subroutine, public qmmm_forces(qs_env, qmmm_env, mm_particles, calc_force, mm_cell)
General driver to Compute the contribution to the forces due to the QM/MM potential.
Definition: qmmm_gpw_forces.F:100
qmmm_se_forces
Calculation of the derivative of the QMMM Hamiltonian integral matrix <a|\sum_i q_i|b> for semi-empir...
Definition: qmmm_se_forces.F:13
qmmm_se_forces::deriv_se_qmmm_matrix
subroutine, public deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, Forces, Forces_added_charges)
Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian QMMM terms.
Definition: qmmm_se_forces.F:78
qmmm_tb_methods
TB methods used with QMMM.
Definition: qmmm_tb_methods.F:12
qmmm_tb_methods::deriv_tb_qmmm_matrix_pc
subroutine, public deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, Forces, Forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition: qmmm_tb_methods.F:830
qmmm_tb_methods::deriv_tb_qmmm_matrix
subroutine, public deriv_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, Forces, Forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition: qmmm_tb_methods.F:593
qmmm_types_low
Definition: qmmm_types_low.F:13
qmmm_util
Definition: qmmm_util.F:13
qmmm_util::spherical_cutoff_factor
subroutine, public spherical_cutoff_factor(spherical_cutoff, rij, factor)
Computes a spherical cutoff factor for the QMMM interactions.
Definition: qmmm_util.F:615
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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.
Definition: qs_environment_types.F:502
qs_ks_qmmm_types
Definition: qs_ks_qmmm_types.F:13
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229