(git:15c1bfc)
Loading...
Searching...
No Matches
qmmm_gpw_forces.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines to compute energy and forces in a QM/MM calculation
10!> \par History
11!> 05.2004 created [tlaino]
12!> \author Teodoro Laino
13! **************************************************************************************************
14MODULE qmmm_gpw_forces
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_print_key_finished_output,&
21 cp_print_key_unit_nr
22 USE cp_spline_utils, ONLY: pw_restrict_s3,&
23 spline3_nopbc_interp,&
24 spline3_pbc_interp
25 USE cube_utils, ONLY: cube_info_type
26 USE input_constants, ONLY: do_par_atom,&
27 do_qmmm_coulomb,&
28 do_qmmm_gauss,&
29 do_qmmm_none,&
30 do_qmmm_pcharge,&
31 do_qmmm_swave
32 USE input_section_types, ONLY: section_vals_get_subs_vals,&
33 section_vals_type,&
34 section_vals_val_get
35 USE kinds, ONLY: dp
36 USE message_passing, ONLY: mp_comm_type,&
37 mp_para_env_type,&
38 mp_request_type
39 USE mm_collocate_potential, ONLY: collocate_gf_rspace_nopbc,&
40 integrate_gf_rspace_nopbc
41 USE particle_types, ONLY: particle_type
42 USE pw_env_types, ONLY: pw_env_get,&
43 pw_env_type
44 USE pw_methods, ONLY: pw_axpy,&
45 pw_integral_ab,&
46 pw_transfer,&
47 pw_zero
48 USE pw_pool_types, ONLY: pw_pool_p_type,&
49 pw_pool_type,&
50 pw_pools_create_pws,&
51 pw_pools_give_back_pws
52 USE pw_types, ONLY: pw_c1d_gs_type,&
53 pw_r3d_rs_type
54 USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,&
55 qmmm_gaussian_type
56 USE qmmm_gpw_energy, ONLY: qmmm_elec_with_gaussian,&
57 qmmm_elec_with_gaussian_lg,&
58 qmmm_elec_with_gaussian_lr
59 USE qmmm_se_forces, ONLY: deriv_se_qmmm_matrix
60 USE qmmm_tb_methods, ONLY: deriv_tb_qmmm_matrix,&
61 deriv_tb_qmmm_matrix_pc
62 USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
63 qmmm_per_pot_p_type,&
64 qmmm_per_pot_type,&
65 qmmm_pot_p_type,&
66 qmmm_pot_type
67 USE qmmm_util, ONLY: spherical_cutoff_factor
68 USE qs_energy_types, ONLY: qs_energy_type
69 USE qs_environment_types, ONLY: get_qs_env,&
70 qs_environment_type
71 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
72 USE qs_rho_types, ONLY: qs_rho_get,&
73 qs_rho_type
74#include "./base/base_uses.f90"
75
76 IMPLICIT NONE
77
78 PRIVATE
79 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
80 REAL(KIND=dp), PARAMETER, PRIVATE :: dx = 0.01_dp ! Debug Variables
81 REAL(KIND=dp), PARAMETER, PRIVATE :: maxerr = 10.0_dp ! Debug Variables
82 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gpw_forces'
83 PUBLIC :: qmmm_forces
84
85CONTAINS
86
87! **************************************************************************************************
88!> \brief General driver to Compute the contribution
89!> to the forces due to the QM/MM potential
90!> \param qs_env ...
91!> \param qmmm_env ...
92!> \param mm_particles ...
93!> \param calc_force ...
94!> \param mm_cell ...
95!> \par History
96!> 06.2004 created [tlaino]
97!> \author Teodoro Laino
98! **************************************************************************************************
99 SUBROUTINE qmmm_forces(qs_env, qmmm_env, mm_particles, calc_force, mm_cell)
100 TYPE(qs_environment_type), POINTER :: qs_env
101 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
102 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
103 LOGICAL, INTENT(in), OPTIONAL :: calc_force
104 TYPE(cell_type), POINTER :: mm_cell
105
106 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces'
107
108 INTEGER :: handle, iatom, image_indmm, imm, indmm, &
109 ispin, iw
110 LOGICAL :: gapw, need_f, periodic
111 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges, &
112 forces_added_shells
113 TYPE(cp_logger_type), POINTER :: logger
114 TYPE(dft_control_type), POINTER :: dft_control
115 TYPE(mp_para_env_type), POINTER :: para_env
116 TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
117 TYPE(pw_env_type), POINTER :: pw_env
118 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
119 TYPE(pw_pool_type), POINTER :: auxbas_pool
120 TYPE(pw_r3d_rs_type) :: rho_tot_r, rho_tot_r2, rho_tot_r3
121 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
122 TYPE(qs_energy_type), POINTER :: energy
123 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
124 TYPE(qs_rho_type), POINTER :: rho
125 TYPE(section_vals_type), POINTER :: input_section, interp_section, &
126 print_section
127
128 CALL timeset(routinen, handle)
129 need_f = .true.
130 periodic = qmmm_env%periodic
131 IF (PRESENT(calc_force)) need_f = calc_force
132 NULLIFY (dft_control, ks_qmmm_env_loc, rho, pw_env, energy, forces, &
133 forces_added_charges, input_section, rho0_s_gs, rhoz_cneo_s_gs, rho_r)
134 CALL get_qs_env(qs_env=qs_env, &
135 rho=rho, &
136 rho_core=rho_core, &
137 pw_env=pw_env, &
138 energy=energy, &
139 para_env=para_env, &
140 input=input_section, &
141 rho0_s_gs=rho0_s_gs, &
142 rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
143 dft_control=dft_control)
144
145 CALL qs_rho_get(rho, rho_r=rho_r)
146
147 logger => cp_get_default_logger()
148 ks_qmmm_env_loc => qs_env%ks_qmmm_env
149 interp_section => section_vals_get_subs_vals(input_section, "QMMM%INTERPOLATOR")
150 print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
151 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
152 extension=".qmmmLog")
153 gapw = dft_control%qs_control%gapw
154 ! If forces are required allocate these temporary arrays
155 IF (need_f) THEN
156 ALLOCATE (forces(3, qmmm_env%num_mm_atoms))
157 ALLOCATE (forces_added_charges(3, qmmm_env%added_charges%num_mm_atoms))
158 ALLOCATE (forces_added_shells(3, qmmm_env%added_shells%num_mm_atoms))
159 forces(:, :) = 0.0_dp
160 forces_added_charges(:, :) = 0.0_dp
161 forces_added_shells(:, :) = 0.0_dp
162 END IF
163 IF (dft_control%qs_control%semi_empirical) THEN
164 ! SEMIEMPIRICAL
165 SELECT CASE (qmmm_env%qmmm_coupl_type)
166 CASE (do_qmmm_coulomb)
167 CALL deriv_se_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
168 need_f, forces, forces_added_charges)
169 CASE (do_qmmm_pcharge)
170 cpabort("Point Charge QM/MM electrostatic coupling not yet implemented for SE.")
171 CASE (do_qmmm_gauss, do_qmmm_swave)
172 cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.")
173 CASE (do_qmmm_none)
174 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
175 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
176 CASE DEFAULT
177 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
178 cpabort("")
179 END SELECT
180 ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
181 ! DFTB
182 SELECT CASE (qmmm_env%qmmm_coupl_type)
183 CASE (do_qmmm_none)
184 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
185 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
186 CASE (do_qmmm_coulomb)
187 CALL deriv_tb_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
188 need_f, forces, forces_added_charges)
189 CASE (do_qmmm_pcharge)
190 CALL deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
191 need_f, forces, forces_added_charges)
192 CASE (do_qmmm_gauss, do_qmmm_swave)
193 cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for DFTB.")
194 CASE DEFAULT
195 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
196 cpabort("")
197 END SELECT
198 IF (need_f) THEN
199 forces(:, :) = forces(:, :)/real(para_env%num_pe, kind=dp)
200 forces_added_charges(:, :) = forces_added_charges(:, :)/real(para_env%num_pe, kind=dp)
201 END IF
202 ELSE
203 ! GPW/GAPW
204 CALL pw_env_get(pw_env=pw_env, &
205 pw_pools=pw_pools, &
206 auxbas_pw_pool=auxbas_pool)
207 CALL auxbas_pool%create_pw(rho_tot_r)
208 ! IF GAPW the core charge is replaced by the compensation charge
209 IF (gapw) THEN
210 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
211 CALL pw_transfer(rho_core, rho_tot_r)
212 energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
213 CALL auxbas_pool%create_pw(rho_tot_r2)
214 CALL pw_transfer(rho0_s_gs, rho_tot_r2)
215 CALL pw_axpy(rho_tot_r2, rho_tot_r)
216 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
217 CALL auxbas_pool%create_pw(rho_tot_r3)
218 CALL pw_transfer(rhoz_cneo_s_gs, rho_tot_r3)
219 CALL pw_axpy(rho_tot_r3, rho_tot_r)
220 CALL auxbas_pool%give_back_pw(rho_tot_r3)
221 END IF
222 CALL auxbas_pool%give_back_pw(rho_tot_r2)
223 ELSE
224 CALL pw_transfer(rho0_s_gs, rho_tot_r)
225 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
226 CALL auxbas_pool%create_pw(rho_tot_r3)
227 CALL pw_transfer(rhoz_cneo_s_gs, rho_tot_r3)
228 CALL pw_axpy(rho_tot_r3, rho_tot_r)
229 CALL auxbas_pool%give_back_pw(rho_tot_r3)
230 END IF
231 !
232 ! QM/MM Nuclear Electrostatic Potential already included through rho0
233 !
234 energy%qmmm_nu = 0.0_dp
235 END IF
236 ELSE
237 CALL pw_transfer(rho_core, rho_tot_r)
238 !
239 ! Computes the QM/MM Nuclear Electrostatic Potential
240 !
241 energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
242 END IF
243 IF (need_f) THEN
244 !
245 DO ispin = 1, SIZE(rho_r)
246 CALL pw_axpy(rho_r(ispin), rho_tot_r)
247 END DO
248 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Evaluating forces on MM atoms due to the:"
249 ! Electrostatic Interaction type...
250 SELECT CASE (qmmm_env%qmmm_coupl_type)
251 CASE (do_qmmm_coulomb)
252 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
253 CASE (do_qmmm_pcharge)
254 cpabort("Point Charge QM/MM electrostatic coupling not yet implemented for GPW/GAPW.")
255 CASE (do_qmmm_gauss, do_qmmm_swave)
256 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
257 "- QM/MM Coupling computed collocating the Gaussian Potential Functions."
258 CALL qmmm_forces_with_gaussian(rho=rho_tot_r, &
259 qmmm_env=qmmm_env, &
260 mm_particles=mm_particles, &
261 aug_pools=qmmm_env%aug_pools, &
262 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
263 coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
264 para_env=para_env, &
265 pw_pools=pw_pools, &
266 eps_mm_rspace=qmmm_env%eps_mm_rspace, &
267 cube_info=ks_qmmm_env_loc%cube_info, &
268 forces=forces, &
269 forces_added_charges=forces_added_charges, &
270 forces_added_shells=forces_added_shells, &
271 interp_section=interp_section, &
272 iw=iw, &
273 mm_cell=mm_cell)
274 CASE (do_qmmm_none)
275 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
276 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
277 CASE DEFAULT
278 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "- Unknown Coupling..."
279 cpabort("")
280 END SELECT
281 END IF
282 END IF
283 ! Correct Total Energy adding the contribution of the QM/MM nuclear interaction
284 energy%total = energy%total + energy%qmmm_nu
285 ! Proceed if gradients are requested..
286 IF (need_f) THEN
287 !ikuo Temporary change to alleviate compiler problems on Intel with
288 !array dimension of 0
289 IF (qmmm_env%num_mm_atoms /= 0) CALL para_env%sum(forces)
290 IF (qmmm_env%added_charges%num_mm_atoms /= 0) CALL para_env%sum(forces_added_charges)
291 IF (qmmm_env%added_shells%num_mm_atoms /= 0) CALL para_env%sum(forces_added_shells)
292 ! Debug Forces
293 IF (debug_this_module) THEN
294 IF (dft_control%qs_control%semi_empirical .OR. &
295 dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
296 WRITE (iw, *) "NO DEBUG AVAILABLE in module"//trim(routinen)
297 ELSE
298 ! Print Out Forces
299 IF (iw > 0) THEN
300 DO imm = 1, SIZE(qmmm_env%mm_atom_index)
301 WRITE (iw, *) "ANALYTICAL FORCES:"
302 indmm = qmmm_env%mm_atom_index(imm)
303 WRITE (iw, '(I6,3F15.9)') indmm, forces(:, imm)
304 END DO
305 END IF
306 CALL qmmm_debug_forces(rho=rho_tot_r, &
307 qs_env=qs_env, &
308 qmmm_env=qmmm_env, &
309 analytical_forces=forces, &
310 mm_particles=mm_particles, &
311 mm_atom_index=qmmm_env%mm_atom_index, &
312 num_mm_atoms=qmmm_env%num_mm_atoms, &
313 interp_section=interp_section, &
314 mm_cell=mm_cell)
315 END IF
316 END IF
317 END IF
318 ! Give back rho_tot_t to auxbas_pool only for GPW/GAPW
319 IF ((.NOT. dft_control%qs_control%semi_empirical) .AND. &
320 (.NOT. dft_control%qs_control%dftb) .AND. (.NOT. dft_control%qs_control%xtb)) THEN
321 CALL auxbas_pool%give_back_pw(rho_tot_r)
322 END IF
323 IF (iw > 0) THEN
324 IF (.NOT. gapw) WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
325 "QM/MM Nuclear Electrostatic Potential :", energy%qmmm_nu
326 WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
327 "QMMM Total Energy (QM + QMMM electronic + QMMM nuclear):", energy%total
328 WRITE (iw, '(T2,"QMMM|",1X,A)') "MM energy NOT included in the above term!"// &
329 " Check for: FORCE_EVAL ( QMMM )"
330 WRITE (iw, '(T2,"QMMM|",1X,A)') "that includes both QM, QMMM and MM energy terms!"
331 END IF
332 IF (need_f) THEN
333 ! Transfer Forces
334 DO imm = 1, qmmm_env%num_mm_atoms
335 indmm = qmmm_env%mm_atom_index(imm)
336
337 !add image forces to Forces
338 IF (qmmm_env%image_charge) THEN
339 DO iatom = 1, qmmm_env%num_image_mm_atoms
340 image_indmm = qmmm_env%image_charge_pot%image_mm_list(iatom)
341 IF (image_indmm .EQ. indmm) THEN
342 forces(:, imm) = forces(:, imm) &
343 + qmmm_env%image_charge_pot%image_forcesMM(:, iatom)
344 END IF
345 END DO
346 END IF
347
348 ! Hack: In Forces there the gradients indeed...
349 ! Minux sign to take care of this misunderstanding...
350 mm_particles(indmm)%f(:) = -forces(:, imm) + mm_particles(indmm)%f(:)
351 END DO
352 DEALLOCATE (forces)
353 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
354 DO imm = 1, qmmm_env%added_charges%num_mm_atoms
355 indmm = qmmm_env%added_charges%mm_atom_index(imm)
356 ! Hack: In Forces there the gradients indeed...
357 ! Minux sign to take care of this misunderstanding...
358 qmmm_env%added_charges%added_particles(indmm)%f(:) = -forces_added_charges(:, imm)
359 END DO
360 END IF
361 DEALLOCATE (forces_added_charges)
362 IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN
363 DO imm = 1, qmmm_env%added_shells%num_mm_atoms
364 indmm = qmmm_env%added_shells%mm_core_index(imm)
365 ! Hack: In Forces there the gradients indeed...
366 ! Minux sign to take care of this misunderstanding...
367 qmmm_env%added_shells%added_particles(imm)%f(:) = qmmm_env%added_shells%added_particles(imm)%f(:) - &
368 forces_added_shells(:, imm)
369
370 END DO
371 END IF
372 DEALLOCATE (forces_added_shells)
373 END IF
374 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
375 CALL timestop(handle)
376
377 END SUBROUTINE qmmm_forces
378
379! **************************************************************************************************
380!> \brief Evaluates the contribution to the forces due to the
381!> QM/MM potential computed collocating the Electrostatic
382!> Gaussian Potential.
383!> \param rho ...
384!> \param qmmm_env ...
385!> \param mm_particles ...
386!> \param aug_pools ...
387!> \param auxbas_grid ...
388!> \param coarser_grid ...
389!> \param cube_info ...
390!> \param para_env ...
391!> \param eps_mm_rspace ...
392!> \param pw_pools ...
393!> \param Forces ...
394!> \param Forces_added_charges ...
395!> \param Forces_added_shells ...
396!> \param interp_section ...
397!> \param iw ...
398!> \param mm_cell ...
399!> \par History
400!> 06.2004 created [tlaino]
401!> \author Teodoro Laino
402! **************************************************************************************************
403 SUBROUTINE qmmm_forces_with_gaussian(rho, qmmm_env, mm_particles, &
404 aug_pools, auxbas_grid, coarser_grid, cube_info, para_env, &
405 eps_mm_rspace, pw_pools, Forces, Forces_added_charges, Forces_added_shells, &
406 interp_section, iw, mm_cell)
407 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
408 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
409 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
410 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
411 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
412 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
413 TYPE(mp_para_env_type), POINTER :: para_env
414 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
415 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
416 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges, &
417 forces_added_shells
418 TYPE(section_vals_type), POINTER :: interp_section
419 INTEGER, INTENT(IN) :: iw
420 TYPE(cell_type), POINTER :: mm_cell
421
422 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian'
423
424 INTEGER :: handle, i, igrid, j, k, kind_interp, me, &
425 ngrids
426 INTEGER, DIMENSION(3) :: glb, gub, lb, ub
427 INTEGER, DIMENSION(:), POINTER :: pos_of_x
428 LOGICAL :: shells
429 REAL(kind=dp), DIMENSION(:, :), POINTER :: tmp
430 TYPE(mp_comm_type) :: group
431 TYPE(mp_request_type) :: request
432 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
433
434! Statements
435
436 CALL timeset(routinen, handle)
437 NULLIFY (tmp)
438 cpassert(ASSOCIATED(mm_particles))
439 cpassert(ASSOCIATED(qmmm_env%mm_atom_chrg))
440 cpassert(ASSOCIATED(qmmm_env%mm_atom_index))
441 cpassert(ASSOCIATED(forces))
442 !Statements
443 ngrids = SIZE(pw_pools)
444 CALL pw_pools_create_pws(aug_pools, grids)
445 DO igrid = 1, ngrids
446 CALL pw_zero(grids(igrid))
447 END DO
448 ! Collocate Density on multigrids
449 lb = rho%pw_grid%bounds_local(1, :)
450 ub = rho%pw_grid%bounds_local(2, :)
451 grids(auxbas_grid)%array(lb(1):ub(1), &
452 lb(2):ub(2), &
453 lb(3):ub(3)) = rho%array
454 ! copy the boundaries
455 DO i = lb(1), ub(1)
456 grids(auxbas_grid)%array(i, ub(2) + 1, ub(3) + 1) = rho%array(i, lb(2), lb(3))
457 END DO
458 DO k = lb(3), ub(3)
459 DO i = lb(1), ub(1)
460 grids(auxbas_grid)%array(i, ub(2) + 1, k) = rho%array(i, lb(2), k)
461 END DO
462 END DO
463 DO j = lb(2), ub(2)
464 DO i = lb(1), ub(1)
465 grids(auxbas_grid)%array(i, j, ub(3) + 1) = rho%array(i, j, lb(3))
466 END DO
467 END DO
468 pos_of_x => grids(auxbas_grid)%pw_grid%para%pos_of_x
469 group = grids(auxbas_grid)%pw_grid%para%group
470 me = grids(auxbas_grid)%pw_grid%para%group%mepos
471 glb = rho%pw_grid%bounds(1, :)
472 gub = rho%pw_grid%bounds(2, :)
473 IF ((pos_of_x(glb(1)) .EQ. me) .AND. (pos_of_x(gub(1)) .EQ. me)) THEN
474 DO k = lb(3), ub(3)
475 DO j = lb(2), ub(2)
476 grids(auxbas_grid)%array(ub(1) + 1, j, k) = rho%array(lb(1), j, k)
477 END DO
478 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, k) = rho%array(lb(1), lb(2), k)
479 END DO
480 DO j = lb(2), ub(2)
481 grids(auxbas_grid)%array(ub(1) + 1, j, ub(3) + 1) = rho%array(lb(1), j, lb(3))
482 END DO
483 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, ub(3) + 1) = rho%array(lb(1), lb(2), lb(3))
484 ELSE IF (pos_of_x(glb(1)) .EQ. me) THEN
485 ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
486 rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)))
487 tmp = rho%array(lb(1), :, :)
488 CALL group%isend(msgin=tmp, dest=pos_of_x(rho%pw_grid%bounds(2, 1)), &
489 request=request, tag=112)
490 CALL request%wait()
491 ELSE IF (pos_of_x(gub(1)) .EQ. me) THEN
492 ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
493 rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)))
494 CALL group%irecv(msgout=tmp, source=pos_of_x(rho%pw_grid%bounds(1, 1)), &
495 request=request, tag=112)
496 CALL request%wait()
497
498 DO k = lb(3), ub(3)
499 DO j = lb(2), ub(2)
500 grids(auxbas_grid)%array(ub(1) + 1, j, k) = tmp(j, k)
501 END DO
502 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, k) = tmp(lb(2), k)
503 END DO
504 DO j = lb(2), ub(2)
505 grids(auxbas_grid)%array(ub(1) + 1, j, ub(3) + 1) = tmp(j, lb(3))
506 END DO
507 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, ub(3) + 1) = tmp(lb(2), lb(3))
508 END IF
509 IF (ASSOCIATED(tmp)) THEN
510 DEALLOCATE (tmp)
511 END IF
512 ! Further setup of parallelization scheme
513 IF (qmmm_env%par_scheme == do_par_atom) THEN
514 CALL para_env%sum(grids(auxbas_grid)%array)
515 END IF
516 ! RealSpace Interpolation
517 CALL section_vals_val_get(interp_section, "kind", i_val=kind_interp)
518 SELECT CASE (kind_interp)
519 CASE (spline3_nopbc_interp, spline3_pbc_interp)
520 ! Spline Interpolator
521 DO igrid = auxbas_grid, SIZE(grids) - 1
522 CALL pw_restrict_s3(grids(igrid), &
523 grids(igrid + 1), &
524 aug_pools(igrid + 1)%pool, &
525 param_section=interp_section)
526 END DO
527 CASE DEFAULT
528 cpabort("")
529 END SELECT
530
531 shells = .false.
532 CALL qmmm_force_with_gaussian_low(grids, mm_particles, &
533 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
534 qmmm_env%num_mm_atoms, cube_info, para_env, eps_mm_rspace, auxbas_grid, &
535 coarser_grid, qmmm_env%pgfs, qmmm_env%potentials, forces, aug_pools, &
536 mm_cell, qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%per_potentials, &
537 iw, qmmm_env%par_scheme, qmmm_env%spherical_cutoff, shells)
538
539 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
540 CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_charges%added_particles, &
541 qmmm_env%added_charges%mm_atom_chrg, &
542 qmmm_env%added_charges%mm_atom_index, qmmm_env%added_charges%num_mm_atoms, &
543 cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_charges%pgfs, &
544 qmmm_env%added_charges%potentials, forces_added_charges, aug_pools, mm_cell, &
545 qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%added_charges%per_potentials, iw, qmmm_env%par_scheme, &
546 qmmm_env%spherical_cutoff, shells)
547 END IF
548
549 IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN
550 shells = .true.
551 CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_shells%added_particles, &
552 qmmm_env%added_shells%mm_core_chrg, &
553 qmmm_env%added_shells%mm_core_index, qmmm_env%added_shells%num_mm_atoms, &
554 cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_shells%pgfs, &
555 qmmm_env%added_shells%potentials, forces_added_shells, aug_pools, mm_cell, &
556 qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%added_shells%per_potentials, iw, qmmm_env%par_scheme, &
557 qmmm_env%spherical_cutoff, shells)
558 END IF
559
560 CALL pw_pools_give_back_pws(aug_pools, grids)
561 CALL timestop(handle)
562
563 END SUBROUTINE qmmm_forces_with_gaussian
564
565! **************************************************************************************************
566!> \brief Evaluates the contribution to the forces due to the
567!> QM/MM potential computed collocating the Electrostatic
568!> Gaussian Potential. Low Level
569!> \param grids ...
570!> \param mm_particles ...
571!> \param mm_charges ...
572!> \param mm_atom_index ...
573!> \param num_mm_atoms ...
574!> \param cube_info ...
575!> \param para_env ...
576!> \param eps_mm_rspace ...
577!> \param auxbas_grid ...
578!> \param coarser_grid ...
579!> \param pgfs ...
580!> \param potentials ...
581!> \param Forces ...
582!> \param aug_pools ...
583!> \param mm_cell ...
584!> \param dOmmOqm ...
585!> \param periodic ...
586!> \param per_potentials ...
587!> \param iw ...
588!> \param par_scheme ...
589!> \param qmmm_spherical_cutoff ...
590!> \param shells ...
591!> \par History
592!> 06.2004 created [tlaino]
593!> \author Teodoro Laino
594! **************************************************************************************************
595 SUBROUTINE qmmm_force_with_gaussian_low(grids, mm_particles, mm_charges, &
596 mm_atom_index, num_mm_atoms, cube_info, para_env, &
597 eps_mm_rspace, auxbas_grid, coarser_grid, pgfs, potentials, Forces, &
598 aug_pools, mm_cell, dOmmOqm, periodic, per_potentials, iw, par_scheme, &
599 qmmm_spherical_cutoff, shells)
600 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: grids
601 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
602 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
603 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
604 INTEGER, INTENT(IN) :: num_mm_atoms
605 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
606 TYPE(mp_para_env_type), POINTER :: para_env
607 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
608 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
609 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
610 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
611 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
612 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
613 TYPE(cell_type), POINTER :: mm_cell
614 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
615 LOGICAL, INTENT(in) :: periodic
616 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
617 INTEGER, INTENT(IN) :: iw, par_scheme
618 REAL(kind=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
619 LOGICAL, INTENT(in) :: shells
620
621 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_force_with_gaussian_low', &
622 routinenb = 'qmmm_forces_gaussian_low'
623
624 INTEGER :: handle, handle2, igauss, ilevel, imm, &
625 indmm, iradtyp, lindmm, myind, &
626 n_rep_real(3)
627 INTEGER, DIMENSION(2, 3) :: bo
628 REAL(kind=dp) :: alpha, dvol, height, sph_chrg_factor, w
629 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xdat, ydat, zdat
630 REAL(kind=dp), DIMENSION(3) :: force, ra
631 TYPE(qmmm_gaussian_type), POINTER :: pgf
632 TYPE(qmmm_per_pot_type), POINTER :: per_pot
633 TYPE(qmmm_pot_type), POINTER :: pot
634
635 CALL timeset(routinen, handle)
636 CALL timeset(routinenb//"_G", handle2)
637 NULLIFY (pgf, pot, per_pot)
638 IF (par_scheme == do_par_atom) myind = 0
639 radius: DO iradtyp = 1, SIZE(pgfs)
640 pgf => pgfs(iradtyp)%pgf
641 pot => potentials(iradtyp)%pot
642 n_rep_real = 0
643 IF (periodic) THEN
644 per_pot => per_potentials(iradtyp)%pot
645 n_rep_real = per_pot%n_rep_real
646 END IF
647 gaussian: DO igauss = 1, pgf%Number_of_Gaussians
648 alpha = 1.0_dp/pgf%Gk(igauss)
649 alpha = alpha*alpha
650 height = pgf%Ak(igauss)
651 ilevel = pgf%grid_level(igauss)
652 dvol = grids(ilevel)%pw_grid%dvol
653 bo = grids(ilevel)%pw_grid%bounds_local
654 ALLOCATE (xdat(2, bo(1, 1):bo(2, 1)))
655 ALLOCATE (ydat(2, bo(1, 2):bo(2, 2)))
656 ALLOCATE (zdat(2, bo(1, 3):bo(2, 3)))
657 !$OMP PARALLEL DO DEFAULT(NONE) &
658 !$OMP SHARED(pot, par_scheme, dvol, alpha, para_env, mm_atom_index, shells) &
659 !$OMP SHARED(mm_particles, dOmmOqm, mm_cell, height, mm_charges, qmmm_spherical_cutoff) &
660 !$OMP SHARED(grids, cube_info, bo, n_rep_real, eps_mm_rspace, Forces, ilevel) &
661 !$OMP SHARED(IGauss, pgf, IRadTyp, iw, aug_pools, auxbas_grid) &
662 !$OMP PRIVATE(xdat, ydat, zdat) &
663 !$OMP PRIVATE(Imm, LIndMM, IndMM, ra, W, force, sph_chrg_factor, myind)
664 atoms: DO imm = 1, SIZE(pot%mm_atom_index)
665 IF (par_scheme == do_par_atom) THEN
666 myind = imm + (igauss - 1)*SIZE(pot%mm_atom_index) + (iradtyp - 1)*pgf%Number_of_Gaussians
667 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
668 END IF
669 lindmm = pot%mm_atom_index(imm)
670 indmm = mm_atom_index(lindmm)
671 IF (shells) THEN
672 ra(:) = pbc(mm_particles(imm)%r - dommoqm, mm_cell) + dommoqm
673 ELSE
674 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
675 END IF
676 w = mm_charges(lindmm)*height
677 force = 0.0_dp
678 ! Possible Spherical Cutoff
679 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
680 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
681 w = w*sph_chrg_factor
682 END IF
683 IF (abs(w) <= epsilon(0.0_dp)) cycle atoms
684 CALL integrate_gf_rspace_nopbc(zetp=alpha, &
685 rp=ra, &
686 scale=-1.0_dp, &
687 w=w, &
688 pwgrid=grids(ilevel), &
689 cube_info=cube_info(ilevel), &
690 eps_mm_rspace=eps_mm_rspace, &
691 xdat=xdat, &
692 ydat=ydat, &
693 zdat=zdat, &
694 bo=bo, &
695 force=force, &
696 n_rep_real=n_rep_real, &
697 mm_cell=mm_cell)
698 force = force*dvol
699 forces(:, lindmm) = forces(:, lindmm) + force(:)
700 !
701 ! Debug Statement
702 !
703 IF (debug_this_module) THEN
704 CALL debug_integrate_gf_rspace_nopbc(ilevel=ilevel, &
705 zetp=alpha, &
706 rp=ra, &
707 w=w, &
708 pwgrid=grids(ilevel), &
709 cube_info=cube_info(ilevel), &
710 eps_mm_rspace=eps_mm_rspace, &
711 aug_pools=aug_pools, &
712 debug_force=force, &
713 mm_cell=mm_cell, &
714 auxbas_grid=auxbas_grid, &
715 n_rep_real=n_rep_real, &
716 iw=iw)
717 END IF
718 END DO atoms
719 !$OMP END PARALLEL DO
720 DEALLOCATE (xdat)
721 DEALLOCATE (ydat)
722 DEALLOCATE (zdat)
723 END DO gaussian
724 END DO radius
725 CALL timestop(handle2)
726 CALL timeset(routinenb//"_R", handle2)
727 IF (periodic) THEN
728 CALL qmmm_forces_with_gaussian_lg(pgfs=pgfs, &
729 cgrid=grids(coarser_grid), &
730 num_mm_atoms=num_mm_atoms, &
731 mm_charges=mm_charges, &
732 mm_atom_index=mm_atom_index, &
733 mm_particles=mm_particles, &
734 para_env=para_env, &
735 coarser_grid_level=coarser_grid, &
736 forces=forces, &
737 per_potentials=per_potentials, &
738 aug_pools=aug_pools, &
739 mm_cell=mm_cell, &
740 dommoqm=dommoqm, &
741 iw=iw, &
742 par_scheme=par_scheme, &
743 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
744 shells=shells)
745 ELSE
746 CALL qmmm_forces_with_gaussian_lr(pgfs=pgfs, &
747 cgrid=grids(coarser_grid), &
748 num_mm_atoms=num_mm_atoms, &
749 mm_charges=mm_charges, &
750 mm_atom_index=mm_atom_index, &
751 mm_particles=mm_particles, &
752 para_env=para_env, &
753 coarser_grid_level=coarser_grid, &
754 forces=forces, &
755 potentials=potentials, &
756 aug_pools=aug_pools, &
757 mm_cell=mm_cell, &
758 dommoqm=dommoqm, &
759 iw=iw, &
760 par_scheme=par_scheme, &
761 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
762 shells=shells)
763 END IF
764 CALL timestop(handle2)
765 CALL timestop(handle)
766 END SUBROUTINE qmmm_force_with_gaussian_low
767
768! **************************************************************************************************
769!> \brief Evaluates the contribution to the forces due to the Long Range
770!> part of the QM/MM potential computed collocating the Electrostatic
771!> Gaussian Potential.
772!> \param pgfs ...
773!> \param cgrid ...
774!> \param num_mm_atoms ...
775!> \param mm_charges ...
776!> \param mm_atom_index ...
777!> \param mm_particles ...
778!> \param para_env ...
779!> \param coarser_grid_level ...
780!> \param Forces ...
781!> \param per_potentials ...
782!> \param aug_pools ...
783!> \param mm_cell ...
784!> \param dOmmOqm ...
785!> \param iw ...
786!> \param par_scheme ...
787!> \param qmmm_spherical_cutoff ...
788!> \param shells ...
789!> \par History
790!> 08.2004 created [tlaino]
791!> \author Teodoro Laino
792! **************************************************************************************************
793 SUBROUTINE qmmm_forces_with_gaussian_lg(pgfs, cgrid, num_mm_atoms, mm_charges, mm_atom_index, &
794 mm_particles, para_env, coarser_grid_level, Forces, per_potentials, &
795 aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
796 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
797 TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
798 INTEGER, INTENT(IN) :: num_mm_atoms
799 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
800 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
801 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
802 TYPE(mp_para_env_type), POINTER :: para_env
803 INTEGER, INTENT(IN) :: coarser_grid_level
804 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
805 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
806 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
807 TYPE(cell_type), POINTER :: mm_cell
808 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
809 INTEGER, INTENT(IN) :: iw, par_scheme
810 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
811 LOGICAL :: shells
812
813 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian_LG'
814
815 INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, ij3, ij4, ik1, ik2, ik3, ik4, imm, &
816 indmm, iradtyp, ivec(3), j, k, lindmm, my_i, my_j, my_k, myind, npts(3)
817 INTEGER, DIMENSION(2, 3) :: bo, gbo
818 REAL(kind=dp) :: a1, a2, a3, abc_x(4, 4), abc_x_y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, &
819 dr1, dr1c, dr1i, dr2, dr2c, dr2i, dr3, dr3c, dr3i, dvol, e1, e2, e3, f1, f2, f3, fac, &
820 ft1, ft2, ft3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, &
821 rt3, rv1, rv2, rv3, s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, &
822 sph_chrg_factor, t1, t1d, t1o, t2, t2d, t2o, t3, t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, &
823 v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, v4o, xd1, xd2, xd3, xs1, xs2, xs3
824 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: lforces
825 REAL(kind=dp), DIMENSION(3) :: ra, val, vec
826 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid, grid2
827 TYPE(pw_r3d_rs_type), POINTER :: pw
828 TYPE(qmmm_per_pot_type), POINTER :: per_pot
829
830 CALL timeset(routinen, handle)
831 NULLIFY (grid)
832 ALLOCATE (lforces(3, num_mm_atoms))
833 lforces = 0.0_dp
834 dr1c = cgrid%pw_grid%dr(1)
835 dr2c = cgrid%pw_grid%dr(2)
836 dr3c = cgrid%pw_grid%dr(3)
837 dvol = cgrid%pw_grid%dvol
838 gbo = cgrid%pw_grid%bounds
839 bo = cgrid%pw_grid%bounds_local
840 grid => cgrid%array
841 IF (par_scheme == do_par_atom) myind = 0
842 radius: DO iradtyp = 1, SIZE(pgfs)
843 per_pot => per_potentials(iradtyp)%pot
844 pw => per_pot%TabLR
845 grid2 => pw%array(:, :, :)
846 npts = pw%pw_grid%npts
847 dr1 = pw%pw_grid%dr(1)
848 dr2 = pw%pw_grid%dr(2)
849 dr3 = pw%pw_grid%dr(3)
850 dr1i = 1.0_dp/dr1
851 dr2i = 1.0_dp/dr2
852 dr3i = 1.0_dp/dr3
853
854 !$OMP PARALLEL DO DEFAULT(NONE) &
855 !$OMP SHARED(bo, grid, grid2, pw, npts, gbo, per_pot, mm_atom_index) &
856 !$OMP SHARED(dr1, dr2, dr3, dr1i, dr2i, dr3i, dr1c, dr2c, dr3c, par_scheme, mm_charges) &
857 !$OMP SHARED(mm_cell, dOmmOqm, dvol, shells, para_env, IRadTyp) &
858 !$OMP SHARED(qmmm_spherical_cutoff, mm_particles, Forces, LForces) &
859 !$OMP PRIVATE(qt, Imm, LIndMM, IndMM, sph_chrg_factor, ra, myind) &
860 !$OMP PRIVATE(rt1, rt2, rt3, ft1, ft2, ft3, my_k, my_j, my_i, xs3, xs2, xs1) &
861 !$OMP PRIVATE(rv3, rv2, rv1, vec, ivec, ik1, ik2, ik3, ik4, xd3, xd2, xd1) &
862 !$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, u1, u2, u3, v1o, v2o, v3o, v4o) &
863 !$OMP PRIVATE(v1d, v2d, v3d, v4d, ij1, ij2, ij3, ij4, e1, e2, e3, f1, f2, f3) &
864 !$OMP PRIVATE(g1, g2, g3, h1, h2, h3, s1o, s2o, s3o, s4o, s1d, s2d, s3d, s4d) &
865 !$OMP PRIVATE(ii1, ii2, ii3, ii4, a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3) &
866 !$OMP PRIVATE(t1o, t2o, t3o, t4o, t1d, t2d, t3d, t4d, t1, t2, t3, t4, s1, s2, s3, s4) &
867 !$OMP PRIVATE(v1, v2, v3, v4, abc_x, abc_x_y, val, fac)
868 atoms: DO imm = 1, SIZE(per_pot%mm_atom_index)
869 IF (par_scheme == do_par_atom) THEN
870 myind = imm + (iradtyp - 1)*SIZE(per_pot%mm_atom_index)
871 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle
872 END IF
873 lindmm = per_pot%mm_atom_index(imm)
874 indmm = mm_atom_index(lindmm)
875 IF (shells) THEN
876 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
877 ELSE
878 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
879 END IF
880 qt = mm_charges(lindmm)
881 ! Possible Spherical Cutoff
882 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
883 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
884 qt = qt*sph_chrg_factor
885 END IF
886 IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
887 rt1 = ra(1)
888 rt2 = ra(2)
889 rt3 = ra(3)
890 ft1 = 0.0_dp
891 ft2 = 0.0_dp
892 ft3 = 0.0_dp
893 loopongrid: DO k = bo(1, 3), bo(2, 3)
894 my_k = k - gbo(1, 3)
895 xs3 = real(my_k, dp)*dr3c
896 my_j = bo(1, 2) - gbo(1, 2)
897 xs2 = real(my_j, dp)*dr2c
898 rv3 = rt3 - xs3
899 vec(3) = rv3
900 ivec(3) = floor(vec(3)/pw%pw_grid%dr(3))
901 ik1 = modulo(ivec(3) - 1, npts(3)) + 1
902 ik2 = modulo(ivec(3), npts(3)) + 1
903 ik3 = modulo(ivec(3) + 1, npts(3)) + 1
904 ik4 = modulo(ivec(3) + 2, npts(3)) + 1
905 xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
906 p1 = 3.0_dp + xd3
907 p2 = p1*p1
908 p3 = p2*p1
909 q1 = 2.0_dp + xd3
910 q2 = q1*q1
911 q3 = q2*q1
912 r1 = 1.0_dp + xd3
913 r2 = r1*r1
914 r3 = r2*r1
915 u1 = xd3
916 u2 = u1*u1
917 u3 = u2*u1
918 v1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
919 v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
920 v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
921 v4o = 1.0_dp/6.0_dp*u3
922 v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2
923 v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2
924 v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2
925 v4d = 0.5_dp*u2
926 DO j = bo(1, 2), bo(2, 2)
927 my_i = bo(1, 1) - gbo(1, 1)
928 xs1 = real(my_i, dp)*dr1c
929 rv2 = rt2 - xs2
930 vec(2) = rv2
931 ivec(2) = floor(vec(2)/pw%pw_grid%dr(2))
932 ij1 = modulo(ivec(2) - 1, npts(2)) + 1
933 ij2 = modulo(ivec(2), npts(2)) + 1
934 ij3 = modulo(ivec(2) + 1, npts(2)) + 1
935 ij4 = modulo(ivec(2) + 2, npts(2)) + 1
936 xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
937 e1 = 3.0_dp + xd2
938 e2 = e1*e1
939 e3 = e2*e1
940 f1 = 2.0_dp + xd2
941 f2 = f1*f1
942 f3 = f2*f1
943 g1 = 1.0_dp + xd2
944 g2 = g1*g1
945 g3 = g2*g1
946 h1 = xd2
947 h2 = h1*h1
948 h3 = h2*h1
949 s1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
950 s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
951 s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
952 s4o = 1.0_dp/6.0_dp*h3
953 s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2
954 s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2
955 s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2
956 s4d = 0.5_dp*h2
957 DO i = bo(1, 1), bo(2, 1)
958 rv1 = rt1 - xs1
959 vec(1) = rv1
960 ivec(1) = floor(vec(1)/pw%pw_grid%dr(1))
961 ii1 = modulo(ivec(1) - 1, npts(1)) + 1
962 ii2 = modulo(ivec(1), npts(1)) + 1
963 ii3 = modulo(ivec(1) + 1, npts(1)) + 1
964 ii4 = modulo(ivec(1) + 2, npts(1)) + 1
965 xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
966 a1 = 3.0_dp + xd1
967 a2 = a1*a1
968 a3 = a2*a1
969 b1 = 2.0_dp + xd1
970 b2 = b1*b1
971 b3 = b2*b1
972 c1 = 1.0_dp + xd1
973 c2 = c1*c1
974 c3 = c2*c1
975 d1 = xd1
976 d2 = d1*d1
977 d3 = d2*d1
978 t1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
979 t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
980 t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
981 t4o = 1.0_dp/6.0_dp*d3
982 t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2
983 t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2
984 t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2
985 t4d = 0.5_dp*d2
986
987 t1 = t1d*dr1i
988 t2 = t2d*dr1i
989 t3 = t3d*dr1i
990 t4 = t4d*dr1i
991 s1 = s1o
992 s2 = s2o
993 s3 = s3o
994 s4 = s4o
995 v1 = v1o
996 v2 = v2o
997 v3 = v3o
998 v4 = v4o
999
1000 abc_x(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
1001 abc_x(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
1002 abc_x(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
1003 abc_x(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
1004 abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1005
1006 abc_x(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
1007 abc_x(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
1008 abc_x(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
1009 abc_x(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
1010 abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1011
1012 abc_x(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
1013 abc_x(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
1014 abc_x(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
1015 abc_x(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
1016 abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1017
1018 abc_x(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
1019 abc_x(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
1020 abc_x(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
1021 abc_x(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
1022 abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1023
1024 val(1) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1025
1026 t1 = t1o
1027 t2 = t2o
1028 t3 = t3o
1029 t4 = t4o
1030 s1 = s1d*dr2i
1031 s2 = s2d*dr2i
1032 s3 = s3d*dr2i
1033 s4 = s4d*dr2i
1034
1035 abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1036 abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1037 abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1038 abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1039
1040 val(2) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1041
1042 t1 = t1o
1043 t2 = t2o
1044 t3 = t3o
1045 t4 = t4o
1046 s1 = s1o
1047 s2 = s2o
1048 s3 = s3o
1049 s4 = s4o
1050 v1 = v1d*dr3i
1051 v2 = v2d*dr3i
1052 v3 = v3d*dr3i
1053 v4 = v4d*dr3i
1054
1055 abc_x(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
1056 abc_x(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
1057 abc_x(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
1058 abc_x(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
1059 abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1060 abc_x(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
1061 abc_x(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
1062 abc_x(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
1063 abc_x(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
1064 abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1065 abc_x(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
1066 abc_x(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
1067 abc_x(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
1068 abc_x(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
1069 abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1070 abc_x(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
1071 abc_x(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
1072 abc_x(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
1073 abc_x(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
1074 abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1075
1076 val(3) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1077
1078 fac = grid(i, j, k)
1079 ft1 = ft1 + val(1)*fac
1080 ft2 = ft2 + val(2)*fac
1081 ft3 = ft3 + val(3)*fac
1082 xs1 = xs1 + dr1c
1083 END DO
1084 xs2 = xs2 + dr2c
1085 END DO
1086 END DO loopongrid
1087 qt = -qt*dvol
1088 lforces(1, lindmm) = ft1*qt
1089 lforces(2, lindmm) = ft2*qt
1090 lforces(3, lindmm) = ft3*qt
1091
1092 forces(1, lindmm) = forces(1, lindmm) + lforces(1, lindmm)
1093 forces(2, lindmm) = forces(2, lindmm) + lforces(2, lindmm)
1094 forces(3, lindmm) = forces(3, lindmm) + lforces(3, lindmm)
1095 END DO atoms
1096 !$OMP END PARALLEL DO
1097 END DO radius
1098 !
1099 ! Debug Statement
1100 !
1101 IF (debug_this_module) THEN
1102 CALL debug_qmmm_forces_with_gauss_lg(pgfs=pgfs, &
1103 aug_pools=aug_pools, &
1104 rho=cgrid, &
1105 num_mm_atoms=num_mm_atoms, &
1106 mm_charges=mm_charges, &
1107 mm_atom_index=mm_atom_index, &
1108 mm_particles=mm_particles, &
1109 coarser_grid_level=coarser_grid_level, &
1110 debug_force=lforces, &
1111 per_potentials=per_potentials, &
1112 para_env=para_env, &
1113 mm_cell=mm_cell, &
1114 dommoqm=dommoqm, &
1115 iw=iw, &
1116 par_scheme=par_scheme, &
1117 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1118 shells=shells)
1119 END IF
1120 DEALLOCATE (lforces)
1121 CALL timestop(handle)
1122 END SUBROUTINE qmmm_forces_with_gaussian_lg
1123
1124! **************************************************************************************************
1125!> \brief Evaluates the contribution to the forces due to the Long Range
1126!> part of the QM/MM potential computed collocating the Electrostatic
1127!> Gaussian Potential.
1128!> \param pgfs ...
1129!> \param cgrid ...
1130!> \param num_mm_atoms ...
1131!> \param mm_charges ...
1132!> \param mm_atom_index ...
1133!> \param mm_particles ...
1134!> \param para_env ...
1135!> \param coarser_grid_level ...
1136!> \param Forces ...
1137!> \param potentials ...
1138!> \param aug_pools ...
1139!> \param mm_cell ...
1140!> \param dOmmOqm ...
1141!> \param iw ...
1142!> \param par_scheme ...
1143!> \param qmmm_spherical_cutoff ...
1144!> \param shells ...
1145!> \par History
1146!> 08.2004 created [tlaino]
1147!> \author Teodoro Laino
1148! **************************************************************************************************
1149 SUBROUTINE qmmm_forces_with_gaussian_lr(pgfs, cgrid, num_mm_atoms, mm_charges, mm_atom_index, &
1150 mm_particles, para_env, coarser_grid_level, Forces, potentials, &
1151 aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1152 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1153 TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
1154 INTEGER, INTENT(IN) :: num_mm_atoms
1155 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1156 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1157 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1158 TYPE(mp_para_env_type), POINTER :: para_env
1159 INTEGER, INTENT(IN) :: coarser_grid_level
1160 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
1161 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1162 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1163 TYPE(cell_type), POINTER :: mm_cell
1164 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1165 INTEGER, INTENT(IN) :: iw, par_scheme
1166 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1167 LOGICAL :: shells
1168
1169 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian_LR'
1170
1171 INTEGER :: handle, i, imm, indmm, iradtyp, ix, j, &
1172 k, lindmm, my_i, my_j, my_k, myind, &
1173 n1, n2, n3
1174 INTEGER, DIMENSION(2, 3) :: bo, gbo
1175 REAL(kind=dp) :: dr1, dr2, dr3, dvol, dx, fac, ft1, ft2, &
1176 ft3, qt, r, r2, rd1, rd2, rd3, rt1, &
1177 rt2, rt3, rv1, rv2, rv3, rx, rx2, &
1178 sph_chrg_factor, term, xs1, xs2, xs3
1179 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: lforces
1180 REAL(kind=dp), DIMENSION(3) :: ra
1181 REAL(kind=dp), DIMENSION(:, :), POINTER :: pot0_2
1182 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
1183 TYPE(qmmm_pot_type), POINTER :: pot
1184
1185 CALL timeset(routinen, handle)
1186 ALLOCATE (lforces(3, num_mm_atoms))
1187 lforces = 0.0_dp
1188 n1 = cgrid%pw_grid%npts(1)
1189 n2 = cgrid%pw_grid%npts(2)
1190 n3 = cgrid%pw_grid%npts(3)
1191 dr1 = cgrid%pw_grid%dr(1)
1192 dr2 = cgrid%pw_grid%dr(2)
1193 dr3 = cgrid%pw_grid%dr(3)
1194 dvol = cgrid%pw_grid%dvol
1195 gbo = cgrid%pw_grid%bounds
1196 bo = cgrid%pw_grid%bounds_local
1197 grid => cgrid%array
1198 IF (par_scheme == do_par_atom) myind = 0
1199 radius: DO iradtyp = 1, SIZE(pgfs)
1200 pot => potentials(iradtyp)%pot
1201 dx = pot%dx
1202 pot0_2 => pot%pot0_2
1203 !$OMP PARALLEL DO DEFAULT(NONE) &
1204 !$OMP SHARED(pot, par_scheme, para_env, dvol, mm_atom_index, mm_particles, dOmmOqm) &
1205 !$OMP SHARED(mm_cell, mm_charges, dx, LForces, Forces, qmmm_spherical_cutoff, shells, dr1, dr2, dr3, gbo, bo) &
1206 !$OMP SHARED(IRadTyp, pot0_2, grid) &
1207 !$OMP PRIVATE(Imm, myind, ra, LIndMM, IndMM, qt, rt1, rt2, rt3, ft1, ft2, ft3, i, j, k, sph_chrg_factor) &
1208 !$OMP PRIVATE(my_k, my_j, my_i, xs3, xs2, xs1, rv1, rv2, rv3, r, ix, rx, rx2, r2, Term, fac) &
1209 !$OMP PRIVATE(rd1, rd2, rd3)
1210 atoms: DO imm = 1, SIZE(pot%mm_atom_index)
1211 IF (par_scheme == do_par_atom) THEN
1212 myind = imm + (iradtyp - 1)*SIZE(pot%mm_atom_index)
1213 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle
1214 END IF
1215 lindmm = pot%mm_atom_index(imm)
1216 indmm = mm_atom_index(lindmm)
1217 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
1218 IF (shells) &
1219 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
1220 qt = mm_charges(lindmm)
1221 ! Possible Spherical Cutoff
1222 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1223 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
1224 qt = qt*sph_chrg_factor
1225 END IF
1226 IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
1227 rt1 = ra(1)
1228 rt2 = ra(2)
1229 rt3 = ra(3)
1230 ft1 = 0.0_dp
1231 ft2 = 0.0_dp
1232 ft3 = 0.0_dp
1233 loopongrid: DO k = bo(1, 3), bo(2, 3)
1234 my_k = k - gbo(1, 3)
1235 xs3 = real(my_k, dp)*dr3
1236 my_j = bo(1, 2) - gbo(1, 2)
1237 xs2 = real(my_j, dp)*dr2
1238 rv3 = rt3 - xs3
1239 DO j = bo(1, 2), bo(2, 2)
1240 my_i = bo(1, 1) - gbo(1, 1)
1241 xs1 = real(my_i, dp)*dr1
1242 rv2 = rt2 - xs2
1243 DO i = bo(1, 1), bo(2, 1)
1244 rv1 = rt1 - xs1
1245 r2 = rv1*rv1 + rv2*rv2 + rv3*rv3
1246 r = sqrt(r2)
1247 ix = floor(r/dx) + 1
1248 rx = (r - real(ix - 1, dp)*dx)/dx
1249 rx2 = rx*rx
1250 term = pot0_2(1, ix)*(-6.0_dp*(rx - rx2)) &
1251 + pot0_2(2, ix)*(1.0_dp - 4.0_dp*rx + 3.0_dp*rx2) &
1252 + pot0_2(1, ix + 1)*(6.0_dp*(rx - rx2)) &
1253 + pot0_2(2, ix + 1)*(-2.0_dp*rx + 3.0_dp*rx2)
1254 fac = grid(i, j, k)*term
1255 IF (r == 0.0_dp) THEN
1256 rd1 = 1.0_dp
1257 rd2 = 1.0_dp
1258 rd3 = 1.0_dp
1259 ELSE
1260 rd1 = rv1/r
1261 rd2 = rv2/r
1262 rd3 = rv3/r
1263 END IF
1264 ft1 = ft1 + fac*rd1
1265 ft2 = ft2 + fac*rd2
1266 ft3 = ft3 + fac*rd3
1267 xs1 = xs1 + dr1
1268 END DO
1269 xs2 = xs2 + dr2
1270 END DO
1271 END DO loopongrid
1272 qt = -qt*dvol/dx
1273 lforces(1, lindmm) = ft1*qt
1274 lforces(2, lindmm) = ft2*qt
1275 lforces(3, lindmm) = ft3*qt
1276
1277 forces(1, lindmm) = forces(1, lindmm) + lforces(1, lindmm)
1278 forces(2, lindmm) = forces(2, lindmm) + lforces(2, lindmm)
1279 forces(3, lindmm) = forces(3, lindmm) + lforces(3, lindmm)
1280 END DO atoms
1281 !$OMP END PARALLEL DO
1282 END DO radius
1283 !
1284 ! Debug Statement
1285 !
1286 IF (debug_this_module) THEN
1287 CALL debug_qmmm_forces_with_gauss_lr(pgfs=pgfs, &
1288 aug_pools=aug_pools, &
1289 rho=cgrid, &
1290 num_mm_atoms=num_mm_atoms, &
1291 mm_charges=mm_charges, &
1292 mm_atom_index=mm_atom_index, &
1293 mm_particles=mm_particles, &
1294 coarser_grid_level=coarser_grid_level, &
1295 debug_force=lforces, &
1296 potentials=potentials, &
1297 para_env=para_env, &
1298 mm_cell=mm_cell, &
1299 dommoqm=dommoqm, &
1300 iw=iw, &
1301 par_scheme=par_scheme, &
1302 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1303 shells=shells)
1304 END IF
1305
1306 DEALLOCATE (lforces)
1307 CALL timestop(handle)
1308 END SUBROUTINE qmmm_forces_with_gaussian_lr
1309
1310! **************************************************************************************************
1311!> \brief Evaluates numerically QM/MM forces and compares them with
1312!> the analytically computed ones.
1313!> It is evaluated only when debug_this_module is set to .TRUE.
1314!> \param rho ...
1315!> \param qs_env ...
1316!> \param qmmm_env ...
1317!> \param Analytical_Forces ...
1318!> \param mm_particles ...
1319!> \param mm_atom_index ...
1320!> \param num_mm_atoms ...
1321!> \param interp_section ...
1322!> \param mm_cell ...
1323!> \par History
1324!> 08.2004 created [tlaino]
1325!> \author Teodoro Laino
1326! **************************************************************************************************
1327 SUBROUTINE qmmm_debug_forces(rho, qs_env, qmmm_env, Analytical_Forces, &
1328 mm_particles, mm_atom_index, num_mm_atoms, &
1329 interp_section, mm_cell)
1330 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1331 TYPE(qs_environment_type), POINTER :: qs_env
1332 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1333 REAL(kind=dp), DIMENSION(:, :), POINTER :: analytical_forces
1334 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1335 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1336 INTEGER, INTENT(IN) :: num_mm_atoms
1337 TYPE(section_vals_type), POINTER :: interp_section
1338 TYPE(cell_type), POINTER :: mm_cell
1339
1340 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_debug_forces'
1341
1342 INTEGER :: handle, i, indmm, iw, j, k
1343 REAL(kind=dp) :: coord_save
1344 REAL(kind=dp), DIMENSION(2) :: energy
1345 REAL(kind=dp), DIMENSION(3) :: err
1346 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1347 TYPE(cp_logger_type), POINTER :: logger
1348 TYPE(mp_para_env_type), POINTER :: para_env
1349 TYPE(pw_env_type), POINTER :: pw_env
1350 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1351 TYPE(pw_r3d_rs_type) :: v_qmmm_rspace
1352 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
1353 TYPE(section_vals_type), POINTER :: input_section, print_section
1354
1355 CALL timeset(routinen, handle)
1356 NULLIFY (num_forces)
1357 CALL get_qs_env(qs_env=qs_env, &
1358 pw_env=pw_env, &
1359 input=input_section, &
1360 para_env=para_env)
1361
1362 print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
1363 logger => cp_get_default_logger()
1364 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".qmmmLog")
1365 CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools)
1366 CALL pw_pools(1)%pool%create_pw(v_qmmm_rspace)
1367 ALLOCATE (num_forces(3, num_mm_atoms))
1368 ks_qmmm_env_loc => qs_env%ks_qmmm_env
1369 IF (iw > 0) WRITE (iw, '(/A)') "DEBUG SECTION:"
1370 atoms: DO i = 1, num_mm_atoms
1371 indmm = mm_atom_index(i)
1372 coords: DO j = 1, 3
1373 coord_save = mm_particles(indmm)%r(j)
1374 energy = 0.0_dp
1375 diff: DO k = 1, 2
1376 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1377 CALL pw_zero(v_qmmm_rspace)
1378 SELECT CASE (qmmm_env%qmmm_coupl_type)
1379 CASE (do_qmmm_coulomb)
1380 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1381 CASE (do_qmmm_pcharge)
1382 cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1383 CASE (do_qmmm_gauss, do_qmmm_swave)
1384 CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env, &
1385 v_qmmm=v_qmmm_rspace, &
1386 mm_particles=mm_particles, &
1387 aug_pools=qmmm_env%aug_pools, &
1388 para_env=para_env, &
1389 eps_mm_rspace=qmmm_env%eps_mm_rspace, &
1390 cube_info=ks_qmmm_env_loc%cube_info, &
1391 pw_pools=pw_pools, &
1392 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
1393 coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
1394 interp_section=interp_section, &
1395 mm_cell=mm_cell)
1396 CASE (do_qmmm_none)
1397 cycle diff
1398 CASE DEFAULT
1399 IF (iw > 0) WRITE (iw, '(T3,A)') "Unknown Coupling..."
1400 cpabort("")
1401 END SELECT
1402 energy(k) = pw_integral_ab(rho, v_qmmm_rspace)
1403 END DO diff
1404 IF (iw > 0) THEN
1405 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1406 "DEBUG :: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1407 END IF
1408 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1409 mm_particles(indmm)%r(j) = coord_save
1410 END DO coords
1411 END DO atoms
1412
1413 SELECT CASE (qmmm_env%qmmm_coupl_type)
1414 CASE (do_qmmm_coulomb)
1415 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1416 CASE (do_qmmm_pcharge)
1417 cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1418 CASE (do_qmmm_gauss, do_qmmm_swave)
1419 IF (iw > 0) WRITE (iw, '(/A/)') "CHECKING NUMERICAL Vs ANALYTICAL FORCES (Err%):"
1420 DO i = 1, num_mm_atoms
1421 indmm = mm_atom_index(i)
1422 err = 0.0_dp
1423 DO k = 1, 3
1424 IF (abs(num_forces(k, i)) >= 5.0e-5_dp) THEN
1425 err(k) = (analytical_forces(k, i) - num_forces(k, i))/num_forces(k, i)*100.0_dp
1426 END IF
1427 END DO
1428 IF (iw > 0) &
1429 WRITE (iw, 100) indmm, analytical_forces(1, i), num_forces(1, i), err(1), &
1430 analytical_forces(2, i), num_forces(2, i), err(2), &
1431 analytical_forces(3, i), num_forces(3, i), err(3)
1432 cpassert(abs(err(1)) <= maxerr)
1433 cpassert(abs(err(2)) <= maxerr)
1434 cpassert(abs(err(3)) <= maxerr)
1435 END DO
1436 CASE (do_qmmm_none)
1437 IF (iw > 0) WRITE (iw, '(T3,A)') "No QM/MM Derivatives to debug. Just Mechanical Coupling!"
1438 CASE DEFAULT
1439 IF (iw > 0) WRITE (iw, '(T3,A)') "Unknown Coupling..."
1440 cpabort("")
1441 END SELECT
1442 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
1443
1444 CALL pw_pools(1)%pool%give_back_pw(v_qmmm_rspace)
1445 DEALLOCATE (num_forces)
1446 CALL timestop(handle)
1447100 FORMAT(i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1448 END SUBROUTINE qmmm_debug_forces
1449
1450! **************************************************************************************************
1451!> \brief Debugs the integrate_gf_rspace_NoPBC.. It may helps ;-P
1452!> \param ilevel ...
1453!> \param zetp ...
1454!> \param rp ...
1455!> \param W ...
1456!> \param pwgrid ...
1457!> \param cube_info ...
1458!> \param eps_mm_rspace ...
1459!> \param aug_pools ...
1460!> \param debug_force ...
1461!> \param mm_cell ...
1462!> \param auxbas_grid ...
1463!> \param n_rep_real ...
1464!> \param iw ...
1465!> \par History
1466!> 08.2004 created [tlaino]
1467!> \author Teodoro Laino
1468! **************************************************************************************************
1469 SUBROUTINE debug_integrate_gf_rspace_nopbc(ilevel, zetp, rp, W, pwgrid, cube_info, &
1470 eps_mm_rspace, aug_pools, debug_force, &
1471 mm_cell, auxbas_grid, n_rep_real, iw)
1472 INTEGER, INTENT(IN) :: ilevel
1473 REAL(kind=dp), INTENT(IN) :: zetp
1474 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rp
1475 REAL(kind=dp), INTENT(IN) :: w
1476 TYPE(pw_r3d_rs_type), INTENT(IN) :: pwgrid
1477 TYPE(cube_info_type), INTENT(IN) :: cube_info
1478 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
1479 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1480 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: debug_force
1481 TYPE(cell_type), POINTER :: mm_cell
1482 INTEGER, INTENT(IN) :: auxbas_grid
1483 INTEGER, DIMENSION(3), INTENT(IN) :: n_rep_real
1484 INTEGER, INTENT(IN) :: iw
1485
1486 CHARACTER(len=*), PARAMETER :: routinen = 'debug_integrate_gf_rspace_NoPBC'
1487
1488 INTEGER :: handle, i, igrid, k, ngrids
1489 INTEGER, DIMENSION(2, 3) :: bo2
1490 INTEGER, SAVE :: icount
1491 REAL(kind=dp), DIMENSION(2) :: energy
1492 REAL(kind=dp), DIMENSION(3) :: err, force, myrp
1493 REAL(kind=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat
1494 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1495
1496 DATA icount/0/
1497 ! Statements
1498 CALL timeset(routinen, handle)
1499 !Statements
1500 ngrids = SIZE(aug_pools)
1501 CALL pw_pools_create_pws(aug_pools, grids)
1502 DO igrid = 1, ngrids
1503 CALL pw_zero(grids(igrid))
1504 END DO
1505 bo2 = grids(auxbas_grid)%pw_grid%bounds
1506 ALLOCATE (xdat(bo2(1, 1):bo2(2, 1)))
1507 ALLOCATE (ydat(bo2(1, 2):bo2(2, 2)))
1508 ALLOCATE (zdat(bo2(1, 3):bo2(2, 3)))
1509
1510 icount = icount + 1
1511 DO i = 1, 3
1512 DO k = 1, 2
1513 myrp = rp
1514 myrp(i) = myrp(i) + (-1.0_dp)**k*dx
1515 CALL pw_zero(grids(ilevel))
1516 CALL collocate_gf_rspace_nopbc(zetp=zetp, &
1517 rp=myrp, &
1518 scale=-1.0_dp, &
1519 w=w, &
1520 pwgrid=grids(ilevel), &
1521 cube_info=cube_info, &
1522 eps_mm_rspace=eps_mm_rspace, &
1523 xdat=xdat, &
1524 ydat=ydat, &
1525 zdat=zdat, &
1526 bo2=bo2, &
1527 n_rep_real=n_rep_real, &
1528 mm_cell=mm_cell)
1529
1530 energy(k) = pw_integral_ab(pwgrid, grids(ilevel))
1531 END DO
1532 force(i) = (energy(2) - energy(1))/(2.0_dp*dx)
1533 END DO
1534 err = 0.0_dp
1535 IF (all(force /= 0.0_dp)) THEN
1536 err(1) = (debug_force(1) - force(1))/force(1)*100.0_dp
1537 err(2) = (debug_force(2) - force(2))/force(2)*100.0_dp
1538 err(3) = (debug_force(3) - force(3))/force(3)*100.0_dp
1539 END IF
1540 IF (iw > 0) &
1541 WRITE (iw, 100) icount, debug_force(1), force(1), err(1), &
1542 debug_force(2), force(2), err(2), &
1543 debug_force(3), force(3), err(3)
1544 cpassert(abs(err(1)) <= maxerr)
1545 cpassert(abs(err(2)) <= maxerr)
1546 cpassert(abs(err(3)) <= maxerr)
1547
1548 IF (ASSOCIATED(xdat)) THEN
1549 DEALLOCATE (xdat)
1550 END IF
1551 IF (ASSOCIATED(ydat)) THEN
1552 DEALLOCATE (ydat)
1553 END IF
1554 IF (ASSOCIATED(zdat)) THEN
1555 DEALLOCATE (zdat)
1556 END IF
1557
1558 CALL pw_pools_give_back_pws(aug_pools, grids)
1559 CALL timestop(handle)
1560100 FORMAT("Collocation : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1561 END SUBROUTINE debug_integrate_gf_rspace_nopbc
1562
1563! **************************************************************************************************
1564!> \brief Debugs qmmm_forces_with_gaussian_LG.. It may helps too ... ;-]
1565!> \param pgfs ...
1566!> \param aug_pools ...
1567!> \param rho ...
1568!> \param mm_charges ...
1569!> \param mm_atom_index ...
1570!> \param mm_particles ...
1571!> \param num_mm_atoms ...
1572!> \param coarser_grid_level ...
1573!> \param per_potentials ...
1574!> \param debug_force ...
1575!> \param para_env ...
1576!> \param mm_cell ...
1577!> \param dOmmOqm ...
1578!> \param iw ...
1579!> \param par_scheme ...
1580!> \param qmmm_spherical_cutoff ...
1581!> \param shells ...
1582!> \par History
1583!> 08.2004 created [tlaino]
1584!> \author Teodoro Laino
1585! **************************************************************************************************
1586 SUBROUTINE debug_qmmm_forces_with_gauss_lg(pgfs, aug_pools, rho, mm_charges, mm_atom_index, &
1587 mm_particles, num_mm_atoms, coarser_grid_level, per_potentials, &
1588 debug_force, para_env, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1589
1590 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1591 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1592 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1593 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1594 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1595 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1596 INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level
1597 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
1598 REAL(kind=dp), DIMENSION(:, :) :: debug_force
1599 TYPE(mp_para_env_type), POINTER :: para_env
1600 TYPE(cell_type), POINTER :: mm_cell
1601 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1602 INTEGER, INTENT(IN) :: iw, par_scheme
1603 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1604 LOGICAL :: shells
1605
1606 CHARACTER(len=*), PARAMETER :: routinen = 'debug_qmmm_forces_with_gauss_LG'
1607
1608 INTEGER :: handle, i, igrid, indmm, j, k, ngrids
1609 REAL(kind=dp) :: coord_save
1610 REAL(kind=dp), DIMENSION(2) :: energy
1611 REAL(kind=dp), DIMENSION(3) :: err
1612 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1613 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1614
1615 ALLOCATE (num_forces(3, num_mm_atoms))
1616 CALL timeset(routinen, handle)
1617 ngrids = SIZE(aug_pools)
1618 CALL pw_pools_create_pws(aug_pools, grids)
1619 DO igrid = 1, ngrids
1620 CALL pw_zero(grids(igrid))
1621 END DO
1622 atoms: DO i = 1, num_mm_atoms
1623 indmm = mm_atom_index(i)
1624 coords: DO j = 1, 3
1625 coord_save = mm_particles(indmm)%r(j)
1626 energy = 0.0_dp
1627 diff: DO k = 1, 2
1628 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1629 CALL pw_zero(grids(coarser_grid_level))
1630
1631 CALL qmmm_elec_with_gaussian_lg(pgfs=pgfs, &
1632 cgrid=grids(coarser_grid_level), &
1633 mm_charges=mm_charges, &
1634 mm_atom_index=mm_atom_index, &
1635 mm_particles=mm_particles, &
1636 para_env=para_env, &
1637 per_potentials=per_potentials, &
1638 mm_cell=mm_cell, &
1639 dommoqm=dommoqm, &
1640 par_scheme=par_scheme, &
1641 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1642 shells=shells)
1643
1644 energy(k) = pw_integral_ab(rho, grids(coarser_grid_level))
1645 END DO diff
1646 IF (iw > 0) &
1647 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1648 "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1649 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1650 mm_particles(indmm)%r(j) = coord_save
1651 END DO coords
1652 END DO atoms
1653
1654 DO i = 1, num_mm_atoms
1655 indmm = mm_atom_index(i)
1656 err = 0.0_dp
1657 IF (all(num_forces /= 0.0_dp)) THEN
1658 err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1659 err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1660 err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1661 END IF
1662 IF (iw > 0) &
1663 WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1664 debug_force(2, i), num_forces(2, i), err(2), &
1665 debug_force(3, i), num_forces(3, i), err(3)
1666 cpassert(abs(err(1)) <= maxerr)
1667 cpassert(abs(err(2)) <= maxerr)
1668 cpassert(abs(err(3)) <= maxerr)
1669 END DO
1670
1671 DEALLOCATE (num_forces)
1672 CALL pw_pools_give_back_pws(aug_pools, grids)
1673 CALL timestop(handle)
1674100 FORMAT("MM Atom LR : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1675 END SUBROUTINE debug_qmmm_forces_with_gauss_lg
1676
1677! **************************************************************************************************
1678!> \brief Debugs qmmm_forces_with_gaussian_LR.. It may helps too ... ;-]
1679!> \param pgfs ...
1680!> \param aug_pools ...
1681!> \param rho ...
1682!> \param mm_charges ...
1683!> \param mm_atom_index ...
1684!> \param mm_particles ...
1685!> \param num_mm_atoms ...
1686!> \param coarser_grid_level ...
1687!> \param potentials ...
1688!> \param debug_force ...
1689!> \param para_env ...
1690!> \param mm_cell ...
1691!> \param dOmmOqm ...
1692!> \param iw ...
1693!> \param par_scheme ...
1694!> \param qmmm_spherical_cutoff ...
1695!> \param shells ...
1696!> \par History
1697!> 08.2004 created [tlaino]
1698!> \author Teodoro Laino
1699! **************************************************************************************************
1700 SUBROUTINE debug_qmmm_forces_with_gauss_lr(pgfs, aug_pools, rho, mm_charges, mm_atom_index, &
1701 mm_particles, num_mm_atoms, coarser_grid_level, potentials, &
1702 debug_force, para_env, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1703
1704 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1705 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1706 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1707 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1708 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1709 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1710 INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level
1711 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1712 REAL(kind=dp), DIMENSION(:, :) :: debug_force
1713 TYPE(mp_para_env_type), POINTER :: para_env
1714 TYPE(cell_type), POINTER :: mm_cell
1715 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1716 INTEGER, INTENT(IN) :: iw, par_scheme
1717 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1718 LOGICAL :: shells
1719
1720 CHARACTER(len=*), PARAMETER :: routinen = 'debug_qmmm_forces_with_gauss_LR'
1721
1722 INTEGER :: handle, i, igrid, indmm, j, k, ngrids
1723 REAL(kind=dp) :: coord_save
1724 REAL(kind=dp), DIMENSION(2) :: energy
1725 REAL(kind=dp), DIMENSION(3) :: err
1726 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1727 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1728
1729 ALLOCATE (num_forces(3, num_mm_atoms))
1730 CALL timeset(routinen, handle)
1731 ngrids = SIZE(aug_pools)
1732 CALL pw_pools_create_pws(aug_pools, grids)
1733 DO igrid = 1, ngrids
1734 CALL pw_zero(grids(igrid))
1735 END DO
1736 atoms: DO i = 1, num_mm_atoms
1737 indmm = mm_atom_index(i)
1738 coords: DO j = 1, 3
1739 coord_save = mm_particles(indmm)%r(j)
1740 energy = 0.0_dp
1741 diff: DO k = 1, 2
1742 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1743 CALL pw_zero(grids(coarser_grid_level))
1744
1745 CALL qmmm_elec_with_gaussian_lr(pgfs=pgfs, &
1746 grid=grids(coarser_grid_level), &
1747 mm_charges=mm_charges, &
1748 mm_atom_index=mm_atom_index, &
1749 mm_particles=mm_particles, &
1750 para_env=para_env, &
1751 potentials=potentials, &
1752 mm_cell=mm_cell, &
1753 dommoqm=dommoqm, &
1754 par_scheme=par_scheme, &
1755 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1756 shells=shells)
1757
1758 energy(k) = pw_integral_ab(rho, grids(coarser_grid_level))
1759 END DO diff
1760 IF (iw > 0) &
1761 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1762 "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1763 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1764 mm_particles(indmm)%r(j) = coord_save
1765 END DO coords
1766 END DO atoms
1767
1768 DO i = 1, num_mm_atoms
1769 indmm = mm_atom_index(i)
1770 err = 0.0_dp
1771 IF (all(num_forces(:, i) /= 0.0_dp)) THEN
1772 err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1773 err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1774 err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1775 END IF
1776 IF (iw > 0) &
1777 WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1778 debug_force(2, i), num_forces(2, i), err(2), &
1779 debug_force(3, i), num_forces(3, i), err(3)
1780 cpassert(abs(err(1)) <= maxerr)
1781 cpassert(abs(err(2)) <= maxerr)
1782 cpassert(abs(err(3)) <= maxerr)
1783 END DO
1784
1785 DEALLOCATE (num_forces)
1786 CALL pw_pools_give_back_pws(aug_pools, grids)
1787 CALL timestop(handle)
1788100 FORMAT("MM Atom LR : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1789 END SUBROUTINE debug_qmmm_forces_with_gauss_lr
1790
1791END MODULE qmmm_gpw_forces
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:51
cell_types::pbc
Definition cell_types.F:92
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_integral_ab
Definition pw_methods.F:221
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
pw_pool_types::pw_pools_create_pws
Definition pw_pool_types.F:139
pw_pool_types::pw_pools_give_back_pws
Definition pw_pool_types.F:150
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_spline_utils
utils to manipulate splines on the regular grid of a pw
Definition cp_spline_utils.F:14
cp_spline_utils::spline3_nopbc_interp
integer, parameter, public spline3_nopbc_interp
Definition cp_spline_utils.F:38
cp_spline_utils::spline3_pbc_interp
integer, parameter, public spline3_pbc_interp
Definition cp_spline_utils.F:38
cp_spline_utils::pw_restrict_s3
subroutine, public pw_restrict_s3(pw_fine_in, pw_coarse_out, coarse_pool, param_section)
restricts the function from a fine grid to a coarse one
Definition cp_spline_utils.F:55
cube_utils
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
Definition cube_utils.F:18
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_par_atom
integer, parameter, public do_par_atom
Definition input_constants.F:287
input_constants::do_qmmm_none
integer, parameter, public do_qmmm_none
Definition input_constants.F:385
input_constants::do_qmmm_pcharge
integer, parameter, public do_qmmm_pcharge
Definition input_constants.F:385
input_constants::do_qmmm_coulomb
integer, parameter, public do_qmmm_coulomb
Definition input_constants.F:385
input_constants::do_qmmm_swave
integer, parameter, public do_qmmm_swave
Definition input_constants.F:385
input_constants::do_qmmm_gauss
integer, parameter, public do_qmmm_gauss
Definition input_constants.F:385
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:153
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mm_collocate_potential
Calculate the MM potential by collocating the primitive Gaussian functions (pgf)
Definition mm_collocate_potential.F:15
mm_collocate_potential::integrate_gf_rspace_nopbc
subroutine, public integrate_gf_rspace_nopbc(zetp, rp, scale, w, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo, force, n_rep_real, mm_cell)
Main driver to integrate gaussian functions on a grid function without using periodic boundary condit...
Definition mm_collocate_potential.F:296
mm_collocate_potential::collocate_gf_rspace_nopbc
subroutine, public collocate_gf_rspace_nopbc(zetp, rp, scale, w, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo2, n_rep_real, mm_cell)
Main driver to collocate gaussian functions on grid without using periodic boundary conditions (NoPBC...
Definition mm_collocate_potential.F:173
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qmmm_gaussian_types
Sets the typo for the gaussian treatment of the qm/mm interaction.
Definition qmmm_gaussian_types.F:14
qmmm_gpw_energy
A collection of methods to treat the QM/MM electrostatic coupling.
Definition qmmm_gpw_energy.F:14
qmmm_gpw_energy::qmmm_elec_with_gaussian
subroutine, public qmmm_elec_with_gaussian(qmmm_env, v_qmmm, mm_particles, aug_pools, cube_info, para_env, eps_mm_rspace, pw_pools, auxbas_grid, coarser_grid, interp_section, mm_cell)
Compute the QM/MM electrostatic Interaction collocating the gaussian Electrostatic Potential.
Definition qmmm_gpw_energy.F:243
qmmm_gpw_energy::qmmm_elec_with_gaussian_lr
subroutine, public qmmm_elec_with_gaussian_lr(pgfs, grid, mm_charges, mm_atom_index, mm_particles, para_env, potentials, mm_cell, dommoqm, par_scheme, qmmm_spherical_cutoff, shells)
Compute the QM/MM electrostatic Interaction collocating (1/R - Sum_NG Gaussians) on the coarser grid ...
Definition qmmm_gpw_energy.F:757
qmmm_gpw_energy::qmmm_elec_with_gaussian_lg
subroutine, public qmmm_elec_with_gaussian_lg(pgfs, cgrid, mm_charges, mm_atom_index, mm_particles, para_env, per_potentials, mm_cell, dommoqm, par_scheme, qmmm_spherical_cutoff, shells)
Compute the QM/MM electrostatic Interaction collocating (1/R - Sum_NG Gaussians) on the coarser grid ...
Definition qmmm_gpw_energy.F:532
qmmm_gpw_forces
Routines to compute energy and forces in a QM/MM calculation.
Definition qmmm_gpw_forces.F:14
qmmm_gpw_forces::qmmm_forces
subroutine, public qmmm_forces(qs_env, qmmm_env, mm_particles, calc_force, mm_cell)
General driver to Compute the contribution to the forces due to the QM/MM potential.
Definition qmmm_gpw_forces.F:100
qmmm_se_forces
Calculation of the derivative of the QMMM Hamiltonian integral matrix <a|\sum_i q_i|b> for semi-empir...
Definition qmmm_se_forces.F:13
qmmm_se_forces::deriv_se_qmmm_matrix
subroutine, public deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian QMMM terms.
Definition qmmm_se_forces.F:78
qmmm_tb_methods
TB methods used with QMMM.
Definition qmmm_tb_methods.F:12
qmmm_tb_methods::deriv_tb_qmmm_matrix
subroutine, public deriv_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition qmmm_tb_methods.F:593
qmmm_tb_methods::deriv_tb_qmmm_matrix_pc
subroutine, public deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition qmmm_tb_methods.F:830
qmmm_types_low
Definition qmmm_types_low.F:13
qmmm_util
Definition qmmm_util.F:13
qmmm_util::spherical_cutoff_factor
subroutine, public spherical_cutoff_factor(spherical_cutoff, rij, factor)
Computes a spherical cutoff factor for the QMMM interactions.
Definition qmmm_util.F:615
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:526
qs_ks_qmmm_types
Definition qs_ks_qmmm_types.F:13
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:598
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cube_utils::cube_info_type
Definition cube_utils.F:38
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qmmm_gaussian_types::qmmm_gaussian_p_type
represent a pointer to a qmmm_gaussian_type, to be able to create arrays of pointers
Definition qmmm_gaussian_types.F:42
qmmm_gaussian_types::qmmm_gaussian_type
Definition qmmm_gaussian_types.F:27
qmmm_types_low::qmmm_env_qm_type
...
Definition qmmm_types_low.F:195
qmmm_types_low::qmmm_per_pot_p_type
Definition qmmm_types_low.F:100
qmmm_types_low::qmmm_per_pot_type
Periodic Potential.
Definition qmmm_types_low.F:90
qmmm_types_low::qmmm_pot_p_type
Definition qmmm_types_low.F:83
qmmm_types_low::qmmm_pot_type
Real Space Potential.
Definition qmmm_types_low.F:76
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_ks_qmmm_types::qs_ks_qmmm_env_type
calculation environment to calculate the ks_qmmm matrix, holds the QM/MM potential and all the needed...
Definition qs_ks_qmmm_types.F:47
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62