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