(git:ccc2433)
qmmm_gpw_energy.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 A collection of methods to treat the QM/MM electrostatic coupling
10 !> \par History
11 !> 5.2004 created [tlaino]
12 !> \author Teodoro Laino
13 ! **************************************************************************************************
14 MODULE qmmm_gpw_energy
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_p_file,&
21  cp_print_key_finished_output,&
22  cp_print_key_should_output,&
23  cp_print_key_unit_nr
24  USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
25  USE cp_spline_utils, ONLY: pw_prolongate_s3,&
26  spline3_nopbc_interp,&
27  spline3_pbc_interp
28  USE cube_utils, ONLY: cube_info_type
29  USE input_constants, ONLY: do_par_atom,&
30  do_qmmm_coulomb,&
31  do_qmmm_gauss,&
32  do_qmmm_none,&
33  do_qmmm_pcharge,&
34  do_qmmm_swave
35  USE input_section_types, ONLY: section_get_ivals,&
36  section_vals_get_subs_vals,&
37  section_vals_type,&
38  section_vals_val_get
39  USE kinds, ONLY: dp
40  USE message_passing, ONLY: mp_para_env_type
41  USE mm_collocate_potential, ONLY: collocate_gf_rspace_nopbc
42  USE particle_list_types, ONLY: particle_list_type
43  USE particle_types, ONLY: particle_type
44  USE pw_env_types, ONLY: pw_env_get,&
45  pw_env_type
46  USE pw_methods, ONLY: pw_zero
47  USE pw_pool_types, ONLY: pw_pool_p_type,&
48  pw_pools_create_pws,&
49  pw_pools_give_back_pws
50  USE pw_types, ONLY: pw_r3d_rs_type
51  USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,&
52  qmmm_gaussian_type
53  USE qmmm_se_energy, ONLY: build_se_qmmm_matrix
54  USE qmmm_tb_methods, ONLY: build_tb_qmmm_matrix,&
55  build_tb_qmmm_matrix_pc,&
56  build_tb_qmmm_matrix_zero
57  USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
58  qmmm_per_pot_p_type,&
59  qmmm_per_pot_type,&
60  qmmm_pot_p_type,&
61  qmmm_pot_type
62  USE qmmm_util, ONLY: spherical_cutoff_factor
63  USE qs_environment_types, ONLY: get_qs_env,&
64  qs_environment_type
65  USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
66  USE qs_subsys_types, ONLY: qs_subsys_get,&
67  qs_subsys_type
68 #include "./base/base_uses.f90"
69 
70  IMPLICIT NONE
71  PRIVATE
72 
73  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
74  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gpw_energy'
75 
76  PUBLIC :: qmmm_el_coupling
77  PUBLIC :: qmmm_elec_with_gaussian, &
78  qmmm_elec_with_gaussian_lr, qmmm_elec_with_gaussian_lg
79 !***
80 CONTAINS
81 
82 ! **************************************************************************************************
83 !> \brief Main Driver to compute the QM/MM Electrostatic Coupling
84 !> \param qs_env ...
85 !> \param qmmm_env ...
86 !> \param mm_particles ...
87 !> \param mm_cell ...
88 !> \par History
89 !> 05.2004 created [tlaino]
90 !> \author Teodoro Laino
91 ! **************************************************************************************************
92  SUBROUTINE qmmm_el_coupling(qs_env, qmmm_env, mm_particles, mm_cell)
93  TYPE(qs_environment_type), POINTER :: qs_env
94  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
95  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
96  TYPE(cell_type), POINTER :: mm_cell
97 
98  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_el_coupling'
99 
100  INTEGER :: handle, iw, iw2
101  LOGICAL :: mpi_io
102  TYPE(cp_logger_type), POINTER :: logger
103  TYPE(dft_control_type), POINTER :: dft_control
104  TYPE(mp_para_env_type), POINTER :: para_env
105  TYPE(particle_list_type), POINTER :: particles
106  TYPE(pw_env_type), POINTER :: pw_env
107  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
108  TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
109  TYPE(qs_subsys_type), POINTER :: subsys
110  TYPE(section_vals_type), POINTER :: input_section, interp_section, &
111  print_section
112 
113  CALL timeset(routinen, handle)
114  logger => cp_get_default_logger()
115  NULLIFY (ks_qmmm_env_loc, pw_pools, pw_env, input_section, dft_control)
116  CALL get_qs_env(qs_env=qs_env, &
117  pw_env=pw_env, &
118  para_env=para_env, &
119  input=input_section, &
120  ks_qmmm_env=ks_qmmm_env_loc, &
121  subsys=subsys, &
122  dft_control=dft_control)
123  CALL qs_subsys_get(subsys, particles=particles)
124 
125  CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools)
126  print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
127  iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
128  extension=".qmmmLog")
129  IF (iw > 0) &
130  WRITE (iw, '(T2,"QMMM|",1X,A)') "Information on the QM/MM Electrostatic Potential:"
131  !
132  ! Initializing vectors:
133  ! Zeroing v_qmmm_rspace
134  CALL pw_zero(ks_qmmm_env_loc%v_qmmm_rspace)
135  IF (dft_control%qs_control%semi_empirical) THEN
136  ! SEMIEMPIRICAL
137  SELECT CASE (qmmm_env%qmmm_coupl_type)
138  CASE (do_qmmm_coulomb, do_qmmm_none)
139  CALL build_se_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
140  IF (qmmm_env%qmmm_coupl_type == do_qmmm_none) THEN
141  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
142  "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
143  END IF
144  CASE (do_qmmm_pcharge)
145  cpabort("Point charge QM/MM electrostatic coupling not yet implemented for SE.")
146  CASE (do_qmmm_gauss, do_qmmm_swave)
147  cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.")
148  CASE DEFAULT
149  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
150  cpabort("")
151  END SELECT
152  ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
153  ! DFTB
154  SELECT CASE (qmmm_env%qmmm_coupl_type)
155  CASE (do_qmmm_none)
156  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
157  "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
158  CALL build_tb_qmmm_matrix_zero(qs_env, para_env)
159  CASE (do_qmmm_coulomb)
160  CALL build_tb_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
161  CASE (do_qmmm_pcharge)
162  CALL build_tb_qmmm_matrix_pc(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
163  CASE (do_qmmm_gauss, do_qmmm_swave)
164  cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not implemented for DFTB.")
165  CASE DEFAULT
166  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
167  cpabort("")
168  END SELECT
169  ELSE
170  ! QS
171  SELECT CASE (qmmm_env%qmmm_coupl_type)
172  CASE (do_qmmm_coulomb)
173  cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
174  CASE (do_qmmm_pcharge)
175  cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
176  CASE (do_qmmm_gauss, do_qmmm_swave)
177  IF (iw > 0) &
178  WRITE (iw, '(T2,"QMMM|",1X,A)') &
179  "QM/MM Coupling computed collocating the Gaussian Potential Functions."
180  interp_section => section_vals_get_subs_vals(input_section, &
181  "QMMM%INTERPOLATOR")
182  CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env, &
183  v_qmmm=ks_qmmm_env_loc%v_qmmm_rspace, &
184  mm_particles=mm_particles, &
185  aug_pools=qmmm_env%aug_pools, &
186  para_env=para_env, &
187  eps_mm_rspace=qmmm_env%eps_mm_rspace, &
188  cube_info=ks_qmmm_env_loc%cube_info, &
189  pw_pools=pw_pools, &
190  auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
191  coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
192  interp_section=interp_section, &
193  mm_cell=mm_cell)
194  CASE (do_qmmm_none)
195  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
196  "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
197  CASE DEFAULT
198  IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
199  cpabort("")
200  END SELECT
201  ! Dump info on the electrostatic potential if requested
202  IF (btest(cp_print_key_should_output(logger%iter_info, print_section, &
203  "POTENTIAL"), cp_p_file)) THEN
204  mpi_io = .true.
205  iw2 = cp_print_key_unit_nr(logger, print_section, "POTENTIAL", &
206  extension=".qmmmLog", mpi_io=mpi_io)
207  CALL cp_pw_to_cube(ks_qmmm_env_loc%v_qmmm_rspace, iw2, &
208  particles=particles, &
209  stride=section_get_ivals(print_section, "POTENTIAL%STRIDE"), &
210  title="QM/MM: MM ELECTROSTATIC POTENTIAL ", &
211  mpi_io=mpi_io)
212  CALL cp_print_key_finished_output(iw2, logger, print_section, &
213  "POTENTIAL", mpi_io=mpi_io)
214  END IF
215  END IF
216  CALL cp_print_key_finished_output(iw, logger, print_section, &
217  "PROGRAM_RUN_INFO")
218  CALL timestop(handle)
219  END SUBROUTINE qmmm_el_coupling
220 
221 ! **************************************************************************************************
222 !> \brief Compute the QM/MM electrostatic Interaction collocating the gaussian
223 !> Electrostatic Potential
224 !> \param qmmm_env ...
225 !> \param v_qmmm ...
226 !> \param mm_particles ...
227 !> \param aug_pools ...
228 !> \param cube_info ...
229 !> \param para_env ...
230 !> \param eps_mm_rspace ...
231 !> \param pw_pools ...
232 !> \param auxbas_grid ...
233 !> \param coarser_grid ...
234 !> \param interp_section ...
235 !> \param mm_cell ...
236 !> \par History
237 !> 06.2004 created [tlaino]
238 !> \author Teodoro Laino
239 ! **************************************************************************************************
240  SUBROUTINE qmmm_elec_with_gaussian(qmmm_env, v_qmmm, mm_particles, &
241  aug_pools, cube_info, para_env, eps_mm_rspace, pw_pools, &
242  auxbas_grid, coarser_grid, interp_section, mm_cell)
243  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
244  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_qmmm
245  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
246  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
247  TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
248  TYPE(mp_para_env_type), POINTER :: para_env
249  REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
250  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
251  INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
252  TYPE(section_vals_type), POINTER :: interp_section
253  TYPE(cell_type), POINTER :: mm_cell
254 
255  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian'
256 
257  INTEGER :: handle, handle2, igrid, ilevel, &
258  kind_interp, lb(3), ngrids, ub(3)
259  LOGICAL :: shells
260  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
261 
262  cpassert(ASSOCIATED(mm_particles))
263  cpassert(ASSOCIATED(qmmm_env%mm_atom_chrg))
264  cpassert(ASSOCIATED(qmmm_env%mm_atom_index))
265  cpassert(ASSOCIATED(aug_pools))
266  cpassert(ASSOCIATED(pw_pools))
267  !Statements
268  CALL timeset(routinen, handle)
269  ngrids = SIZE(pw_pools)
270  CALL pw_pools_create_pws(aug_pools, grids)
271  DO igrid = 1, ngrids
272  CALL pw_zero(grids(igrid))
273  END DO
274 
275  shells = .false.
276 
277  CALL qmmm_elec_with_gaussian_low(grids, mm_particles, &
278  qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
279  cube_info, para_env, eps_mm_rspace, qmmm_env%pgfs, &
280  auxbas_grid, coarser_grid, qmmm_env%potentials, &
281  mm_cell=mm_cell, dommoqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
282  per_potentials=qmmm_env%per_potentials, par_scheme=qmmm_env%par_scheme, &
283  qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)
284 
285  IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
286  CALL qmmm_elec_with_gaussian_low(grids, qmmm_env%added_charges%added_particles, &
287  qmmm_env%added_charges%mm_atom_chrg, &
288  qmmm_env%added_charges%mm_atom_index, &
289  cube_info, para_env, eps_mm_rspace, qmmm_env%added_charges%pgfs, auxbas_grid, &
290  coarser_grid, qmmm_env%added_charges%potentials, &
291  mm_cell=mm_cell, dommoqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
292  per_potentials=qmmm_env%added_charges%per_potentials, par_scheme=qmmm_env%par_scheme, &
293  qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)
294  END IF
295  IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN
296  shells = .true.
297  CALL qmmm_elec_with_gaussian_low(grids, qmmm_env%added_shells%added_particles, &
298  qmmm_env%added_shells%mm_core_chrg, &
299  qmmm_env%added_shells%mm_core_index, &
300  cube_info, para_env, eps_mm_rspace, qmmm_env%added_shells%pgfs, auxbas_grid, &
301  coarser_grid, qmmm_env%added_shells%potentials, &
302  mm_cell=mm_cell, dommoqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
303  per_potentials=qmmm_env%added_shells%per_potentials, &
304  par_scheme=qmmm_env%par_scheme, qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, &
305  shells=shells)
306  END IF
307  ! Sumup all contributions according the parallelization scheme
308  IF (qmmm_env%par_scheme == do_par_atom) THEN
309  DO ilevel = 1, SIZE(grids)
310  CALL para_env%sum(grids(ilevel)%array)
311  END DO
312  END IF
313  ! RealSpace Interpolation
314  CALL section_vals_val_get(interp_section, "kind", i_val=kind_interp)
315  SELECT CASE (kind_interp)
316  CASE (spline3_nopbc_interp, spline3_pbc_interp)
317  ! Spline Iterpolator
318  CALL para_env%sync()
319  CALL timeset(trim(routinen)//":spline3Int", handle2)
320  DO ilevel = coarser_grid, auxbas_grid + 1, -1
321  CALL pw_prolongate_s3(grids(ilevel), &
322  grids(ilevel - 1), &
323  aug_pools(ilevel)%pool, &
324  param_section=interp_section)
325  END DO
326  CALL timestop(handle2)
327  CASE DEFAULT
328  cpabort("")
329  END SELECT
330  lb = v_qmmm%pw_grid%bounds_local(1, :)
331  ub = v_qmmm%pw_grid%bounds_local(2, :)
332 
333  v_qmmm%array = grids(auxbas_grid)%array(lb(1):ub(1), &
334  lb(2):ub(2), &
335  lb(3):ub(3))
336 
337  CALL pw_pools_give_back_pws(aug_pools, grids)
338 
339  CALL timestop(handle)
340  END SUBROUTINE qmmm_elec_with_gaussian
341 
342 ! **************************************************************************************************
343 !> \brief Compute the QM/MM electrostatic Interaction collocating the gaussian
344 !> Electrostatic Potential - Low Level
345 !> \param tmp_grid ...
346 !> \param mm_particles ...
347 !> \param mm_charges ...
348 !> \param mm_atom_index ...
349 !> \param cube_info ...
350 !> \param para_env ...
351 !> \param eps_mm_rspace ...
352 !> \param pgfs ...
353 !> \param auxbas_grid ...
354 !> \param coarser_grid ...
355 !> \param potentials ...
356 !> \param mm_cell ...
357 !> \param dOmmOqm ...
358 !> \param periodic ...
359 !> \param per_potentials ...
360 !> \param par_scheme ...
361 !> \param qmmm_spherical_cutoff ...
362 !> \param shells ...
363 !> \par History
364 !> 06.2004 created [tlaino]
365 !> \author Teodoro Laino
366 ! **************************************************************************************************
367  SUBROUTINE qmmm_elec_with_gaussian_low(tmp_grid, mm_particles, mm_charges, &
368  mm_atom_index, cube_info, para_env, &
369  eps_mm_rspace, pgfs, auxbas_grid, coarser_grid, &
370  potentials, mm_cell, dOmmOqm, periodic, per_potentials, par_scheme, &
371  qmmm_spherical_cutoff, shells)
372  TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: tmp_grid
373  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
374  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
375  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
376  TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
377  TYPE(mp_para_env_type), POINTER :: para_env
378  REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
379  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
380  INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
381  TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
382  TYPE(cell_type), POINTER :: mm_cell
383  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
384  LOGICAL, INTENT(IN) :: periodic
385  TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
386  INTEGER, INTENT(IN) :: par_scheme
387  REAL(kind=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
388  LOGICAL, INTENT(IN) :: shells
389 
390  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian_low', &
391  routinenb = 'qmmm_elec_gaussian_low'
392 
393  INTEGER :: handle, handle2, igauss, ilevel, imm, &
394  indmm, iradtyp, lindmm, myind, &
395  n_rep_real(3)
396  INTEGER, DIMENSION(2, 3) :: bo2
397  REAL(kind=dp) :: alpha, height, sph_chrg_factor, w
398  REAL(kind=dp), DIMENSION(3) :: ra
399  REAL(kind=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat
400  TYPE(qmmm_gaussian_type), POINTER :: pgf
401  TYPE(qmmm_per_pot_type), POINTER :: per_pot
402  TYPE(qmmm_pot_type), POINTER :: pot
403 
404  NULLIFY (pgf, pot, per_pot, xdat, ydat, zdat)
405  CALL timeset(routinen, handle)
406  CALL timeset(routinenb//"_G", handle2)
407  bo2 = tmp_grid(auxbas_grid)%pw_grid%bounds
408  ALLOCATE (xdat(bo2(1, 1):bo2(2, 1)))
409  ALLOCATE (ydat(bo2(1, 2):bo2(2, 2)))
410  ALLOCATE (zdat(bo2(1, 3):bo2(2, 3)))
411  IF (par_scheme == do_par_atom) myind = 0
412  radius: DO iradtyp = 1, SIZE(pgfs)
413  pgf => pgfs(iradtyp)%pgf
414  pot => potentials(iradtyp)%pot
415  n_rep_real = 0
416  IF (periodic) THEN
417  per_pot => per_potentials(iradtyp)%pot
418  n_rep_real = per_pot%n_rep_real
419  END IF
420  gaussian: DO igauss = 1, pgf%Number_of_Gaussians
421  alpha = 1.0_dp/pgf%Gk(igauss)
422  alpha = alpha*alpha
423  height = pgf%Ak(igauss)
424  ilevel = pgf%grid_level(igauss)
425  atoms: DO imm = 1, SIZE(pot%mm_atom_index)
426  IF (par_scheme == do_par_atom) THEN
427  myind = myind + 1
428  IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
429  END IF
430  lindmm = pot%mm_atom_index(imm)
431  indmm = mm_atom_index(lindmm)
432  IF (shells) THEN
433  ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
434  ELSE
435  ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
436  END IF
437  w = mm_charges(lindmm)*height
438  ! Possible Spherical Cutoff
439  IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
440  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
441  w = w*sph_chrg_factor
442  END IF
443  IF (abs(w) <= epsilon(0.0_dp)) cycle atoms
444  CALL collocate_gf_rspace_nopbc(zetp=alpha, &
445  rp=ra, &
446  scale=-1.0_dp, &
447  w=w, &
448  pwgrid=tmp_grid(ilevel), &
449  cube_info=cube_info(ilevel), &
450  eps_mm_rspace=eps_mm_rspace, &
451  xdat=xdat, &
452  ydat=ydat, &
453  zdat=zdat, &
454  bo2=bo2, &
455  n_rep_real=n_rep_real, &
456  mm_cell=mm_cell)
457  END DO atoms
458  END DO gaussian
459  END DO radius
460  IF (ASSOCIATED(xdat)) THEN
461  DEALLOCATE (xdat)
462  END IF
463  IF (ASSOCIATED(ydat)) THEN
464  DEALLOCATE (ydat)
465  END IF
466  IF (ASSOCIATED(zdat)) THEN
467  DEALLOCATE (zdat)
468  END IF
469  CALL timestop(handle2)
470  CALL timeset(routinenb//"_R", handle2)
471  IF (periodic) THEN
472  ! Long Range Part of the QM/MM Potential with Gaussians With Periodic Boundary Conditions
473  CALL qmmm_elec_with_gaussian_lg(pgfs=pgfs, &
474  cgrid=tmp_grid(coarser_grid), &
475  mm_charges=mm_charges, &
476  mm_atom_index=mm_atom_index, &
477  mm_particles=mm_particles, &
478  para_env=para_env, &
479  per_potentials=per_potentials, &
480  mm_cell=mm_cell, &
481  dommoqm=dommoqm, &
482  par_scheme=par_scheme, &
483  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
484  shells=shells)
485  ELSE
486  ! Long Range Part of the QM/MM Potential with Gaussians
487  CALL qmmm_elec_with_gaussian_lr(pgfs=pgfs, &
488  grid=tmp_grid(coarser_grid), &
489  mm_charges=mm_charges, &
490  mm_atom_index=mm_atom_index, &
491  mm_particles=mm_particles, &
492  para_env=para_env, &
493  potentials=potentials, &
494  mm_cell=mm_cell, &
495  dommoqm=dommoqm, &
496  par_scheme=par_scheme, &
497  qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
498  shells=shells)
499  END IF
500  CALL timestop(handle2)
501  CALL timestop(handle)
502 
503  END SUBROUTINE qmmm_elec_with_gaussian_low
504 
505 ! **************************************************************************************************
506 !> \brief Compute the QM/MM electrostatic Interaction collocating
507 !> (1/R - Sum_NG Gaussians) on the coarser grid level in G-SPACE
508 !> Long Range QM/MM Electrostatic Potential with Gaussian - Low Level
509 !> PERIODIC BOUNDARY CONDITION VERSION
510 !> \param pgfs ...
511 !> \param cgrid ...
512 !> \param mm_charges ...
513 !> \param mm_atom_index ...
514 !> \param mm_particles ...
515 !> \param para_env ...
516 !> \param per_potentials ...
517 !> \param mm_cell ...
518 !> \param dOmmOqm ...
519 !> \param par_scheme ...
520 !> \param qmmm_spherical_cutoff ...
521 !> \param shells ...
522 !> \par History
523 !> 07.2004 created [tlaino]
524 !> \author Teodoro Laino
525 !> \note
526 !> This version includes the explicit code of Eval_Interp_Spl3_pbc
527 !> in order to achieve better performance
528 ! **************************************************************************************************
529  SUBROUTINE qmmm_elec_with_gaussian_lg(pgfs, cgrid, mm_charges, mm_atom_index, &
530  mm_particles, para_env, per_potentials, &
531  mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
532  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
533  TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
534  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
535  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
536  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
537  TYPE(mp_para_env_type), POINTER :: para_env
538  TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
539  TYPE(cell_type), POINTER :: mm_cell
540  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
541  INTEGER, INTENT(IN) :: par_scheme
542  REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
543  LOGICAL :: shells
544 
545  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian_LG'
546 
547  INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, &
548  ij3, ij4, ik1, ik2, ik3, ik4, imm, &
549  indmm, iradtyp, ivec(3), j, k, lindmm, &
550  my_j, my_k, myind, npts(3)
551  INTEGER, DIMENSION(2, 3) :: bo, gbo
552  REAL(kind=dp) :: a1, a2, a3, abc_x(4, 4), abc_x_y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, &
553  dr1, dr1c, dr2, dr2c, dr3, dr3c, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, &
554  p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, rt3, rv1, rv2, rv3, s1, s2, s3, s4, &
555  sph_chrg_factor, t1, t2, t3, t4, u1, u2, u3, v1, v2, v3, v4, val, xd1, xd2, xd3, xs1, &
556  xs2, xs3
557  REAL(kind=dp), DIMENSION(3) :: ra, vec
558  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid, grid2
559  TYPE(pw_r3d_rs_type), POINTER :: pw
560  TYPE(qmmm_per_pot_type), POINTER :: per_pot
561 
562  CALL timeset(routinen, handle)
563  NULLIFY (grid, pw)
564  dr1c = cgrid%pw_grid%dr(1)
565  dr2c = cgrid%pw_grid%dr(2)
566  dr3c = cgrid%pw_grid%dr(3)
567  gbo = cgrid%pw_grid%bounds
568  bo = cgrid%pw_grid%bounds_local
569  grid2 => cgrid%array
570  IF (par_scheme == do_par_atom) myind = 0
571  radius: DO iradtyp = 1, SIZE(pgfs)
572  per_pot => per_potentials(iradtyp)%pot
573  pw => per_pot%TabLR
574  npts = pw%pw_grid%npts
575  dr1 = pw%pw_grid%dr(1)
576  dr2 = pw%pw_grid%dr(2)
577  dr3 = pw%pw_grid%dr(3)
578  grid => pw%array(:, :, :)
579  !$OMP PARALLEL DO DEFAULT(NONE) &
580  !$OMP SHARED(bo, gbo, grid, grid2, pw, npts, per_pot, mm_atom_index) &
581  !$OMP SHARED(dr1, dr2, dr3, dr1c, dr2c, dr3c, par_scheme, mm_charges, mm_particles) &
582  !$OMP SHARED(mm_cell, dOmmOqm, shells, para_env, IRadTyp, qmmm_spherical_cutoff) &
583  !$OMP PRIVATE(Imm, LIndMM, IndMM, qt, sph_chrg_factor, ra, myind) &
584  !$OMP PRIVATE(rt1, rt2, rt3, k, vec, ivec, xd1, xd2, xd3, ik1, ik2, ik3, ik4) &
585  !$OMP PRIVATE(ij1, ij2, ij3, ij4, ii1, ii2, ii3, ii4, my_k, my_j, xs1, xs2, xs3) &
586  !$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, v1, v2, v3, v4, e1, e2, e3) &
587  !$OMP PRIVATE(f1, f2, f3, g1, g2, g3, h1, h2, h3, s1, s2, s3, s4, a1, a2, a3) &
588  !$OMP PRIVATE(b1, b2, b3, c1, c2, c3, d1, d2, d3, t1, t2, t3, t4, u1, u2, u3, val) &
589  !$OMP PRIVATE(rv1, rv2, rv3, abc_X, abc_X_Y)
590  atoms: DO imm = 1, SIZE(per_pot%mm_atom_index)
591  IF (par_scheme == do_par_atom) THEN
592  myind = imm + (iradtyp - 1)*SIZE(per_pot%mm_atom_index)
593  IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle
594  END IF
595  lindmm = per_pot%mm_atom_index(imm)
596  indmm = mm_atom_index(lindmm)
597  qt = mm_charges(lindmm)
598  IF (shells) THEN
599  ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
600  ELSE
601  ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
602  END IF
603  ! Possible Spherical Cutoff
604  IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
605  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
606  qt = qt*sph_chrg_factor
607  END IF
608  IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
609  rt1 = ra(1)
610  rt2 = ra(2)
611  rt3 = ra(3)
612  loopongrid: DO k = bo(1, 3), bo(2, 3)
613  my_k = k - gbo(1, 3)
614  xs3 = real(my_k, dp)*dr3c
615  my_j = bo(1, 2) - gbo(1, 2)
616  xs2 = real(my_j, dp)*dr2c
617  rv3 = rt3 - xs3
618  vec(3) = rv3
619  ivec(3) = floor(vec(3)/pw%pw_grid%dr(3))
620  xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
621  ik1 = modulo(ivec(3) - 1, npts(3)) + 1
622  ik2 = modulo(ivec(3), npts(3)) + 1
623  ik3 = modulo(ivec(3) + 1, npts(3)) + 1
624  ik4 = modulo(ivec(3) + 2, npts(3)) + 1
625  p1 = 3.0_dp + xd3
626  p2 = p1*p1
627  p3 = p2*p1
628  q1 = 2.0_dp + xd3
629  q2 = q1*q1
630  q3 = q2*q1
631  r1 = 1.0_dp + xd3
632  r2 = r1*r1
633  r3 = r2*r1
634  u1 = xd3
635  u2 = u1*u1
636  u3 = u2*u1
637  v1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
638  v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
639  v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
640  v4 = 1.0_dp/6.0_dp*u3
641  DO j = bo(1, 2), bo(2, 2)
642  xs1 = (bo(1, 1) - gbo(1, 1))*dr1c
643  rv2 = rt2 - xs2
644  vec(2) = rv2
645  ivec(2) = floor(vec(2)/pw%pw_grid%dr(2))
646  xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
647  ij1 = modulo(ivec(2) - 1, npts(2)) + 1
648  ij2 = modulo(ivec(2), npts(2)) + 1
649  ij3 = modulo(ivec(2) + 1, npts(2)) + 1
650  ij4 = modulo(ivec(2) + 2, npts(2)) + 1
651  e1 = 3.0_dp + xd2
652  e2 = e1*e1
653  e3 = e2*e1
654  f1 = 2.0_dp + xd2
655  f2 = f1*f1
656  f3 = f2*f1
657  g1 = 1.0_dp + xd2
658  g2 = g1*g1
659  g3 = g2*g1
660  h1 = xd2
661  h2 = h1*h1
662  h3 = h2*h1
663  s1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
664  s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
665  s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
666  s4 = 1.0_dp/6.0_dp*h3
667  DO i = bo(1, 1), bo(2, 1)
668  rv1 = rt1 - xs1
669  vec(1) = rv1
670  ivec(1) = floor(vec(1)/pw%pw_grid%dr(1))
671  xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
672  ii1 = modulo(ivec(1) - 1, npts(1)) + 1
673  ii2 = modulo(ivec(1), npts(1)) + 1
674  ii3 = modulo(ivec(1) + 1, npts(1)) + 1
675  ii4 = modulo(ivec(1) + 2, npts(1)) + 1
676  !
677  ! Spline Interpolation
678  !
679 
680  a1 = 3.0_dp + xd1
681  a2 = a1*a1
682  a3 = a2*a1
683  b1 = 2.0_dp + xd1
684  b2 = b1*b1
685  b3 = b2*b1
686  c1 = 1.0_dp + xd1
687  c2 = c1*c1
688  c3 = c2*c1
689  d1 = xd1
690  d2 = d1*d1
691  d3 = d2*d1
692  t1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
693  t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
694  t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
695  t4 = 1.0_dp/6.0_dp*d3
696 
697  abc_x(1, 1) = grid(ii1, ij1, ik1)*v1 + grid(ii1, ij1, ik2)*v2 + grid(ii1, ij1, ik3)*v3 + grid(ii1, ij1, ik4)*v4
698  abc_x(1, 2) = grid(ii1, ij2, ik1)*v1 + grid(ii1, ij2, ik2)*v2 + grid(ii1, ij2, ik3)*v3 + grid(ii1, ij2, ik4)*v4
699  abc_x(1, 3) = grid(ii1, ij3, ik1)*v1 + grid(ii1, ij3, ik2)*v2 + grid(ii1, ij3, ik3)*v3 + grid(ii1, ij3, ik4)*v4
700  abc_x(1, 4) = grid(ii1, ij4, ik1)*v1 + grid(ii1, ij4, ik2)*v2 + grid(ii1, ij4, ik3)*v3 + grid(ii1, ij4, ik4)*v4
701  abc_x(2, 1) = grid(ii2, ij1, ik1)*v1 + grid(ii2, ij1, ik2)*v2 + grid(ii2, ij1, ik3)*v3 + grid(ii2, ij1, ik4)*v4
702  abc_x(2, 2) = grid(ii2, ij2, ik1)*v1 + grid(ii2, ij2, ik2)*v2 + grid(ii2, ij2, ik3)*v3 + grid(ii2, ij2, ik4)*v4
703  abc_x(2, 3) = grid(ii2, ij3, ik1)*v1 + grid(ii2, ij3, ik2)*v2 + grid(ii2, ij3, ik3)*v3 + grid(ii2, ij3, ik4)*v4
704  abc_x(2, 4) = grid(ii2, ij4, ik1)*v1 + grid(ii2, ij4, ik2)*v2 + grid(ii2, ij4, ik3)*v3 + grid(ii2, ij4, ik4)*v4
705  abc_x(3, 1) = grid(ii3, ij1, ik1)*v1 + grid(ii3, ij1, ik2)*v2 + grid(ii3, ij1, ik3)*v3 + grid(ii3, ij1, ik4)*v4
706  abc_x(3, 2) = grid(ii3, ij2, ik1)*v1 + grid(ii3, ij2, ik2)*v2 + grid(ii3, ij2, ik3)*v3 + grid(ii3, ij2, ik4)*v4
707  abc_x(3, 3) = grid(ii3, ij3, ik1)*v1 + grid(ii3, ij3, ik2)*v2 + grid(ii3, ij3, ik3)*v3 + grid(ii3, ij3, ik4)*v4
708  abc_x(3, 4) = grid(ii3, ij4, ik1)*v1 + grid(ii3, ij4, ik2)*v2 + grid(ii3, ij4, ik3)*v3 + grid(ii3, ij4, ik4)*v4
709  abc_x(4, 1) = grid(ii4, ij1, ik1)*v1 + grid(ii4, ij1, ik2)*v2 + grid(ii4, ij1, ik3)*v3 + grid(ii4, ij1, ik4)*v4
710  abc_x(4, 2) = grid(ii4, ij2, ik1)*v1 + grid(ii4, ij2, ik2)*v2 + grid(ii4, ij2, ik3)*v3 + grid(ii4, ij2, ik4)*v4
711  abc_x(4, 3) = grid(ii4, ij3, ik1)*v1 + grid(ii4, ij3, ik2)*v2 + grid(ii4, ij3, ik3)*v3 + grid(ii4, ij3, ik4)*v4
712  abc_x(4, 4) = grid(ii4, ij4, ik1)*v1 + grid(ii4, ij4, ik2)*v2 + grid(ii4, ij4, ik3)*v3 + grid(ii4, ij4, ik4)*v4
713 
714  abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
715  abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
716  abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
717  abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
718 
719  val = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
720  !$OMP ATOMIC
721  grid2(i, j, k) = grid2(i, j, k) - val*qt
722  !$OMP END ATOMIC
723  xs1 = xs1 + dr1c
724  END DO
725  xs2 = xs2 + dr2c
726  END DO
727  END DO loopongrid
728  END DO atoms
729  !$OMP END PARALLEL DO
730  END DO radius
731  CALL timestop(handle)
732  END SUBROUTINE qmmm_elec_with_gaussian_lg
733 
734 ! **************************************************************************************************
735 !> \brief Compute the QM/MM electrostatic Interaction collocating
736 !> (1/R - Sum_NG Gaussians) on the coarser grid level.
737 !> Long Range QM/MM Electrostatic Potential with Gaussian - Low Level
738 !> \param pgfs ...
739 !> \param grid ...
740 !> \param mm_charges ...
741 !> \param mm_atom_index ...
742 !> \param mm_particles ...
743 !> \param para_env ...
744 !> \param potentials ...
745 !> \param mm_cell ...
746 !> \param dOmmOqm ...
747 !> \param par_scheme ...
748 !> \param qmmm_spherical_cutoff ...
749 !> \param shells ...
750 !> \par History
751 !> 07.2004 created [tlaino]
752 !> \author Teodoro Laino
753 ! **************************************************************************************************
754  SUBROUTINE qmmm_elec_with_gaussian_lr(pgfs, grid, mm_charges, mm_atom_index, &
755  mm_particles, para_env, potentials, &
756  mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
757  TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
758  TYPE(pw_r3d_rs_type), INTENT(IN) :: grid
759  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
760  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
761  TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
762  TYPE(mp_para_env_type), POINTER :: para_env
763  TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
764  TYPE(cell_type), POINTER :: mm_cell
765  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
766  INTEGER, INTENT(IN) :: par_scheme
767  REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
768  LOGICAL :: shells
769 
770  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian_LR'
771 
772  INTEGER :: handle, i, imm, indmm, iradtyp, ix, j, &
773  k, lindmm, my_j, my_k, myind, n1, n2, &
774  n3
775  INTEGER, DIMENSION(2, 3) :: bo, gbo
776  REAL(kind=dp) :: dr1, dr2, dr3, dx, qt, r, r2, rt1, rt2, &
777  rt3, rv1, rv2, rv3, rx, rx2, rx3, &
778  sph_chrg_factor, term, xs1, xs2, xs3
779  REAL(kind=dp), DIMENSION(3) :: ra
780  REAL(kind=dp), DIMENSION(:, :), POINTER :: pot0_2
781  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid2
782  TYPE(qmmm_pot_type), POINTER :: pot
783 
784  CALL timeset(routinen, handle)
785  n1 = grid%pw_grid%npts(1)
786  n2 = grid%pw_grid%npts(2)
787  n3 = grid%pw_grid%npts(3)
788  dr1 = grid%pw_grid%dr(1)
789  dr2 = grid%pw_grid%dr(2)
790  dr3 = grid%pw_grid%dr(3)
791  gbo = grid%pw_grid%bounds
792  bo = grid%pw_grid%bounds_local
793  grid2 => grid%array
794  IF (par_scheme == do_par_atom) myind = 0
795  radius: DO iradtyp = 1, SIZE(pgfs)
796  pot => potentials(iradtyp)%pot
797  dx = pot%dx
798  pot0_2 => pot%pot0_2
799  !$OMP PARALLEL DO DEFAULT(NONE) &
800  !$OMP SHARED(pot, par_scheme, para_env, mm_atom_index, mm_particles, dOmmOqm, mm_cell, qmmm_spherical_cutoff) &
801  !$OMP SHARED(bo, gbo, dr1, dr2, dr3, grid2, shells, pot0_2, dx, mm_charges, IRadTyp) &
802  !$OMP PRIVATE(myind, Imm, LIndMM, IndMM, ra, qt, sph_chrg_factor, rt1, rt2, rt3, my_k, my_j) &
803  !$OMP PRIVATE(rv1, rv2, rv3, rx2, rx3, r, r2, rx, Term, xs1, xs2, xs3, i, j, k, ix)
804  atoms: DO imm = 1, SIZE(pot%mm_atom_index)
805  IF (par_scheme == do_par_atom) THEN
806  myind = imm + (iradtyp - 1)*SIZE(pot%mm_atom_index)
807  IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle
808  END IF
809  lindmm = pot%mm_atom_index(imm)
810  indmm = mm_atom_index(lindmm)
811  ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
812  qt = mm_charges(lindmm)
813  IF (shells) &
814  ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
815  ! Possible Spherical Cutoff
816  IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
817  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
818  qt = qt*sph_chrg_factor
819  END IF
820  IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
821  rt1 = ra(1)
822  rt2 = ra(2)
823  rt3 = ra(3)
824  loopongrid: DO k = bo(1, 3), bo(2, 3)
825  my_k = k - gbo(1, 3)
826  xs3 = real(my_k, dp)*dr3
827  my_j = bo(1, 2) - gbo(1, 2)
828  xs2 = real(my_j, dp)*dr2
829  rv3 = rt3 - xs3
830  DO j = bo(1, 2), bo(2, 2)
831  xs1 = (bo(1, 1) - gbo(1, 1))*dr1
832  rv2 = rt2 - xs2
833  DO i = bo(1, 1), bo(2, 1)
834  rv1 = rt1 - xs1
835  r2 = rv1*rv1 + rv2*rv2 + rv3*rv3
836  r = sqrt(r2)
837  ix = floor(r/dx) + 1
838  rx = (r - real(ix - 1, dp)*dx)/dx
839  rx2 = rx*rx
840  rx3 = rx2*rx
841  term = pot0_2(1, ix)*(1.0_dp - 3.0_dp*rx2 + 2.0_dp*rx3) &
842  + pot0_2(2, ix)*(rx - 2.0_dp*rx2 + rx3) &
843  + pot0_2(1, ix + 1)*(3.0_dp*rx2 - 2.0_dp*rx3) &
844  + pot0_2(2, ix + 1)*(-rx2 + rx3)
845  !$OMP ATOMIC
846  grid2(i, j, k) = grid2(i, j, k) - term*qt
847  !$OMP END ATOMIC
848  xs1 = xs1 + dr1
849  END DO
850  xs2 = xs2 + dr2
851  END DO
852  END DO loopongrid
853  END DO atoms
854  !$OMP END PARALLEL DO
855  END DO radius
856  CALL timestop(handle)
857  END SUBROUTINE qmmm_elec_with_gaussian_lr
858 
859 END MODULE qmmm_gpw_energy
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
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_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition: cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Definition: cp_realspace_grid_cube.F:45
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::pw_prolongate_s3
subroutine, public pw_prolongate_s3(pw_coarse_in, pw_fine_out, coarse_pool, param_section)
prolongates a function from a coarse grid into a fine one
Definition: cp_spline_utils.F:138
cp_spline_utils::spline3_pbc_interp
integer, parameter, public spline3_pbc_interp
Definition: cp_spline_utils.F:38
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_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition: input_section_types.F:989
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
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
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_el_coupling
subroutine, public qmmm_el_coupling(qs_env, qmmm_env, mm_particles, mm_cell)
Main Driver to compute the QM/MM Electrostatic Coupling.
Definition: qmmm_gpw_energy.F:93
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_se_energy
Calculation of the QMMM Hamiltonian integral matrix <a|\sum_i q_i|b> for semi-empirical methods.
Definition: qmmm_se_energy.F:13
qmmm_se_energy::build_se_qmmm_matrix
subroutine, public build_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
Constructs the 1-el semi-empirical hamiltonian.
Definition: qmmm_se_energy.F:88
qmmm_tb_methods
TB methods used with QMMM.
Definition: qmmm_tb_methods.F:12
qmmm_tb_methods::build_tb_qmmm_matrix_zero
subroutine, public build_tb_qmmm_matrix_zero(qs_env, para_env)
Constructs an empty 1-el DFTB hamiltonian.
Definition: qmmm_tb_methods.F:244
qmmm_tb_methods::build_tb_qmmm_matrix_pc
subroutine, public build_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
Constructs the 1-el DFTB hamiltonian.
Definition: qmmm_tb_methods.F:310
qmmm_tb_methods::build_tb_qmmm_matrix
subroutine, public build_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env)
Constructs the 1-el DFTB hamiltonian.
Definition: qmmm_tb_methods.F:102
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_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_subsys_types
types that represent a quickstep subsys
Definition: qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition: qs_subsys_types.F:134