(git:ed6f26b)
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-2025 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 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
859END 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:120
cell_types::pbc
Definition cell_types.F:92
pw_methods::pw_zero
Definition pw_methods.F:82
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:1063
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:287
input_constants::do_qmmm_none
integer, parameter, public do_qmmm_none
Definition input_constants.F:385
input_constants::do_qmmm_pcharge
integer, parameter, public do_qmmm_pcharge
Definition input_constants.F:385
input_constants::do_qmmm_coulomb
integer, parameter, public do_qmmm_coulomb
Definition input_constants.F:385
input_constants::do_qmmm_swave
integer, parameter, public do_qmmm_swave
Definition input_constants.F:385
input_constants::do_qmmm_gauss
integer, parameter, public do_qmmm_gauss
Definition input_constants.F:385
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:153
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:243
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:757
qmmm_gpw_energy::qmmm_elec_with_gaussian_lg
subroutine, public qmmm_elec_with_gaussian_lg(pgfs, cgrid, mm_charges, mm_atom_index, mm_particles, para_env, per_potentials, mm_cell, dommoqm, par_scheme, qmmm_spherical_cutoff, shells)
Compute the QM/MM electrostatic Interaction collocating (1/R - Sum_NG Gaussians) on the coarser grid ...
Definition qmmm_gpw_energy.F:532
qmmm_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_pp, 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, 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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
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:763
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:217
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