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_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 : !***
81 : CONTAINS
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 4022 : 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 4022 : 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 4022 : CALL timeset(routineN, handle)
115 4022 : logger => cp_get_default_logger()
116 4022 : 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 4022 : dft_control=dft_control)
124 4022 : CALL qs_subsys_get(subsys, particles=particles)
125 :
126 4022 : CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools)
127 4022 : 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 4022 : extension=".qmmmLog")
130 4022 : IF (iw > 0) THEN
131 1023 : 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 4022 : CALL pw_zero(ks_qmmm_env_loc%v_qmmm_rspace)
137 4022 : IF (dft_control%qs_control%semi_empirical) THEN
138 : ! SEMIEMPIRICAL
139 2892 : SELECT CASE (qmmm_env%qmmm_coupl_type)
140 : CASE (do_qmmm_coulomb, do_qmmm_none)
141 1446 : CALL build_se_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
142 1446 : IF (qmmm_env%qmmm_coupl_type == do_qmmm_none) THEN
143 510 : IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
144 176 : "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
145 : END IF
146 : CASE (do_qmmm_pcharge)
147 0 : CPABORT("Point charge QM/MM electrostatic coupling not yet implemented for SE.")
148 : CASE (do_qmmm_gauss, do_qmmm_swave)
149 0 : CPABORT("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.")
150 : CASE DEFAULT
151 1446 : CPABORT("Unknown QM/MM coupling")
152 : END SELECT
153 2576 : ELSE IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
154 : ! DFTB
155 1800 : SELECT CASE (qmmm_env%qmmm_coupl_type)
156 : CASE (do_qmmm_none)
157 8 : IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
158 4 : "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
159 8 : CALL build_tb_qmmm_matrix_zero(qs_env, para_env)
160 : CASE (do_qmmm_coulomb)
161 448 : CALL build_tb_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
162 : CASE (do_qmmm_pcharge)
163 1116 : CALL build_tb_qmmm_matrix_pc(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
164 : CASE (do_qmmm_gauss)
165 220 : CALL build_tb_qmmm_matrix_gauss(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
166 : CASE (do_qmmm_swave)
167 0 : CPABORT("SWAVE QM/MM electrostatic coupling not implemented for tight-binding methods.")
168 : CASE DEFAULT
169 1792 : CPABORT("Unknown QM/MM coupling")
170 : END SELECT
171 : ELSE
172 : ! QS
173 784 : SELECT CASE (qmmm_env%qmmm_coupl_type)
174 : CASE (do_qmmm_coulomb)
175 0 : CPABORT("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
176 : CASE (do_qmmm_pcharge)
177 0 : CPABORT("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
178 : CASE (do_qmmm_gauss, do_qmmm_swave)
179 744 : IF (iw > 0) THEN
180 : WRITE (iw, '(T2,"QMMM|",1X,A)') &
181 372 : "QM/MM Coupling computed collocating the Gaussian Potential Functions."
182 : END IF
183 : interp_section => section_vals_get_subs_vals(input_section, &
184 744 : "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 744 : mm_cell=mm_cell)
197 : CASE (do_qmmm_none)
198 40 : IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
199 20 : "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
200 : CASE DEFAULT
201 784 : CPABORT("Unknown QM/MM coupling")
202 : END SELECT
203 : ! Dump info on the electrostatic potential if requested
204 784 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, &
205 : "POTENTIAL"), cp_p_file)) THEN
206 24 : mpi_io = .TRUE.
207 : iw2 = cp_print_key_unit_nr(logger, print_section, "POTENTIAL", &
208 24 : 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 24 : mpi_io=mpi_io)
214 : CALL cp_print_key_finished_output(iw2, logger, print_section, &
215 24 : "POTENTIAL", mpi_io=mpi_io)
216 : END IF
217 : END IF
218 : CALL cp_print_key_finished_output(iw, logger, print_section, &
219 4022 : "PROGRAM_RUN_INFO")
220 4022 : CALL timestop(handle)
221 4022 : 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 744 : 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 744 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
263 :
264 744 : CPASSERT(ASSOCIATED(mm_particles))
265 744 : CPASSERT(ASSOCIATED(qmmm_env%mm_atom_chrg))
266 744 : CPASSERT(ASSOCIATED(qmmm_env%mm_atom_index))
267 744 : CPASSERT(ASSOCIATED(aug_pools))
268 744 : CPASSERT(ASSOCIATED(pw_pools))
269 : !Statements
270 744 : CALL timeset(routineN, handle)
271 744 : ngrids = SIZE(pw_pools)
272 744 : CALL pw_pools_create_pws(aug_pools, grids)
273 3748 : DO igrid = 1, ngrids
274 3748 : CALL pw_zero(grids(igrid))
275 : END DO
276 :
277 744 : 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 744 : qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)
286 :
287 744 : 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 32 : qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)
296 : END IF
297 744 : IF (qmmm_env%added_shells%num_mm_atoms > 0) THEN
298 2 : 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 2 : shells=shells)
308 : END IF
309 : ! Sumup all contributions according the parallelization scheme
310 744 : IF (qmmm_env%par_scheme == do_par_atom) THEN
311 3708 : DO ilevel = 1, SIZE(grids)
312 3708 : CALL para_env%sum(grids(ilevel)%array)
313 : END DO
314 : END IF
315 : ! RealSpace Interpolation
316 744 : CALL section_vals_val_get(interp_section, "kind", i_val=kind_interp)
317 744 : SELECT CASE (kind_interp)
318 : CASE (spline3_nopbc_interp, spline3_pbc_interp)
319 : ! Spline Iterpolator
320 744 : CALL para_env%sync()
321 744 : CALL timeset(TRIM(routineN)//":spline3Int", handle2)
322 3004 : 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 3004 : param_section=interp_section)
327 : END DO
328 744 : CALL timestop(handle2)
329 : CASE DEFAULT
330 1488 : CPABORT("Unknown kind interpolation")
331 : END SELECT
332 2976 : lb = v_qmmm%pw_grid%bounds_local(1, :)
333 2976 : 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 33223912 : lb(3):ub(3))
338 :
339 744 : CALL pw_pools_give_back_pws(aug_pools, grids)
340 :
341 744 : CALL timestop(handle)
342 744 : 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 778 : 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 778 : 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 778 : NULLIFY (pgf, pot, per_pot, xdat, ydat, zdat)
407 778 : CALL timeset(routineN, handle)
408 778 : CALL timeset(routineNb//"_G", handle2)
409 7780 : bo2 = tmp_grid(auxbas_grid)%pw_grid%bounds
410 2334 : ALLOCATE (xdat(bo2(1, 1):bo2(2, 1)))
411 2334 : ALLOCATE (ydat(bo2(1, 2):bo2(2, 2)))
412 2334 : ALLOCATE (zdat(bo2(1, 3):bo2(2, 3)))
413 : IF (par_scheme == do_par_atom) myind = 0
414 2176 : Radius: DO IRadTyp = 1, SIZE(pgfs)
415 1398 : pgf => pgfs(IRadTyp)%pgf
416 1398 : pot => potentials(IRadTyp)%pot
417 1398 : n_rep_real = 0
418 1398 : IF (periodic) THEN
419 102 : per_pot => per_potentials(IRadTyp)%pot
420 408 : n_rep_real = per_pot%n_rep_real
421 : END IF
422 12088 : Gaussian: DO IGauss = 1, pgf%Number_of_Gaussians
423 9912 : alpha = 1.0_dp/pgf%Gk(IGauss)
424 9912 : alpha = alpha*alpha
425 9912 : height = pgf%Ak(IGauss)
426 9912 : ilevel = pgf%grid_level(IGauss)
427 60616 : Atoms: DO Imm = 1, SIZE(pot%mm_atom_index)
428 49306 : IF (par_scheme == do_par_atom) THEN
429 48690 : myind = myind + 1
430 48690 : IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE Atoms
431 : END IF
432 25009 : LIndMM = pot%mm_atom_index(Imm)
433 25009 : IndMM = mm_atom_index(LIndMM)
434 25009 : IF (shells) THEN
435 1344 : ra(:) = pbc(mm_particles(LIndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
436 : ELSE
437 198728 : ra(:) = pbc(mm_particles(IndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
438 : END IF
439 25009 : W = mm_charges(LIndMM)*height
440 : ! Possible Spherical Cutoff
441 25009 : IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
442 0 : CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
443 0 : W = W*sph_chrg_factor
444 : END IF
445 25009 : 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 59218 : mm_cell=mm_cell)
459 : END DO Atoms
460 : END DO Gaussian
461 : END DO Radius
462 778 : IF (ASSOCIATED(xdat)) THEN
463 778 : DEALLOCATE (xdat)
464 : END IF
465 778 : IF (ASSOCIATED(ydat)) THEN
466 778 : DEALLOCATE (ydat)
467 : END IF
468 778 : IF (ASSOCIATED(zdat)) THEN
469 778 : DEALLOCATE (zdat)
470 : END IF
471 778 : CALL timestop(handle2)
472 778 : CALL timeset(routineNb//"_R", handle2)
473 778 : 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 64 : 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 714 : shells=shells)
501 : END IF
502 778 : CALL timestop(handle2)
503 778 : CALL timestop(handle)
504 :
505 1556 : 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 64 : 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 64 : 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 64 : CALL timeset(routineN, handle)
565 64 : NULLIFY (grid, pw)
566 64 : dr1c = cgrid%pw_grid%dr(1)
567 64 : dr2c = cgrid%pw_grid%dr(2)
568 64 : dr3c = cgrid%pw_grid%dr(3)
569 640 : gbo = cgrid%pw_grid%bounds
570 640 : bo = cgrid%pw_grid%bounds_local
571 64 : grid2 => cgrid%array
572 64 : IF (par_scheme == do_par_atom) myind = 0
573 166 : Radius: DO IRadTyp = 1, SIZE(pgfs)
574 102 : per_pot => per_potentials(IRadTyp)%pot
575 102 : pw => per_pot%TabLR
576 408 : npts = pw%pw_grid%npts
577 102 : dr1 = pw%pw_grid%dr(1)
578 102 : dr2 = pw%pw_grid%dr(2)
579 102 : dr3 = pw%pw_grid%dr(3)
580 102 : 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 166 : !$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 64 : CALL timestop(handle)
734 64 : 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 714 : 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 714 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pot0_2
783 714 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: grid2
784 : TYPE(qmmm_pot_type), POINTER :: pot
785 :
786 714 : CALL timeset(routineN, handle)
787 714 : n1 = grid%pw_grid%npts(1)
788 714 : n2 = grid%pw_grid%npts(2)
789 714 : n3 = grid%pw_grid%npts(3)
790 714 : dr1 = grid%pw_grid%dr(1)
791 714 : dr2 = grid%pw_grid%dr(2)
792 714 : dr3 = grid%pw_grid%dr(3)
793 7140 : gbo = grid%pw_grid%bounds
794 7140 : bo = grid%pw_grid%bounds_local
795 714 : grid2 => grid%array
796 714 : IF (par_scheme == do_par_atom) myind = 0
797 2010 : Radius: DO IRadTyp = 1, SIZE(pgfs)
798 1296 : pot => potentials(IRadTyp)%pot
799 1296 : dx = Pot%dx
800 1296 : 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 2010 : !$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 714 : CALL timestop(handle)
860 714 : END SUBROUTINE qmmm_elec_with_gaussian_LR
861 :
862 : END MODULE qmmm_gpw_energy
|
