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