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