(git:d4bfb39)
Loading...
Searching...
No Matches
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-2026 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! **************************************************************************************************
14MODULE 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!***
80CONTAINS
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 cpabort("Unknown QM/MM coupling")
150 END SELECT
151 ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
152 ! DFTB
153 SELECT CASE (qmmm_env%qmmm_coupl_type)
154 CASE (do_qmmm_none)
155 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
156 "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
157 CALL build_tb_qmmm_matrix_zero(qs_env, para_env)
158 CASE (do_qmmm_coulomb)
159 CALL build_tb_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
160 CASE (do_qmmm_pcharge)
161 CALL build_tb_qmmm_matrix_pc(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
162 CASE (do_qmmm_gauss, do_qmmm_swave)
163 cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not implemented for DFTB.")
164 CASE DEFAULT
165 cpabort("Unknown QM/MM coupling")
166 END SELECT
167 ELSE
168 ! QS
169 SELECT CASE (qmmm_env%qmmm_coupl_type)
170 CASE (do_qmmm_coulomb)
171 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
172 CASE (do_qmmm_pcharge)
173 cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
174 CASE (do_qmmm_gauss, do_qmmm_swave)
175 IF (iw > 0) &
176 WRITE (iw, '(T2,"QMMM|",1X,A)') &
177 "QM/MM Coupling computed collocating the Gaussian Potential Functions."
178 interp_section => section_vals_get_subs_vals(input_section, &
179 "QMMM%INTERPOLATOR")
180 CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env, &
181 v_qmmm=ks_qmmm_env_loc%v_qmmm_rspace, &
182 mm_particles=mm_particles, &
183 aug_pools=qmmm_env%aug_pools, &
184 para_env=para_env, &
185 eps_mm_rspace=qmmm_env%eps_mm_rspace, &
186 cube_info=ks_qmmm_env_loc%cube_info, &
187 pw_pools=pw_pools, &
188 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
189 coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
190 interp_section=interp_section, &
191 mm_cell=mm_cell)
192 CASE (do_qmmm_none)
193 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
194 "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
195 CASE DEFAULT
196 cpabort("Unknown QM/MM coupling")
197 END SELECT
198 ! Dump info on the electrostatic potential if requested
199 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, &
200 "POTENTIAL"), cp_p_file)) THEN
201 mpi_io = .true.
202 iw2 = cp_print_key_unit_nr(logger, print_section, "POTENTIAL", &
203 extension=".qmmmLog", mpi_io=mpi_io)
204 CALL cp_pw_to_cube(ks_qmmm_env_loc%v_qmmm_rspace, iw2, &
205 particles=particles, &
206 stride=section_get_ivals(print_section, "POTENTIAL%STRIDE"), &
207 title="QM/MM: MM ELECTROSTATIC POTENTIAL ", &
208 mpi_io=mpi_io)
209 CALL cp_print_key_finished_output(iw2, logger, print_section, &
210 "POTENTIAL", mpi_io=mpi_io)
211 END IF
212 END IF
213 CALL cp_print_key_finished_output(iw, logger, print_section, &
214 "PROGRAM_RUN_INFO")
215 CALL timestop(handle)
216 END SUBROUTINE qmmm_el_coupling
217
218! **************************************************************************************************
219!> \brief Compute the QM/MM electrostatic Interaction collocating the gaussian
220!> Electrostatic Potential
221!> \param qmmm_env ...
222!> \param v_qmmm ...
223!> \param mm_particles ...
224!> \param aug_pools ...
225!> \param cube_info ...
226!> \param para_env ...
227!> \param eps_mm_rspace ...
228!> \param pw_pools ...
229!> \param auxbas_grid ...
230!> \param coarser_grid ...
231!> \param interp_section ...
232!> \param mm_cell ...
233!> \par History
234!> 06.2004 created [tlaino]
235!> \author Teodoro Laino
236! **************************************************************************************************
237 SUBROUTINE qmmm_elec_with_gaussian(qmmm_env, v_qmmm, mm_particles, &
238 aug_pools, cube_info, para_env, eps_mm_rspace, pw_pools, &
239 auxbas_grid, coarser_grid, interp_section, mm_cell)
240 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
241 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_qmmm
242 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
243 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
244 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
245 TYPE(mp_para_env_type), POINTER :: para_env
246 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
247 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
248 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
249 TYPE(section_vals_type), POINTER :: interp_section
250 TYPE(cell_type), POINTER :: mm_cell
251
252 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian'
253
254 INTEGER :: handle, handle2, igrid, ilevel, &
255 kind_interp, lb(3), ngrids, ub(3)
256 LOGICAL :: shells
257 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
258
259 cpassert(ASSOCIATED(mm_particles))
260 cpassert(ASSOCIATED(qmmm_env%mm_atom_chrg))
261 cpassert(ASSOCIATED(qmmm_env%mm_atom_index))
262 cpassert(ASSOCIATED(aug_pools))
263 cpassert(ASSOCIATED(pw_pools))
264 !Statements
265 CALL timeset(routinen, handle)
266 ngrids = SIZE(pw_pools)
267 CALL pw_pools_create_pws(aug_pools, grids)
268 DO igrid = 1, ngrids
269 CALL pw_zero(grids(igrid))
270 END DO
271
272 shells = .false.
273
274 CALL qmmm_elec_with_gaussian_low(grids, mm_particles, &
275 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
276 cube_info, para_env, eps_mm_rspace, qmmm_env%pgfs, &
277 auxbas_grid, coarser_grid, qmmm_env%potentials, &
278 mm_cell=mm_cell, dommoqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
279 per_potentials=qmmm_env%per_potentials, par_scheme=qmmm_env%par_scheme, &
280 qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)
281
282 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
283 CALL qmmm_elec_with_gaussian_low(grids, qmmm_env%added_charges%added_particles, &
284 qmmm_env%added_charges%mm_atom_chrg, &
285 qmmm_env%added_charges%mm_atom_index, &
286 cube_info, para_env, eps_mm_rspace, qmmm_env%added_charges%pgfs, auxbas_grid, &
287 coarser_grid, qmmm_env%added_charges%potentials, &
288 mm_cell=mm_cell, dommoqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
289 per_potentials=qmmm_env%added_charges%per_potentials, par_scheme=qmmm_env%par_scheme, &
290 qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)
291 END IF
292 IF (qmmm_env%added_shells%num_mm_atoms > 0) THEN
293 shells = .true.
294 CALL qmmm_elec_with_gaussian_low(grids, qmmm_env%added_shells%added_particles, &
295 qmmm_env%added_shells%mm_core_chrg, &
296 qmmm_env%added_shells%mm_core_index, &
297 cube_info, para_env, eps_mm_rspace, qmmm_env%added_shells%pgfs, auxbas_grid, &
298 coarser_grid, qmmm_env%added_shells%potentials, &
299 mm_cell=mm_cell, dommoqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
300 per_potentials=qmmm_env%added_shells%per_potentials, &
301 par_scheme=qmmm_env%par_scheme, qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, &
302 shells=shells)
303 END IF
304 ! Sumup all contributions according the parallelization scheme
305 IF (qmmm_env%par_scheme == do_par_atom) THEN
306 DO ilevel = 1, SIZE(grids)
307 CALL para_env%sum(grids(ilevel)%array)
308 END DO
309 END IF
310 ! RealSpace Interpolation
311 CALL section_vals_val_get(interp_section, "kind", i_val=kind_interp)
312 SELECT CASE (kind_interp)
313 CASE (spline3_nopbc_interp, spline3_pbc_interp)
314 ! Spline Iterpolator
315 CALL para_env%sync()
316 CALL timeset(trim(routinen)//":spline3Int", handle2)
317 DO ilevel = coarser_grid, auxbas_grid + 1, -1
318 CALL pw_prolongate_s3(grids(ilevel), &
319 grids(ilevel - 1), &
320 aug_pools(ilevel)%pool, &
321 param_section=interp_section)
322 END DO
323 CALL timestop(handle2)
324 CASE DEFAULT
325 cpabort("Unknown kind interpolation")
326 END SELECT
327 lb = v_qmmm%pw_grid%bounds_local(1, :)
328 ub = v_qmmm%pw_grid%bounds_local(2, :)
329
330 v_qmmm%array = grids(auxbas_grid)%array(lb(1):ub(1), &
331 lb(2):ub(2), &
332 lb(3):ub(3))
333
334 CALL pw_pools_give_back_pws(aug_pools, grids)
335
336 CALL timestop(handle)
337 END SUBROUTINE qmmm_elec_with_gaussian
338
339! **************************************************************************************************
340!> \brief Compute the QM/MM electrostatic Interaction collocating the gaussian
341!> Electrostatic Potential - Low Level
342!> \param tmp_grid ...
343!> \param mm_particles ...
344!> \param mm_charges ...
345!> \param mm_atom_index ...
346!> \param cube_info ...
347!> \param para_env ...
348!> \param eps_mm_rspace ...
349!> \param pgfs ...
350!> \param auxbas_grid ...
351!> \param coarser_grid ...
352!> \param potentials ...
353!> \param mm_cell ...
354!> \param dOmmOqm ...
355!> \param periodic ...
356!> \param per_potentials ...
357!> \param par_scheme ...
358!> \param qmmm_spherical_cutoff ...
359!> \param shells ...
360!> \par History
361!> 06.2004 created [tlaino]
362!> \author Teodoro Laino
363! **************************************************************************************************
364 SUBROUTINE qmmm_elec_with_gaussian_low(tmp_grid, mm_particles, mm_charges, &
365 mm_atom_index, cube_info, para_env, &
366 eps_mm_rspace, pgfs, auxbas_grid, coarser_grid, &
367 potentials, mm_cell, dOmmOqm, periodic, per_potentials, par_scheme, &
368 qmmm_spherical_cutoff, shells)
369 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: tmp_grid
370 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
371 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
372 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
373 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
374 TYPE(mp_para_env_type), POINTER :: para_env
375 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
376 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
377 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
378 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
379 TYPE(cell_type), POINTER :: mm_cell
380 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
381 LOGICAL, INTENT(IN) :: periodic
382 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
383 INTEGER, INTENT(IN) :: par_scheme
384 REAL(kind=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
385 LOGICAL, INTENT(IN) :: shells
386
387 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian_low', &
388 routinenb = 'qmmm_elec_gaussian_low'
389
390 INTEGER :: handle, handle2, igauss, ilevel, imm, &
391 indmm, iradtyp, lindmm, myind, &
392 n_rep_real(3)
393 INTEGER, DIMENSION(2, 3) :: bo2
394 REAL(kind=dp) :: alpha, height, sph_chrg_factor, w
395 REAL(kind=dp), DIMENSION(3) :: ra
396 REAL(kind=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat
397 TYPE(qmmm_gaussian_type), POINTER :: pgf
398 TYPE(qmmm_per_pot_type), POINTER :: per_pot
399 TYPE(qmmm_pot_type), POINTER :: pot
400
401 NULLIFY (pgf, pot, per_pot, xdat, ydat, zdat)
402 CALL timeset(routinen, handle)
403 CALL timeset(routinenb//"_G", handle2)
404 bo2 = tmp_grid(auxbas_grid)%pw_grid%bounds
405 ALLOCATE (xdat(bo2(1, 1):bo2(2, 1)))
406 ALLOCATE (ydat(bo2(1, 2):bo2(2, 2)))
407 ALLOCATE (zdat(bo2(1, 3):bo2(2, 3)))
408 IF (par_scheme == do_par_atom) myind = 0
409 radius: DO iradtyp = 1, SIZE(pgfs)
410 pgf => pgfs(iradtyp)%pgf
411 pot => potentials(iradtyp)%pot
412 n_rep_real = 0
413 IF (periodic) THEN
414 per_pot => per_potentials(iradtyp)%pot
415 n_rep_real = per_pot%n_rep_real
416 END IF
417 gaussian: DO igauss = 1, pgf%Number_of_Gaussians
418 alpha = 1.0_dp/pgf%Gk(igauss)
419 alpha = alpha*alpha
420 height = pgf%Ak(igauss)
421 ilevel = pgf%grid_level(igauss)
422 atoms: DO imm = 1, SIZE(pot%mm_atom_index)
423 IF (par_scheme == do_par_atom) THEN
424 myind = myind + 1
425 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
426 END IF
427 lindmm = pot%mm_atom_index(imm)
428 indmm = mm_atom_index(lindmm)
429 IF (shells) THEN
430 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
431 ELSE
432 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
433 END IF
434 w = mm_charges(lindmm)*height
435 ! Possible Spherical Cutoff
436 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
437 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
438 w = w*sph_chrg_factor
439 END IF
440 IF (abs(w) <= epsilon(0.0_dp)) cycle atoms
441 CALL collocate_gf_rspace_nopbc(zetp=alpha, &
442 rp=ra, &
443 scale=-1.0_dp, &
444 w=w, &
445 pwgrid=tmp_grid(ilevel), &
446 cube_info=cube_info(ilevel), &
447 eps_mm_rspace=eps_mm_rspace, &
448 xdat=xdat, &
449 ydat=ydat, &
450 zdat=zdat, &
451 bo2=bo2, &
452 n_rep_real=n_rep_real, &
453 mm_cell=mm_cell)
454 END DO atoms
455 END DO gaussian
456 END DO radius
457 IF (ASSOCIATED(xdat)) THEN
458 DEALLOCATE (xdat)
459 END IF
460 IF (ASSOCIATED(ydat)) THEN
461 DEALLOCATE (ydat)
462 END IF
463 IF (ASSOCIATED(zdat)) THEN
464 DEALLOCATE (zdat)
465 END IF
466 CALL timestop(handle2)
467 CALL timeset(routinenb//"_R", handle2)
468 IF (periodic) THEN
469 ! Long Range Part of the QM/MM Potential with Gaussians With Periodic Boundary Conditions
470 CALL qmmm_elec_with_gaussian_lg(pgfs=pgfs, &
471 cgrid=tmp_grid(coarser_grid), &
472 mm_charges=mm_charges, &
473 mm_atom_index=mm_atom_index, &
474 mm_particles=mm_particles, &
475 para_env=para_env, &
476 per_potentials=per_potentials, &
477 mm_cell=mm_cell, &
478 dommoqm=dommoqm, &
479 par_scheme=par_scheme, &
480 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
481 shells=shells)
482 ELSE
483 ! Long Range Part of the QM/MM Potential with Gaussians
484 CALL qmmm_elec_with_gaussian_lr(pgfs=pgfs, &
485 grid=tmp_grid(coarser_grid), &
486 mm_charges=mm_charges, &
487 mm_atom_index=mm_atom_index, &
488 mm_particles=mm_particles, &
489 para_env=para_env, &
490 potentials=potentials, &
491 mm_cell=mm_cell, &
492 dommoqm=dommoqm, &
493 par_scheme=par_scheme, &
494 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
495 shells=shells)
496 END IF
497 CALL timestop(handle2)
498 CALL timestop(handle)
499
500 END SUBROUTINE qmmm_elec_with_gaussian_low
501
502! **************************************************************************************************
503!> \brief Compute the QM/MM electrostatic Interaction collocating
504!> (1/R - Sum_NG Gaussians) on the coarser grid level in G-SPACE
505!> Long Range QM/MM Electrostatic Potential with Gaussian - Low Level
506!> PERIODIC BOUNDARY CONDITION VERSION
507!> \param pgfs ...
508!> \param cgrid ...
509!> \param mm_charges ...
510!> \param mm_atom_index ...
511!> \param mm_particles ...
512!> \param para_env ...
513!> \param per_potentials ...
514!> \param mm_cell ...
515!> \param dOmmOqm ...
516!> \param par_scheme ...
517!> \param qmmm_spherical_cutoff ...
518!> \param shells ...
519!> \par History
520!> 07.2004 created [tlaino]
521!> \author Teodoro Laino
522!> \note
523!> This version includes the explicit code of Eval_Interp_Spl3_pbc
524!> in order to achieve better performance
525! **************************************************************************************************
526 SUBROUTINE qmmm_elec_with_gaussian_lg(pgfs, cgrid, mm_charges, mm_atom_index, &
527 mm_particles, para_env, per_potentials, &
528 mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
529 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
530 TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
531 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
532 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
533 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
534 TYPE(mp_para_env_type), POINTER :: para_env
535 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
536 TYPE(cell_type), POINTER :: mm_cell
537 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
538 INTEGER, INTENT(IN) :: par_scheme
539 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
540 LOGICAL :: shells
541
542 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian_LG'
543
544 INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, &
545 ij3, ij4, ik1, ik2, ik3, ik4, imm, &
546 indmm, iradtyp, ivec(3), j, k, lindmm, &
547 my_j, my_k, myind, npts(3)
548 INTEGER, DIMENSION(2, 3) :: bo, gbo
549 REAL(kind=dp) :: a1, a2, a3, abc_x(4, 4), abc_x_y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, &
550 dr1, dr1c, dr2, dr2c, dr3, dr3c, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, &
551 p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, rt3, rv1, rv2, rv3, s1, s2, s3, s4, &
552 sph_chrg_factor, t1, t2, t3, t4, u1, u2, u3, v1, v2, v3, v4, val, xd1, xd2, xd3, xs1, &
553 xs2, xs3
554 REAL(kind=dp), DIMENSION(3) :: ra, vec
555 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid, grid2
556 TYPE(pw_r3d_rs_type), POINTER :: pw
557 TYPE(qmmm_per_pot_type), POINTER :: per_pot
558
559 CALL timeset(routinen, handle)
560 NULLIFY (grid, pw)
561 dr1c = cgrid%pw_grid%dr(1)
562 dr2c = cgrid%pw_grid%dr(2)
563 dr3c = cgrid%pw_grid%dr(3)
564 gbo = cgrid%pw_grid%bounds
565 bo = cgrid%pw_grid%bounds_local
566 grid2 => cgrid%array
567 IF (par_scheme == do_par_atom) myind = 0
568 radius: DO iradtyp = 1, SIZE(pgfs)
569 per_pot => per_potentials(iradtyp)%pot
570 pw => per_pot%TabLR
571 npts = pw%pw_grid%npts
572 dr1 = pw%pw_grid%dr(1)
573 dr2 = pw%pw_grid%dr(2)
574 dr3 = pw%pw_grid%dr(3)
575 grid => pw%array(:, :, :)
576 !$OMP PARALLEL DO DEFAULT(NONE) &
577 !$OMP SHARED(bo, gbo, grid, grid2, pw, npts, per_pot, mm_atom_index) &
578 !$OMP SHARED(dr1, dr2, dr3, dr1c, dr2c, dr3c, par_scheme, mm_charges, mm_particles) &
579 !$OMP SHARED(mm_cell, dOmmOqm, shells, para_env, IRadTyp, qmmm_spherical_cutoff) &
580 !$OMP PRIVATE(Imm, LIndMM, IndMM, qt, sph_chrg_factor, ra, myind) &
581 !$OMP PRIVATE(rt1, rt2, rt3, k, vec, ivec, xd1, xd2, xd3, ik1, ik2, ik3, ik4) &
582 !$OMP PRIVATE(ij1, ij2, ij3, ij4, ii1, ii2, ii3, ii4, my_k, my_j, xs1, xs2, xs3) &
583 !$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, v1, v2, v3, v4, e1, e2, e3) &
584 !$OMP PRIVATE(f1, f2, f3, g1, g2, g3, h1, h2, h3, s1, s2, s3, s4, a1, a2, a3) &
585 !$OMP PRIVATE(b1, b2, b3, c1, c2, c3, d1, d2, d3, t1, t2, t3, t4, u1, u2, u3, val) &
586 !$OMP PRIVATE(rv1, rv2, rv3, abc_X, abc_X_Y)
587 atoms: DO imm = 1, SIZE(per_pot%mm_atom_index)
588 IF (par_scheme == do_par_atom) THEN
589 myind = imm + (iradtyp - 1)*SIZE(per_pot%mm_atom_index)
590 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
591 END IF
592 lindmm = per_pot%mm_atom_index(imm)
593 indmm = mm_atom_index(lindmm)
594 qt = mm_charges(lindmm)
595 IF (shells) THEN
596 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
597 ELSE
598 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
599 END IF
600 ! Possible Spherical Cutoff
601 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
602 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
603 qt = qt*sph_chrg_factor
604 END IF
605 IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
606 rt1 = ra(1)
607 rt2 = ra(2)
608 rt3 = ra(3)
609 loopongrid: DO k = bo(1, 3), bo(2, 3)
610 my_k = k - gbo(1, 3)
611 xs3 = real(my_k, dp)*dr3c
612 my_j = bo(1, 2) - gbo(1, 2)
613 xs2 = real(my_j, dp)*dr2c
614 rv3 = rt3 - xs3
615 vec(3) = rv3
616 ivec(3) = floor(vec(3)/pw%pw_grid%dr(3))
617 xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
618 ik1 = modulo(ivec(3) - 1, npts(3)) + 1
619 ik2 = modulo(ivec(3), npts(3)) + 1
620 ik3 = modulo(ivec(3) + 1, npts(3)) + 1
621 ik4 = modulo(ivec(3) + 2, npts(3)) + 1
622 p1 = 3.0_dp + xd3
623 p2 = p1*p1
624 p3 = p2*p1
625 q1 = 2.0_dp + xd3
626 q2 = q1*q1
627 q3 = q2*q1
628 r1 = 1.0_dp + xd3
629 r2 = r1*r1
630 r3 = r2*r1
631 u1 = xd3
632 u2 = u1*u1
633 u3 = u2*u1
634 v1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
635 v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
636 v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
637 v4 = 1.0_dp/6.0_dp*u3
638 DO j = bo(1, 2), bo(2, 2)
639 xs1 = (bo(1, 1) - gbo(1, 1))*dr1c
640 rv2 = rt2 - xs2
641 vec(2) = rv2
642 ivec(2) = floor(vec(2)/pw%pw_grid%dr(2))
643 xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
644 ij1 = modulo(ivec(2) - 1, npts(2)) + 1
645 ij2 = modulo(ivec(2), npts(2)) + 1
646 ij3 = modulo(ivec(2) + 1, npts(2)) + 1
647 ij4 = modulo(ivec(2) + 2, npts(2)) + 1
648 e1 = 3.0_dp + xd2
649 e2 = e1*e1
650 e3 = e2*e1
651 f1 = 2.0_dp + xd2
652 f2 = f1*f1
653 f3 = f2*f1
654 g1 = 1.0_dp + xd2
655 g2 = g1*g1
656 g3 = g2*g1
657 h1 = xd2
658 h2 = h1*h1
659 h3 = h2*h1
660 s1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
661 s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
662 s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
663 s4 = 1.0_dp/6.0_dp*h3
664 DO i = bo(1, 1), bo(2, 1)
665 rv1 = rt1 - xs1
666 vec(1) = rv1
667 ivec(1) = floor(vec(1)/pw%pw_grid%dr(1))
668 xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
669 ii1 = modulo(ivec(1) - 1, npts(1)) + 1
670 ii2 = modulo(ivec(1), npts(1)) + 1
671 ii3 = modulo(ivec(1) + 1, npts(1)) + 1
672 ii4 = modulo(ivec(1) + 2, npts(1)) + 1
673 !
674 ! Spline Interpolation
675 !
676
677 a1 = 3.0_dp + xd1
678 a2 = a1*a1
679 a3 = a2*a1
680 b1 = 2.0_dp + xd1
681 b2 = b1*b1
682 b3 = b2*b1
683 c1 = 1.0_dp + xd1
684 c2 = c1*c1
685 c3 = c2*c1
686 d1 = xd1
687 d2 = d1*d1
688 d3 = d2*d1
689 t1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
690 t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
691 t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
692 t4 = 1.0_dp/6.0_dp*d3
693
694 abc_x(1, 1) = grid(ii1, ij1, ik1)*v1 + grid(ii1, ij1, ik2)*v2 + grid(ii1, ij1, ik3)*v3 + grid(ii1, ij1, ik4)*v4
695 abc_x(1, 2) = grid(ii1, ij2, ik1)*v1 + grid(ii1, ij2, ik2)*v2 + grid(ii1, ij2, ik3)*v3 + grid(ii1, ij2, ik4)*v4
696 abc_x(1, 3) = grid(ii1, ij3, ik1)*v1 + grid(ii1, ij3, ik2)*v2 + grid(ii1, ij3, ik3)*v3 + grid(ii1, ij3, ik4)*v4
697 abc_x(1, 4) = grid(ii1, ij4, ik1)*v1 + grid(ii1, ij4, ik2)*v2 + grid(ii1, ij4, ik3)*v3 + grid(ii1, ij4, ik4)*v4
698 abc_x(2, 1) = grid(ii2, ij1, ik1)*v1 + grid(ii2, ij1, ik2)*v2 + grid(ii2, ij1, ik3)*v3 + grid(ii2, ij1, ik4)*v4
699 abc_x(2, 2) = grid(ii2, ij2, ik1)*v1 + grid(ii2, ij2, ik2)*v2 + grid(ii2, ij2, ik3)*v3 + grid(ii2, ij2, ik4)*v4
700 abc_x(2, 3) = grid(ii2, ij3, ik1)*v1 + grid(ii2, ij3, ik2)*v2 + grid(ii2, ij3, ik3)*v3 + grid(ii2, ij3, ik4)*v4
701 abc_x(2, 4) = grid(ii2, ij4, ik1)*v1 + grid(ii2, ij4, ik2)*v2 + grid(ii2, ij4, ik3)*v3 + grid(ii2, ij4, ik4)*v4
702 abc_x(3, 1) = grid(ii3, ij1, ik1)*v1 + grid(ii3, ij1, ik2)*v2 + grid(ii3, ij1, ik3)*v3 + grid(ii3, ij1, ik4)*v4
703 abc_x(3, 2) = grid(ii3, ij2, ik1)*v1 + grid(ii3, ij2, ik2)*v2 + grid(ii3, ij2, ik3)*v3 + grid(ii3, ij2, ik4)*v4
704 abc_x(3, 3) = grid(ii3, ij3, ik1)*v1 + grid(ii3, ij3, ik2)*v2 + grid(ii3, ij3, ik3)*v3 + grid(ii3, ij3, ik4)*v4
705 abc_x(3, 4) = grid(ii3, ij4, ik1)*v1 + grid(ii3, ij4, ik2)*v2 + grid(ii3, ij4, ik3)*v3 + grid(ii3, ij4, ik4)*v4
706 abc_x(4, 1) = grid(ii4, ij1, ik1)*v1 + grid(ii4, ij1, ik2)*v2 + grid(ii4, ij1, ik3)*v3 + grid(ii4, ij1, ik4)*v4
707 abc_x(4, 2) = grid(ii4, ij2, ik1)*v1 + grid(ii4, ij2, ik2)*v2 + grid(ii4, ij2, ik3)*v3 + grid(ii4, ij2, ik4)*v4
708 abc_x(4, 3) = grid(ii4, ij3, ik1)*v1 + grid(ii4, ij3, ik2)*v2 + grid(ii4, ij3, ik3)*v3 + grid(ii4, ij3, ik4)*v4
709 abc_x(4, 4) = grid(ii4, ij4, ik1)*v1 + grid(ii4, ij4, ik2)*v2 + grid(ii4, ij4, ik3)*v3 + grid(ii4, ij4, ik4)*v4
710
711 abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
712 abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
713 abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
714 abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
715
716 val = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
717 !$OMP ATOMIC
718 grid2(i, j, k) = grid2(i, j, k) - val*qt
719 !$OMP END ATOMIC
720 xs1 = xs1 + dr1c
721 END DO
722 xs2 = xs2 + dr2c
723 END DO
724 END DO loopongrid
725 END DO atoms
726 !$OMP END PARALLEL DO
727 END DO radius
728 CALL timestop(handle)
729 END SUBROUTINE qmmm_elec_with_gaussian_lg
730
731! **************************************************************************************************
732!> \brief Compute the QM/MM electrostatic Interaction collocating
733!> (1/R - Sum_NG Gaussians) on the coarser grid level.
734!> Long Range QM/MM Electrostatic Potential with Gaussian - Low Level
735!> \param pgfs ...
736!> \param grid ...
737!> \param mm_charges ...
738!> \param mm_atom_index ...
739!> \param mm_particles ...
740!> \param para_env ...
741!> \param potentials ...
742!> \param mm_cell ...
743!> \param dOmmOqm ...
744!> \param par_scheme ...
745!> \param qmmm_spherical_cutoff ...
746!> \param shells ...
747!> \par History
748!> 07.2004 created [tlaino]
749!> \author Teodoro Laino
750! **************************************************************************************************
751 SUBROUTINE qmmm_elec_with_gaussian_lr(pgfs, grid, mm_charges, mm_atom_index, &
752 mm_particles, para_env, potentials, &
753 mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
754 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
755 TYPE(pw_r3d_rs_type), INTENT(IN) :: grid
756 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
757 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
758 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
759 TYPE(mp_para_env_type), POINTER :: para_env
760 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
761 TYPE(cell_type), POINTER :: mm_cell
762 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
763 INTEGER, INTENT(IN) :: par_scheme
764 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
765 LOGICAL :: shells
766
767 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_elec_with_gaussian_LR'
768
769 INTEGER :: handle, i, imm, indmm, iradtyp, ix, j, &
770 k, lindmm, my_j, my_k, myind, n1, n2, &
771 n3
772 INTEGER, DIMENSION(2, 3) :: bo, gbo
773 REAL(kind=dp) :: dr1, dr2, dr3, dx, qt, r, r2, rt1, rt2, &
774 rt3, rv1, rv2, rv3, rx, rx2, rx3, &
775 sph_chrg_factor, term, xs1, xs2, xs3
776 REAL(kind=dp), DIMENSION(3) :: ra
777 REAL(kind=dp), DIMENSION(:, :), POINTER :: pot0_2
778 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid2
779 TYPE(qmmm_pot_type), POINTER :: pot
780
781 CALL timeset(routinen, handle)
782 n1 = grid%pw_grid%npts(1)
783 n2 = grid%pw_grid%npts(2)
784 n3 = grid%pw_grid%npts(3)
785 dr1 = grid%pw_grid%dr(1)
786 dr2 = grid%pw_grid%dr(2)
787 dr3 = grid%pw_grid%dr(3)
788 gbo = grid%pw_grid%bounds
789 bo = grid%pw_grid%bounds_local
790 grid2 => grid%array
791 IF (par_scheme == do_par_atom) myind = 0
792 radius: DO iradtyp = 1, SIZE(pgfs)
793 pot => potentials(iradtyp)%pot
794 dx = pot%dx
795 pot0_2 => pot%pot0_2
796 !$OMP PARALLEL DO DEFAULT(NONE) &
797 !$OMP SHARED(pot, par_scheme, para_env, mm_atom_index, mm_particles, dOmmOqm, mm_cell, qmmm_spherical_cutoff) &
798 !$OMP SHARED(bo, gbo, dr1, dr2, dr3, grid2, shells, pot0_2, dx, mm_charges, IRadTyp) &
799 !$OMP PRIVATE(myind, Imm, LIndMM, IndMM, ra, qt, sph_chrg_factor, rt1, rt2, rt3, my_k, my_j) &
800 !$OMP PRIVATE(rv1, rv2, rv3, rx2, rx3, r, r2, rx, Term, xs1, xs2, xs3, i, j, k, ix)
801 atoms: DO imm = 1, SIZE(pot%mm_atom_index)
802 IF (par_scheme == do_par_atom) THEN
803 myind = imm + (iradtyp - 1)*SIZE(pot%mm_atom_index)
804 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
805 END IF
806 lindmm = pot%mm_atom_index(imm)
807 indmm = mm_atom_index(lindmm)
808 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
809 qt = mm_charges(lindmm)
810 IF (shells) &
811 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
812 ! Possible Spherical Cutoff
813 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
814 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
815 qt = qt*sph_chrg_factor
816 END IF
817 IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
818 rt1 = ra(1)
819 rt2 = ra(2)
820 rt3 = ra(3)
821 loopongrid: DO k = bo(1, 3), bo(2, 3)
822 my_k = k - gbo(1, 3)
823 xs3 = real(my_k, dp)*dr3
824 my_j = bo(1, 2) - gbo(1, 2)
825 xs2 = real(my_j, dp)*dr2
826 rv3 = rt3 - xs3
827 DO j = bo(1, 2), bo(2, 2)
828 xs1 = (bo(1, 1) - gbo(1, 1))*dr1
829 rv2 = rt2 - xs2
830 DO i = bo(1, 1), bo(2, 1)
831 rv1 = rt1 - xs1
832 r2 = rv1*rv1 + rv2*rv2 + rv3*rv3
833 r = sqrt(r2)
834 ix = floor(r/dx) + 1
835 rx = (r - real(ix - 1, dp)*dx)/dx
836 rx2 = rx*rx
837 rx3 = rx2*rx
838 term = pot0_2(1, ix)*(1.0_dp - 3.0_dp*rx2 + 2.0_dp*rx3) &
839 + pot0_2(2, ix)*(rx - 2.0_dp*rx2 + rx3) &
840 + pot0_2(1, ix + 1)*(3.0_dp*rx2 - 2.0_dp*rx3) &
841 + pot0_2(2, ix + 1)*(-rx2 + rx3)
842 !$OMP ATOMIC
843 grid2(i, j, k) = grid2(i, j, k) - term*qt
844 !$OMP END ATOMIC
845 xs1 = xs1 + dr1
846 END DO
847 xs2 = xs2 + dr2
848 END DO
849 END DO loopongrid
850 END DO atoms
851 !$OMP END PARALLEL DO
852 END DO radius
853 CALL timestop(handle)
854 END SUBROUTINE qmmm_elec_with_gaussian_lr
855
856END MODULE qmmm_gpw_energy
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:125
cell_types::pbc
Definition cell_types.F:103
pw_methods::pw_zero
Definition pw_methods.F:83
pw_pool_types::pw_pools_create_pws
Definition pw_pool_types.F:139
pw_pool_types::pw_pools_give_back_pws
Definition pw_pool_types.F:150
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:853
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:1067
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, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:48
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:18
input_constants::do_par_atom
integer, parameter, public do_par_atom
Definition input_constants.F:296
input_constants::do_qmmm_none
integer, parameter, public do_qmmm_none
Definition input_constants.F:396
input_constants::do_qmmm_pcharge
integer, parameter, public do_qmmm_pcharge
Definition input_constants.F:396
input_constants::do_qmmm_coulomb
integer, parameter, public do_qmmm_coulomb
Definition input_constants.F:396
input_constants::do_qmmm_swave
integer, parameter, public do_qmmm_swave
Definition input_constants.F:396
input_constants::do_qmmm_gauss
integer, parameter, public do_qmmm_gauss
Definition input_constants.F:396
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:160
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:996
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:731
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:1047
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:240
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_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:754
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:529
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:616
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
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
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:751
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cube_utils::cube_info_type
Definition cube_utils.F:38
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qmmm_gaussian_types::qmmm_gaussian_p_type
represent a pointer to a qmmm_gaussian_type, to be able to create arrays of pointers
Definition qmmm_gaussian_types.F:42
qmmm_gaussian_types::qmmm_gaussian_type
Definition qmmm_gaussian_types.F:27
qmmm_types_low::qmmm_env_qm_type
...
Definition qmmm_types_low.F:195
qmmm_types_low::qmmm_per_pot_p_type
Definition qmmm_types_low.F:100
qmmm_types_low::qmmm_per_pot_type
Periodic Potential.
Definition qmmm_types_low.F:90
qmmm_types_low::qmmm_pot_p_type
Definition qmmm_types_low.F:83
qmmm_types_low::qmmm_pot_type
Real Space Potential.
Definition qmmm_types_low.F:76
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_ks_qmmm_types::qs_ks_qmmm_env_type
calculation environment to calculate the ks_qmmm matrix, holds the QM/MM potential and all the needed...
Definition qs_ks_qmmm_types.F:47
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56