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