(git:71c3ab0)
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_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 cpabort("Unknown QM/MM coupling")
178 END SELECT
179 ELSE IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
180 ! DFTB
181 SELECT CASE (qmmm_env%qmmm_coupl_type)
182 CASE (do_qmmm_none)
183 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
184 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
185 CASE (do_qmmm_coulomb)
186 CALL deriv_tb_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
187 need_f, forces, forces_added_charges)
188 CASE (do_qmmm_pcharge)
189 CALL deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
190 need_f, forces, forces_added_charges)
191 CASE (do_qmmm_gauss, do_qmmm_swave)
192 cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for DFTB.")
193 CASE DEFAULT
194 cpabort("Unknown QM/MM coupling")
195 END SELECT
196 IF (need_f) THEN
197 forces(:, :) = forces(:, :)/real(para_env%num_pe, kind=dp)
198 forces_added_charges(:, :) = forces_added_charges(:, :)/real(para_env%num_pe, kind=dp)
199 END IF
200 ELSE
201 ! GPW/GAPW
202 CALL pw_env_get(pw_env=pw_env, &
203 pw_pools=pw_pools, &
204 auxbas_pw_pool=auxbas_pool)
205 CALL auxbas_pool%create_pw(rho_tot_r)
206 ! IF GAPW the core charge is replaced by the compensation charge
207 IF (gapw) THEN
208 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
209 CALL pw_transfer(rho_core, rho_tot_r)
210 energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
211 CALL auxbas_pool%create_pw(rho_tot_r2)
212 CALL pw_transfer(rho0_s_gs, rho_tot_r2)
213 CALL pw_axpy(rho_tot_r2, rho_tot_r)
214 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
215 CALL auxbas_pool%create_pw(rho_tot_r3)
216 CALL pw_transfer(rhoz_cneo_s_gs, rho_tot_r3)
217 CALL pw_axpy(rho_tot_r3, rho_tot_r)
218 CALL auxbas_pool%give_back_pw(rho_tot_r3)
219 END IF
220 CALL auxbas_pool%give_back_pw(rho_tot_r2)
221 ELSE
222 CALL pw_transfer(rho0_s_gs, rho_tot_r)
223 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
224 CALL auxbas_pool%create_pw(rho_tot_r3)
225 CALL pw_transfer(rhoz_cneo_s_gs, rho_tot_r3)
226 CALL pw_axpy(rho_tot_r3, rho_tot_r)
227 CALL auxbas_pool%give_back_pw(rho_tot_r3)
228 END IF
229 !
230 ! QM/MM Nuclear Electrostatic Potential already included through rho0
231 !
232 energy%qmmm_nu = 0.0_dp
233 END IF
234 ELSE
235 CALL pw_transfer(rho_core, rho_tot_r)
236 !
237 ! Computes the QM/MM Nuclear Electrostatic Potential
238 !
239 energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
240 END IF
241 IF (need_f) THEN
242 !
243 DO ispin = 1, SIZE(rho_r)
244 CALL pw_axpy(rho_r(ispin), rho_tot_r)
245 END DO
246 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Evaluating forces on MM atoms due to the:"
247 ! Electrostatic Interaction type...
248 SELECT CASE (qmmm_env%qmmm_coupl_type)
249 CASE (do_qmmm_coulomb)
250 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
251 CASE (do_qmmm_pcharge)
252 cpabort("Point Charge QM/MM electrostatic coupling not yet implemented for GPW/GAPW.")
253 CASE (do_qmmm_gauss, do_qmmm_swave)
254 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
255 "- QM/MM Coupling computed collocating the Gaussian Potential Functions."
256 CALL qmmm_forces_with_gaussian(rho=rho_tot_r, &
257 qmmm_env=qmmm_env, &
258 mm_particles=mm_particles, &
259 aug_pools=qmmm_env%aug_pools, &
260 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
261 coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
262 para_env=para_env, &
263 pw_pools=pw_pools, &
264 eps_mm_rspace=qmmm_env%eps_mm_rspace, &
265 cube_info=ks_qmmm_env_loc%cube_info, &
266 forces=forces, &
267 forces_added_charges=forces_added_charges, &
268 forces_added_shells=forces_added_shells, &
269 interp_section=interp_section, &
270 iw=iw, &
271 mm_cell=mm_cell)
272 CASE (do_qmmm_none)
273 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
274 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
275 CASE DEFAULT
276 cpabort("Unknown QM/MM coupling")
277 END SELECT
278 END IF
279 END IF
280 ! Correct Total Energy adding the contribution of the QM/MM nuclear interaction
281 energy%total = energy%total + energy%qmmm_nu
282 ! Proceed if gradients are requested..
283 IF (need_f) THEN
284 !ikuo Temporary change to alleviate compiler problems on Intel with
285 !array dimension of 0
286 IF (qmmm_env%num_mm_atoms /= 0) CALL para_env%sum(forces)
287 IF (qmmm_env%added_charges%num_mm_atoms /= 0) CALL para_env%sum(forces_added_charges)
288 IF (qmmm_env%added_shells%num_mm_atoms /= 0) CALL para_env%sum(forces_added_shells)
289 ! Debug Forces
290 IF (debug_this_module) THEN
291 IF (dft_control%qs_control%semi_empirical .OR. &
292 dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
293 WRITE (iw, *) "NO DEBUG AVAILABLE in module"//trim(routinen)
294 ELSE
295 ! Print Out Forces
296 IF (iw > 0) THEN
297 DO imm = 1, SIZE(qmmm_env%mm_atom_index)
298 WRITE (iw, *) "ANALYTICAL FORCES:"
299 indmm = qmmm_env%mm_atom_index(imm)
300 WRITE (iw, '(I6,3F15.9)') indmm, forces(:, imm)
301 END DO
302 END IF
303 CALL qmmm_debug_forces(rho=rho_tot_r, &
304 qs_env=qs_env, &
305 qmmm_env=qmmm_env, &
306 analytical_forces=forces, &
307 mm_particles=mm_particles, &
308 mm_atom_index=qmmm_env%mm_atom_index, &
309 num_mm_atoms=qmmm_env%num_mm_atoms, &
310 interp_section=interp_section, &
311 mm_cell=mm_cell)
312 END IF
313 END IF
314 END IF
315 ! Give back rho_tot_t to auxbas_pool only for GPW/GAPW
316 IF ((.NOT. dft_control%qs_control%semi_empirical) .AND. &
317 (.NOT. dft_control%qs_control%dftb) .AND. (.NOT. dft_control%qs_control%xtb)) THEN
318 CALL auxbas_pool%give_back_pw(rho_tot_r)
319 END IF
320 IF (iw > 0) THEN
321 IF (.NOT. gapw) WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
322 "QM/MM Nuclear Electrostatic Potential :", energy%qmmm_nu
323 WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
324 "QMMM Total Energy (QM + QMMM electronic + QMMM nuclear):", energy%total
325 WRITE (iw, '(T2,"QMMM|",1X,A)') "MM energy NOT included in the above term!"// &
326 " Check for: FORCE_EVAL ( QMMM )"
327 WRITE (iw, '(T2,"QMMM|",1X,A)') "that includes both QM, QMMM and MM energy terms!"
328 END IF
329 IF (need_f) THEN
330 ! Transfer Forces
331 DO imm = 1, qmmm_env%num_mm_atoms
332 indmm = qmmm_env%mm_atom_index(imm)
333
334 !add image forces to Forces
335 IF (qmmm_env%image_charge) THEN
336 DO iatom = 1, qmmm_env%num_image_mm_atoms
337 image_indmm = qmmm_env%image_charge_pot%image_mm_list(iatom)
338 IF (image_indmm == indmm) THEN
339 forces(:, imm) = forces(:, imm) &
340 + qmmm_env%image_charge_pot%image_forcesMM(:, iatom)
341 END IF
342 END DO
343 END IF
344
345 ! Hack: In Forces there the gradients indeed...
346 ! Minux sign to take care of this misunderstanding...
347 mm_particles(indmm)%f(:) = -forces(:, imm) + mm_particles(indmm)%f(:)
348 END DO
349 DEALLOCATE (forces)
350 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
351 DO imm = 1, qmmm_env%added_charges%num_mm_atoms
352 indmm = qmmm_env%added_charges%mm_atom_index(imm)
353 ! Hack: In Forces there the gradients indeed...
354 ! Minux sign to take care of this misunderstanding...
355 qmmm_env%added_charges%added_particles(indmm)%f(:) = -forces_added_charges(:, imm)
356 END DO
357 END IF
358 DEALLOCATE (forces_added_charges)
359 IF (qmmm_env%added_shells%num_mm_atoms > 0) THEN
360 DO imm = 1, qmmm_env%added_shells%num_mm_atoms
361 indmm = qmmm_env%added_shells%mm_core_index(imm)
362 ! Hack: In Forces there the gradients indeed...
363 ! Minux sign to take care of this misunderstanding...
364 qmmm_env%added_shells%added_particles(imm)%f(:) = qmmm_env%added_shells%added_particles(imm)%f(:) - &
365 forces_added_shells(:, imm)
366
367 END DO
368 END IF
369 DEALLOCATE (forces_added_shells)
370 END IF
371 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
372 CALL timestop(handle)
373
374 END SUBROUTINE qmmm_forces
375
376! **************************************************************************************************
377!> \brief Evaluates the contribution to the forces due to the
378!> QM/MM potential computed collocating the Electrostatic
379!> Gaussian Potential.
380!> \param rho ...
381!> \param qmmm_env ...
382!> \param mm_particles ...
383!> \param aug_pools ...
384!> \param auxbas_grid ...
385!> \param coarser_grid ...
386!> \param cube_info ...
387!> \param para_env ...
388!> \param eps_mm_rspace ...
389!> \param pw_pools ...
390!> \param Forces ...
391!> \param Forces_added_charges ...
392!> \param Forces_added_shells ...
393!> \param interp_section ...
394!> \param iw ...
395!> \param mm_cell ...
396!> \par History
397!> 06.2004 created [tlaino]
398!> \author Teodoro Laino
399! **************************************************************************************************
400 SUBROUTINE qmmm_forces_with_gaussian(rho, qmmm_env, mm_particles, &
401 aug_pools, auxbas_grid, coarser_grid, cube_info, para_env, &
402 eps_mm_rspace, pw_pools, Forces, Forces_added_charges, Forces_added_shells, &
403 interp_section, iw, mm_cell)
404 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
405 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
406 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
407 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
408 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
409 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
410 TYPE(mp_para_env_type), POINTER :: para_env
411 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
412 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
413 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges, &
414 forces_added_shells
415 TYPE(section_vals_type), POINTER :: interp_section
416 INTEGER, INTENT(IN) :: iw
417 TYPE(cell_type), POINTER :: mm_cell
418
419 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian'
420
421 INTEGER :: handle, i, igrid, j, k, kind_interp, me, &
422 ngrids
423 INTEGER, DIMENSION(3) :: glb, gub, lb, ub
424 INTEGER, DIMENSION(:), POINTER :: pos_of_x
425 LOGICAL :: shells
426 REAL(kind=dp), DIMENSION(:, :), POINTER :: tmp
427 TYPE(mp_comm_type) :: group
428 TYPE(mp_request_type) :: request
429 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
430
431! Statements
432
433 CALL timeset(routinen, handle)
434 NULLIFY (tmp)
435 cpassert(ASSOCIATED(mm_particles))
436 cpassert(ASSOCIATED(qmmm_env%mm_atom_chrg))
437 cpassert(ASSOCIATED(qmmm_env%mm_atom_index))
438 cpassert(ASSOCIATED(forces))
439 !Statements
440 ngrids = SIZE(pw_pools)
441 CALL pw_pools_create_pws(aug_pools, grids)
442 DO igrid = 1, ngrids
443 CALL pw_zero(grids(igrid))
444 END DO
445 ! Collocate Density on multigrids
446 lb = rho%pw_grid%bounds_local(1, :)
447 ub = rho%pw_grid%bounds_local(2, :)
448 grids(auxbas_grid)%array(lb(1):ub(1), &
449 lb(2):ub(2), &
450 lb(3):ub(3)) = rho%array
451 ! copy the boundaries
452 DO i = lb(1), ub(1)
453 grids(auxbas_grid)%array(i, ub(2) + 1, ub(3) + 1) = rho%array(i, lb(2), lb(3))
454 END DO
455 DO k = lb(3), ub(3)
456 DO i = lb(1), ub(1)
457 grids(auxbas_grid)%array(i, ub(2) + 1, k) = rho%array(i, lb(2), k)
458 END DO
459 END DO
460 DO j = lb(2), ub(2)
461 DO i = lb(1), ub(1)
462 grids(auxbas_grid)%array(i, j, ub(3) + 1) = rho%array(i, j, lb(3))
463 END DO
464 END DO
465 pos_of_x => grids(auxbas_grid)%pw_grid%para%pos_of_x
466 group = grids(auxbas_grid)%pw_grid%para%group
467 me = grids(auxbas_grid)%pw_grid%para%group%mepos
468 glb = rho%pw_grid%bounds(1, :)
469 gub = rho%pw_grid%bounds(2, :)
470 IF ((pos_of_x(glb(1)) == me) .AND. (pos_of_x(gub(1)) == me)) THEN
471 DO k = lb(3), ub(3)
472 DO j = lb(2), ub(2)
473 grids(auxbas_grid)%array(ub(1) + 1, j, k) = rho%array(lb(1), j, k)
474 END DO
475 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, k) = rho%array(lb(1), lb(2), k)
476 END DO
477 DO j = lb(2), ub(2)
478 grids(auxbas_grid)%array(ub(1) + 1, j, ub(3) + 1) = rho%array(lb(1), j, lb(3))
479 END DO
480 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, ub(3) + 1) = rho%array(lb(1), lb(2), lb(3))
481 ELSE IF (pos_of_x(glb(1)) == me) THEN
482 ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
483 rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)))
484 tmp = rho%array(lb(1), :, :)
485 CALL group%isend(msgin=tmp, dest=pos_of_x(rho%pw_grid%bounds(2, 1)), &
486 request=request, tag=112)
487 CALL request%wait()
488 ELSE IF (pos_of_x(gub(1)) == me) THEN
489 ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
490 rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)))
491 CALL group%irecv(msgout=tmp, source=pos_of_x(rho%pw_grid%bounds(1, 1)), &
492 request=request, tag=112)
493 CALL request%wait()
494
495 DO k = lb(3), ub(3)
496 DO j = lb(2), ub(2)
497 grids(auxbas_grid)%array(ub(1) + 1, j, k) = tmp(j, k)
498 END DO
499 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, k) = tmp(lb(2), k)
500 END DO
501 DO j = lb(2), ub(2)
502 grids(auxbas_grid)%array(ub(1) + 1, j, ub(3) + 1) = tmp(j, lb(3))
503 END DO
504 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, ub(3) + 1) = tmp(lb(2), lb(3))
505 END IF
506 IF (ASSOCIATED(tmp)) THEN
507 DEALLOCATE (tmp)
508 END IF
509 ! Further setup of parallelization scheme
510 IF (qmmm_env%par_scheme == do_par_atom) THEN
511 CALL para_env%sum(grids(auxbas_grid)%array)
512 END IF
513 ! RealSpace Interpolation
514 CALL section_vals_val_get(interp_section, "kind", i_val=kind_interp)
515 SELECT CASE (kind_interp)
516 CASE (spline3_nopbc_interp, spline3_pbc_interp)
517 ! Spline Interpolator
518 DO igrid = auxbas_grid, SIZE(grids) - 1
519 CALL pw_restrict_s3(grids(igrid), &
520 grids(igrid + 1), &
521 aug_pools(igrid + 1)%pool, &
522 param_section=interp_section)
523 END DO
524 CASE DEFAULT
525 cpabort("Unknown interpolation method")
526 END SELECT
527
528 shells = .false.
529 CALL qmmm_force_with_gaussian_low(grids, mm_particles, &
530 qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
531 qmmm_env%num_mm_atoms, cube_info, para_env, eps_mm_rspace, auxbas_grid, &
532 coarser_grid, qmmm_env%pgfs, qmmm_env%potentials, forces, aug_pools, &
533 mm_cell, qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%per_potentials, &
534 iw, qmmm_env%par_scheme, qmmm_env%spherical_cutoff, shells)
535
536 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
537 CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_charges%added_particles, &
538 qmmm_env%added_charges%mm_atom_chrg, &
539 qmmm_env%added_charges%mm_atom_index, qmmm_env%added_charges%num_mm_atoms, &
540 cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_charges%pgfs, &
541 qmmm_env%added_charges%potentials, forces_added_charges, aug_pools, mm_cell, &
542 qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%added_charges%per_potentials, iw, qmmm_env%par_scheme, &
543 qmmm_env%spherical_cutoff, shells)
544 END IF
545
546 IF (qmmm_env%added_shells%num_mm_atoms > 0) THEN
547 shells = .true.
548 CALL qmmm_force_with_gaussian_low(grids, qmmm_env%added_shells%added_particles, &
549 qmmm_env%added_shells%mm_core_chrg, &
550 qmmm_env%added_shells%mm_core_index, qmmm_env%added_shells%num_mm_atoms, &
551 cube_info, para_env, eps_mm_rspace, auxbas_grid, coarser_grid, qmmm_env%added_shells%pgfs, &
552 qmmm_env%added_shells%potentials, forces_added_shells, aug_pools, mm_cell, &
553 qmmm_env%dOmmOqm, qmmm_env%periodic, qmmm_env%added_shells%per_potentials, iw, qmmm_env%par_scheme, &
554 qmmm_env%spherical_cutoff, shells)
555 END IF
556
557 CALL pw_pools_give_back_pws(aug_pools, grids)
558 CALL timestop(handle)
559
560 END SUBROUTINE qmmm_forces_with_gaussian
561
562! **************************************************************************************************
563!> \brief Evaluates the contribution to the forces due to the
564!> QM/MM potential computed collocating the Electrostatic
565!> Gaussian Potential. Low Level
566!> \param grids ...
567!> \param mm_particles ...
568!> \param mm_charges ...
569!> \param mm_atom_index ...
570!> \param num_mm_atoms ...
571!> \param cube_info ...
572!> \param para_env ...
573!> \param eps_mm_rspace ...
574!> \param auxbas_grid ...
575!> \param coarser_grid ...
576!> \param pgfs ...
577!> \param potentials ...
578!> \param Forces ...
579!> \param aug_pools ...
580!> \param mm_cell ...
581!> \param dOmmOqm ...
582!> \param periodic ...
583!> \param per_potentials ...
584!> \param iw ...
585!> \param par_scheme ...
586!> \param qmmm_spherical_cutoff ...
587!> \param shells ...
588!> \par History
589!> 06.2004 created [tlaino]
590!> \author Teodoro Laino
591! **************************************************************************************************
592 SUBROUTINE qmmm_force_with_gaussian_low(grids, mm_particles, mm_charges, &
593 mm_atom_index, num_mm_atoms, cube_info, para_env, &
594 eps_mm_rspace, auxbas_grid, coarser_grid, pgfs, potentials, Forces, &
595 aug_pools, mm_cell, dOmmOqm, periodic, per_potentials, iw, par_scheme, &
596 qmmm_spherical_cutoff, shells)
597 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: grids
598 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
599 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
600 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
601 INTEGER, INTENT(IN) :: num_mm_atoms
602 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
603 TYPE(mp_para_env_type), POINTER :: para_env
604 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
605 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
606 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
607 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
608 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
609 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
610 TYPE(cell_type), POINTER :: mm_cell
611 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
612 LOGICAL, INTENT(in) :: periodic
613 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
614 INTEGER, INTENT(IN) :: iw, par_scheme
615 REAL(kind=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
616 LOGICAL, INTENT(in) :: shells
617
618 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_force_with_gaussian_low', &
619 routinenb = 'qmmm_forces_gaussian_low'
620
621 INTEGER :: handle, handle2, igauss, ilevel, imm, &
622 indmm, iradtyp, lindmm, myind, &
623 n_rep_real(3)
624 INTEGER, DIMENSION(2, 3) :: bo
625 REAL(kind=dp) :: alpha, dvol, height, sph_chrg_factor, w
626 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xdat, ydat, zdat
627 REAL(kind=dp), DIMENSION(3) :: force, ra
628 TYPE(qmmm_gaussian_type), POINTER :: pgf
629 TYPE(qmmm_per_pot_type), POINTER :: per_pot
630 TYPE(qmmm_pot_type), POINTER :: pot
631
632 CALL timeset(routinen, handle)
633 CALL timeset(routinenb//"_G", handle2)
634 NULLIFY (pgf, pot, per_pot)
635 IF (par_scheme == do_par_atom) myind = 0
636 radius: DO iradtyp = 1, SIZE(pgfs)
637 pgf => pgfs(iradtyp)%pgf
638 pot => potentials(iradtyp)%pot
639 n_rep_real = 0
640 IF (periodic) THEN
641 per_pot => per_potentials(iradtyp)%pot
642 n_rep_real = per_pot%n_rep_real
643 END IF
644 gaussian: DO igauss = 1, pgf%Number_of_Gaussians
645 alpha = 1.0_dp/pgf%Gk(igauss)
646 alpha = alpha*alpha
647 height = pgf%Ak(igauss)
648 ilevel = pgf%grid_level(igauss)
649 dvol = grids(ilevel)%pw_grid%dvol
650 bo = grids(ilevel)%pw_grid%bounds_local
651 ALLOCATE (xdat(2, bo(1, 1):bo(2, 1)))
652 ALLOCATE (ydat(2, bo(1, 2):bo(2, 2)))
653 ALLOCATE (zdat(2, bo(1, 3):bo(2, 3)))
654 !$OMP PARALLEL DO DEFAULT(NONE) &
655 !$OMP SHARED(pot, par_scheme, dvol, alpha, para_env, mm_atom_index, shells) &
656 !$OMP SHARED(mm_particles, dOmmOqm, mm_cell, height, mm_charges, qmmm_spherical_cutoff) &
657 !$OMP SHARED(grids, cube_info, bo, n_rep_real, eps_mm_rspace, Forces, ilevel) &
658 !$OMP SHARED(IGauss, pgf, IRadTyp, iw, aug_pools, auxbas_grid) &
659 !$OMP PRIVATE(xdat, ydat, zdat) &
660 !$OMP PRIVATE(Imm, LIndMM, IndMM, ra, W, force, sph_chrg_factor, myind)
661 atoms: DO imm = 1, SIZE(pot%mm_atom_index)
662 IF (par_scheme == do_par_atom) THEN
663 myind = imm + (igauss - 1)*SIZE(pot%mm_atom_index) + (iradtyp - 1)*pgf%Number_of_Gaussians
664 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
665 END IF
666 lindmm = pot%mm_atom_index(imm)
667 indmm = mm_atom_index(lindmm)
668 IF (shells) THEN
669 ra(:) = pbc(mm_particles(imm)%r - dommoqm, mm_cell) + dommoqm
670 ELSE
671 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
672 END IF
673 w = mm_charges(lindmm)*height
674 force = 0.0_dp
675 ! Possible Spherical Cutoff
676 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
677 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
678 w = w*sph_chrg_factor
679 END IF
680 IF (abs(w) <= epsilon(0.0_dp)) cycle atoms
681 CALL integrate_gf_rspace_nopbc(zetp=alpha, &
682 rp=ra, &
683 scale=-1.0_dp, &
684 w=w, &
685 pwgrid=grids(ilevel), &
686 cube_info=cube_info(ilevel), &
687 eps_mm_rspace=eps_mm_rspace, &
688 xdat=xdat, &
689 ydat=ydat, &
690 zdat=zdat, &
691 bo=bo, &
692 force=force, &
693 n_rep_real=n_rep_real, &
694 mm_cell=mm_cell)
695 force = force*dvol
696 forces(:, lindmm) = forces(:, lindmm) + force(:)
697 !
698 ! Debug Statement
699 !
700 IF (debug_this_module) THEN
701 CALL debug_integrate_gf_rspace_nopbc(ilevel=ilevel, &
702 zetp=alpha, &
703 rp=ra, &
704 w=w, &
705 pwgrid=grids(ilevel), &
706 cube_info=cube_info(ilevel), &
707 eps_mm_rspace=eps_mm_rspace, &
708 aug_pools=aug_pools, &
709 debug_force=force, &
710 mm_cell=mm_cell, &
711 auxbas_grid=auxbas_grid, &
712 n_rep_real=n_rep_real, &
713 iw=iw)
714 END IF
715 END DO atoms
716 !$OMP END PARALLEL DO
717 DEALLOCATE (xdat)
718 DEALLOCATE (ydat)
719 DEALLOCATE (zdat)
720 END DO gaussian
721 END DO radius
722 CALL timestop(handle2)
723 CALL timeset(routinenb//"_R", handle2)
724 IF (periodic) THEN
725 CALL qmmm_forces_with_gaussian_lg(pgfs=pgfs, &
726 cgrid=grids(coarser_grid), &
727 num_mm_atoms=num_mm_atoms, &
728 mm_charges=mm_charges, &
729 mm_atom_index=mm_atom_index, &
730 mm_particles=mm_particles, &
731 para_env=para_env, &
732 coarser_grid_level=coarser_grid, &
733 forces=forces, &
734 per_potentials=per_potentials, &
735 aug_pools=aug_pools, &
736 mm_cell=mm_cell, &
737 dommoqm=dommoqm, &
738 iw=iw, &
739 par_scheme=par_scheme, &
740 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
741 shells=shells)
742 ELSE
743 CALL qmmm_forces_with_gaussian_lr(pgfs=pgfs, &
744 cgrid=grids(coarser_grid), &
745 num_mm_atoms=num_mm_atoms, &
746 mm_charges=mm_charges, &
747 mm_atom_index=mm_atom_index, &
748 mm_particles=mm_particles, &
749 para_env=para_env, &
750 coarser_grid_level=coarser_grid, &
751 forces=forces, &
752 potentials=potentials, &
753 aug_pools=aug_pools, &
754 mm_cell=mm_cell, &
755 dommoqm=dommoqm, &
756 iw=iw, &
757 par_scheme=par_scheme, &
758 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
759 shells=shells)
760 END IF
761 CALL timestop(handle2)
762 CALL timestop(handle)
763 END SUBROUTINE qmmm_force_with_gaussian_low
764
765! **************************************************************************************************
766!> \brief Evaluates the contribution to the forces due to the Long Range
767!> part of the QM/MM potential computed collocating the Electrostatic
768!> Gaussian Potential.
769!> \param pgfs ...
770!> \param cgrid ...
771!> \param num_mm_atoms ...
772!> \param mm_charges ...
773!> \param mm_atom_index ...
774!> \param mm_particles ...
775!> \param para_env ...
776!> \param coarser_grid_level ...
777!> \param Forces ...
778!> \param per_potentials ...
779!> \param aug_pools ...
780!> \param mm_cell ...
781!> \param dOmmOqm ...
782!> \param iw ...
783!> \param par_scheme ...
784!> \param qmmm_spherical_cutoff ...
785!> \param shells ...
786!> \par History
787!> 08.2004 created [tlaino]
788!> \author Teodoro Laino
789! **************************************************************************************************
790 SUBROUTINE qmmm_forces_with_gaussian_lg(pgfs, cgrid, num_mm_atoms, mm_charges, mm_atom_index, &
791 mm_particles, para_env, coarser_grid_level, Forces, per_potentials, &
792 aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
793 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
794 TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
795 INTEGER, INTENT(IN) :: num_mm_atoms
796 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
797 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
798 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
799 TYPE(mp_para_env_type), POINTER :: para_env
800 INTEGER, INTENT(IN) :: coarser_grid_level
801 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
802 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
803 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
804 TYPE(cell_type), POINTER :: mm_cell
805 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
806 INTEGER, INTENT(IN) :: iw, par_scheme
807 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
808 LOGICAL :: shells
809
810 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian_LG'
811
812 INTEGER :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, ij3, ij4, ik1, ik2, ik3, ik4, imm, &
813 indmm, iradtyp, ivec(3), j, k, lindmm, my_i, my_j, my_k, myind, npts(3)
814 INTEGER, DIMENSION(2, 3) :: bo, gbo
815 REAL(kind=dp) :: a1, a2, a3, abc_x(4, 4), abc_x_y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, &
816 dr1, dr1c, dr1i, dr2, dr2c, dr2i, dr3, dr3c, dr3i, dvol, e1, e2, e3, f1, f2, f3, fac, &
817 ft1, ft2, ft3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, &
818 rt3, rv1, rv2, rv3, s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, &
819 sph_chrg_factor, t1, t1d, t1o, t2, t2d, t2o, t3, t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, &
820 v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, v4o, xd1, xd2, xd3, xs1, xs2, xs3
821 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: lforces
822 REAL(kind=dp), DIMENSION(3) :: ra, val, vec
823 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid, grid2
824 TYPE(pw_r3d_rs_type), POINTER :: pw
825 TYPE(qmmm_per_pot_type), POINTER :: per_pot
826
827 CALL timeset(routinen, handle)
828 NULLIFY (grid)
829 ALLOCATE (lforces(3, num_mm_atoms))
830 lforces = 0.0_dp
831 dr1c = cgrid%pw_grid%dr(1)
832 dr2c = cgrid%pw_grid%dr(2)
833 dr3c = cgrid%pw_grid%dr(3)
834 dvol = cgrid%pw_grid%dvol
835 gbo = cgrid%pw_grid%bounds
836 bo = cgrid%pw_grid%bounds_local
837 grid => cgrid%array
838 IF (par_scheme == do_par_atom) myind = 0
839 radius: DO iradtyp = 1, SIZE(pgfs)
840 per_pot => per_potentials(iradtyp)%pot
841 pw => per_pot%TabLR
842 grid2 => pw%array(:, :, :)
843 npts = pw%pw_grid%npts
844 dr1 = pw%pw_grid%dr(1)
845 dr2 = pw%pw_grid%dr(2)
846 dr3 = pw%pw_grid%dr(3)
847 dr1i = 1.0_dp/dr1
848 dr2i = 1.0_dp/dr2
849 dr3i = 1.0_dp/dr3
850
851 !$OMP PARALLEL DO DEFAULT(NONE) &
852 !$OMP SHARED(bo, grid, grid2, pw, npts, gbo, per_pot, mm_atom_index) &
853 !$OMP SHARED(dr1, dr2, dr3, dr1i, dr2i, dr3i, dr1c, dr2c, dr3c, par_scheme, mm_charges) &
854 !$OMP SHARED(mm_cell, dOmmOqm, dvol, shells, para_env, IRadTyp) &
855 !$OMP SHARED(qmmm_spherical_cutoff, mm_particles, Forces, LForces) &
856 !$OMP PRIVATE(qt, Imm, LIndMM, IndMM, sph_chrg_factor, ra, myind) &
857 !$OMP PRIVATE(rt1, rt2, rt3, ft1, ft2, ft3, my_k, my_j, my_i, xs3, xs2, xs1) &
858 !$OMP PRIVATE(rv3, rv2, rv1, vec, ivec, ik1, ik2, ik3, ik4, xd3, xd2, xd1) &
859 !$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, u1, u2, u3, v1o, v2o, v3o, v4o) &
860 !$OMP PRIVATE(v1d, v2d, v3d, v4d, ij1, ij2, ij3, ij4, e1, e2, e3, f1, f2, f3) &
861 !$OMP PRIVATE(g1, g2, g3, h1, h2, h3, s1o, s2o, s3o, s4o, s1d, s2d, s3d, s4d) &
862 !$OMP PRIVATE(ii1, ii2, ii3, ii4, a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3) &
863 !$OMP PRIVATE(t1o, t2o, t3o, t4o, t1d, t2d, t3d, t4d, t1, t2, t3, t4, s1, s2, s3, s4) &
864 !$OMP PRIVATE(v1, v2, v3, v4, abc_x, abc_x_y, val, fac)
865 atoms: DO imm = 1, SIZE(per_pot%mm_atom_index)
866 IF (par_scheme == do_par_atom) THEN
867 myind = imm + (iradtyp - 1)*SIZE(per_pot%mm_atom_index)
868 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
869 END IF
870 lindmm = per_pot%mm_atom_index(imm)
871 indmm = mm_atom_index(lindmm)
872 IF (shells) THEN
873 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
874 ELSE
875 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
876 END IF
877 qt = mm_charges(lindmm)
878 ! Possible Spherical Cutoff
879 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
880 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
881 qt = qt*sph_chrg_factor
882 END IF
883 IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
884 rt1 = ra(1)
885 rt2 = ra(2)
886 rt3 = ra(3)
887 ft1 = 0.0_dp
888 ft2 = 0.0_dp
889 ft3 = 0.0_dp
890 loopongrid: DO k = bo(1, 3), bo(2, 3)
891 my_k = k - gbo(1, 3)
892 xs3 = real(my_k, dp)*dr3c
893 my_j = bo(1, 2) - gbo(1, 2)
894 xs2 = real(my_j, dp)*dr2c
895 rv3 = rt3 - xs3
896 vec(3) = rv3
897 ivec(3) = floor(vec(3)/pw%pw_grid%dr(3))
898 ik1 = modulo(ivec(3) - 1, npts(3)) + 1
899 ik2 = modulo(ivec(3), npts(3)) + 1
900 ik3 = modulo(ivec(3) + 1, npts(3)) + 1
901 ik4 = modulo(ivec(3) + 2, npts(3)) + 1
902 xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
903 p1 = 3.0_dp + xd3
904 p2 = p1*p1
905 p3 = p2*p1
906 q1 = 2.0_dp + xd3
907 q2 = q1*q1
908 q3 = q2*q1
909 r1 = 1.0_dp + xd3
910 r2 = r1*r1
911 r3 = r2*r1
912 u1 = xd3
913 u2 = u1*u1
914 u3 = u2*u1
915 v1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
916 v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
917 v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
918 v4o = 1.0_dp/6.0_dp*u3
919 v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2
920 v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2
921 v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2
922 v4d = 0.5_dp*u2
923 DO j = bo(1, 2), bo(2, 2)
924 my_i = bo(1, 1) - gbo(1, 1)
925 xs1 = real(my_i, dp)*dr1c
926 rv2 = rt2 - xs2
927 vec(2) = rv2
928 ivec(2) = floor(vec(2)/pw%pw_grid%dr(2))
929 ij1 = modulo(ivec(2) - 1, npts(2)) + 1
930 ij2 = modulo(ivec(2), npts(2)) + 1
931 ij3 = modulo(ivec(2) + 1, npts(2)) + 1
932 ij4 = modulo(ivec(2) + 2, npts(2)) + 1
933 xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
934 e1 = 3.0_dp + xd2
935 e2 = e1*e1
936 e3 = e2*e1
937 f1 = 2.0_dp + xd2
938 f2 = f1*f1
939 f3 = f2*f1
940 g1 = 1.0_dp + xd2
941 g2 = g1*g1
942 g3 = g2*g1
943 h1 = xd2
944 h2 = h1*h1
945 h3 = h2*h1
946 s1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
947 s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
948 s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
949 s4o = 1.0_dp/6.0_dp*h3
950 s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2
951 s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2
952 s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2
953 s4d = 0.5_dp*h2
954 DO i = bo(1, 1), bo(2, 1)
955 rv1 = rt1 - xs1
956 vec(1) = rv1
957 ivec(1) = floor(vec(1)/pw%pw_grid%dr(1))
958 ii1 = modulo(ivec(1) - 1, npts(1)) + 1
959 ii2 = modulo(ivec(1), npts(1)) + 1
960 ii3 = modulo(ivec(1) + 1, npts(1)) + 1
961 ii4 = modulo(ivec(1) + 2, npts(1)) + 1
962 xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
963 a1 = 3.0_dp + xd1
964 a2 = a1*a1
965 a3 = a2*a1
966 b1 = 2.0_dp + xd1
967 b2 = b1*b1
968 b3 = b2*b1
969 c1 = 1.0_dp + xd1
970 c2 = c1*c1
971 c3 = c2*c1
972 d1 = xd1
973 d2 = d1*d1
974 d3 = d2*d1
975 t1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
976 t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
977 t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
978 t4o = 1.0_dp/6.0_dp*d3
979 t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2
980 t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2
981 t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2
982 t4d = 0.5_dp*d2
983
984 t1 = t1d*dr1i
985 t2 = t2d*dr1i
986 t3 = t3d*dr1i
987 t4 = t4d*dr1i
988 s1 = s1o
989 s2 = s2o
990 s3 = s3o
991 s4 = s4o
992 v1 = v1o
993 v2 = v2o
994 v3 = v3o
995 v4 = v4o
996
997 abc_x(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
998 abc_x(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
999 abc_x(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
1000 abc_x(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
1001 abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1002
1003 abc_x(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
1004 abc_x(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
1005 abc_x(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
1006 abc_x(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
1007 abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1008
1009 abc_x(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
1010 abc_x(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
1011 abc_x(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
1012 abc_x(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
1013 abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1014
1015 abc_x(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
1016 abc_x(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
1017 abc_x(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
1018 abc_x(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
1019 abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1020
1021 val(1) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1022
1023 t1 = t1o
1024 t2 = t2o
1025 t3 = t3o
1026 t4 = t4o
1027 s1 = s1d*dr2i
1028 s2 = s2d*dr2i
1029 s3 = s3d*dr2i
1030 s4 = s4d*dr2i
1031
1032 abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1033 abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1034 abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1035 abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1036
1037 val(2) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1038
1039 t1 = t1o
1040 t2 = t2o
1041 t3 = t3o
1042 t4 = t4o
1043 s1 = s1o
1044 s2 = s2o
1045 s3 = s3o
1046 s4 = s4o
1047 v1 = v1d*dr3i
1048 v2 = v2d*dr3i
1049 v3 = v3d*dr3i
1050 v4 = v4d*dr3i
1051
1052 abc_x(1, 1) = grid2(ii1, ij1, ik1)*v1 + grid2(ii1, ij1, ik2)*v2 + grid2(ii1, ij1, ik3)*v3 + grid2(ii1, ij1, ik4)*v4
1053 abc_x(2, 1) = grid2(ii2, ij1, ik1)*v1 + grid2(ii2, ij1, ik2)*v2 + grid2(ii2, ij1, ik3)*v3 + grid2(ii2, ij1, ik4)*v4
1054 abc_x(3, 1) = grid2(ii3, ij1, ik1)*v1 + grid2(ii3, ij1, ik2)*v2 + grid2(ii3, ij1, ik3)*v3 + grid2(ii3, ij1, ik4)*v4
1055 abc_x(4, 1) = grid2(ii4, ij1, ik1)*v1 + grid2(ii4, ij1, ik2)*v2 + grid2(ii4, ij1, ik3)*v3 + grid2(ii4, ij1, ik4)*v4
1056 abc_x_y(1) = abc_x(1, 1)*t1 + abc_x(2, 1)*t2 + abc_x(3, 1)*t3 + abc_x(4, 1)*t4
1057 abc_x(1, 2) = grid2(ii1, ij2, ik1)*v1 + grid2(ii1, ij2, ik2)*v2 + grid2(ii1, ij2, ik3)*v3 + grid2(ii1, ij2, ik4)*v4
1058 abc_x(2, 2) = grid2(ii2, ij2, ik1)*v1 + grid2(ii2, ij2, ik2)*v2 + grid2(ii2, ij2, ik3)*v3 + grid2(ii2, ij2, ik4)*v4
1059 abc_x(3, 2) = grid2(ii3, ij2, ik1)*v1 + grid2(ii3, ij2, ik2)*v2 + grid2(ii3, ij2, ik3)*v3 + grid2(ii3, ij2, ik4)*v4
1060 abc_x(4, 2) = grid2(ii4, ij2, ik1)*v1 + grid2(ii4, ij2, ik2)*v2 + grid2(ii4, ij2, ik3)*v3 + grid2(ii4, ij2, ik4)*v4
1061 abc_x_y(2) = abc_x(1, 2)*t1 + abc_x(2, 2)*t2 + abc_x(3, 2)*t3 + abc_x(4, 2)*t4
1062 abc_x(1, 3) = grid2(ii1, ij3, ik1)*v1 + grid2(ii1, ij3, ik2)*v2 + grid2(ii1, ij3, ik3)*v3 + grid2(ii1, ij3, ik4)*v4
1063 abc_x(2, 3) = grid2(ii2, ij3, ik1)*v1 + grid2(ii2, ij3, ik2)*v2 + grid2(ii2, ij3, ik3)*v3 + grid2(ii2, ij3, ik4)*v4
1064 abc_x(3, 3) = grid2(ii3, ij3, ik1)*v1 + grid2(ii3, ij3, ik2)*v2 + grid2(ii3, ij3, ik3)*v3 + grid2(ii3, ij3, ik4)*v4
1065 abc_x(4, 3) = grid2(ii4, ij3, ik1)*v1 + grid2(ii4, ij3, ik2)*v2 + grid2(ii4, ij3, ik3)*v3 + grid2(ii4, ij3, ik4)*v4
1066 abc_x_y(3) = abc_x(1, 3)*t1 + abc_x(2, 3)*t2 + abc_x(3, 3)*t3 + abc_x(4, 3)*t4
1067 abc_x(1, 4) = grid2(ii1, ij4, ik1)*v1 + grid2(ii1, ij4, ik2)*v2 + grid2(ii1, ij4, ik3)*v3 + grid2(ii1, ij4, ik4)*v4
1068 abc_x(2, 4) = grid2(ii2, ij4, ik1)*v1 + grid2(ii2, ij4, ik2)*v2 + grid2(ii2, ij4, ik3)*v3 + grid2(ii2, ij4, ik4)*v4
1069 abc_x(3, 4) = grid2(ii3, ij4, ik1)*v1 + grid2(ii3, ij4, ik2)*v2 + grid2(ii3, ij4, ik3)*v3 + grid2(ii3, ij4, ik4)*v4
1070 abc_x(4, 4) = grid2(ii4, ij4, ik1)*v1 + grid2(ii4, ij4, ik2)*v2 + grid2(ii4, ij4, ik3)*v3 + grid2(ii4, ij4, ik4)*v4
1071 abc_x_y(4) = abc_x(1, 4)*t1 + abc_x(2, 4)*t2 + abc_x(3, 4)*t3 + abc_x(4, 4)*t4
1072
1073 val(3) = abc_x_y(1)*s1 + abc_x_y(2)*s2 + abc_x_y(3)*s3 + abc_x_y(4)*s4
1074
1075 fac = grid(i, j, k)
1076 ft1 = ft1 + val(1)*fac
1077 ft2 = ft2 + val(2)*fac
1078 ft3 = ft3 + val(3)*fac
1079 xs1 = xs1 + dr1c
1080 END DO
1081 xs2 = xs2 + dr2c
1082 END DO
1083 END DO loopongrid
1084 qt = -qt*dvol
1085 lforces(1, lindmm) = ft1*qt
1086 lforces(2, lindmm) = ft2*qt
1087 lforces(3, lindmm) = ft3*qt
1088
1089 forces(1, lindmm) = forces(1, lindmm) + lforces(1, lindmm)
1090 forces(2, lindmm) = forces(2, lindmm) + lforces(2, lindmm)
1091 forces(3, lindmm) = forces(3, lindmm) + lforces(3, lindmm)
1092 END DO atoms
1093 !$OMP END PARALLEL DO
1094 END DO radius
1095 !
1096 ! Debug Statement
1097 !
1098 IF (debug_this_module) THEN
1099 CALL debug_qmmm_forces_with_gauss_lg(pgfs=pgfs, &
1100 aug_pools=aug_pools, &
1101 rho=cgrid, &
1102 num_mm_atoms=num_mm_atoms, &
1103 mm_charges=mm_charges, &
1104 mm_atom_index=mm_atom_index, &
1105 mm_particles=mm_particles, &
1106 coarser_grid_level=coarser_grid_level, &
1107 debug_force=lforces, &
1108 per_potentials=per_potentials, &
1109 para_env=para_env, &
1110 mm_cell=mm_cell, &
1111 dommoqm=dommoqm, &
1112 iw=iw, &
1113 par_scheme=par_scheme, &
1114 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1115 shells=shells)
1116 END IF
1117 DEALLOCATE (lforces)
1118 CALL timestop(handle)
1119 END SUBROUTINE qmmm_forces_with_gaussian_lg
1120
1121! **************************************************************************************************
1122!> \brief Evaluates the contribution to the forces due to the Long Range
1123!> part of the QM/MM potential computed collocating the Electrostatic
1124!> Gaussian Potential.
1125!> \param pgfs ...
1126!> \param cgrid ...
1127!> \param num_mm_atoms ...
1128!> \param mm_charges ...
1129!> \param mm_atom_index ...
1130!> \param mm_particles ...
1131!> \param para_env ...
1132!> \param coarser_grid_level ...
1133!> \param Forces ...
1134!> \param potentials ...
1135!> \param aug_pools ...
1136!> \param mm_cell ...
1137!> \param dOmmOqm ...
1138!> \param iw ...
1139!> \param par_scheme ...
1140!> \param qmmm_spherical_cutoff ...
1141!> \param shells ...
1142!> \par History
1143!> 08.2004 created [tlaino]
1144!> \author Teodoro Laino
1145! **************************************************************************************************
1146 SUBROUTINE qmmm_forces_with_gaussian_lr(pgfs, cgrid, num_mm_atoms, mm_charges, mm_atom_index, &
1147 mm_particles, para_env, coarser_grid_level, Forces, potentials, &
1148 aug_pools, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1149 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1150 TYPE(pw_r3d_rs_type), INTENT(IN) :: cgrid
1151 INTEGER, INTENT(IN) :: num_mm_atoms
1152 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1153 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1154 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1155 TYPE(mp_para_env_type), POINTER :: para_env
1156 INTEGER, INTENT(IN) :: coarser_grid_level
1157 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
1158 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1159 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1160 TYPE(cell_type), POINTER :: mm_cell
1161 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1162 INTEGER, INTENT(IN) :: iw, par_scheme
1163 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1164 LOGICAL :: shells
1165
1166 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian_LR'
1167
1168 INTEGER :: handle, i, imm, indmm, iradtyp, ix, j, &
1169 k, lindmm, my_i, my_j, my_k, myind, &
1170 n1, n2, n3
1171 INTEGER, DIMENSION(2, 3) :: bo, gbo
1172 REAL(kind=dp) :: dr1, dr2, dr3, dvol, dx, fac, ft1, ft2, &
1173 ft3, qt, r, r2, rd1, rd2, rd3, rt1, &
1174 rt2, rt3, rv1, rv2, rv3, rx, rx2, &
1175 sph_chrg_factor, term, xs1, xs2, xs3
1176 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: lforces
1177 REAL(kind=dp), DIMENSION(3) :: ra
1178 REAL(kind=dp), DIMENSION(:, :), POINTER :: pot0_2
1179 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
1180 TYPE(qmmm_pot_type), POINTER :: pot
1181
1182 CALL timeset(routinen, handle)
1183 ALLOCATE (lforces(3, num_mm_atoms))
1184 lforces = 0.0_dp
1185 n1 = cgrid%pw_grid%npts(1)
1186 n2 = cgrid%pw_grid%npts(2)
1187 n3 = cgrid%pw_grid%npts(3)
1188 dr1 = cgrid%pw_grid%dr(1)
1189 dr2 = cgrid%pw_grid%dr(2)
1190 dr3 = cgrid%pw_grid%dr(3)
1191 dvol = cgrid%pw_grid%dvol
1192 gbo = cgrid%pw_grid%bounds
1193 bo = cgrid%pw_grid%bounds_local
1194 grid => cgrid%array
1195 IF (par_scheme == do_par_atom) myind = 0
1196 radius: DO iradtyp = 1, SIZE(pgfs)
1197 pot => potentials(iradtyp)%pot
1198 dx = pot%dx
1199 pot0_2 => pot%pot0_2
1200 !$OMP PARALLEL DO DEFAULT(NONE) &
1201 !$OMP SHARED(pot, par_scheme, para_env, dvol, mm_atom_index, mm_particles, dOmmOqm) &
1202 !$OMP SHARED(mm_cell, mm_charges, dx, LForces, Forces, qmmm_spherical_cutoff, shells, dr1, dr2, dr3, gbo, bo) &
1203 !$OMP SHARED(IRadTyp, pot0_2, grid) &
1204 !$OMP PRIVATE(Imm, myind, ra, LIndMM, IndMM, qt, rt1, rt2, rt3, ft1, ft2, ft3, i, j, k, sph_chrg_factor) &
1205 !$OMP PRIVATE(my_k, my_j, my_i, xs3, xs2, xs1, rv1, rv2, rv3, r, ix, rx, rx2, r2, Term, fac) &
1206 !$OMP PRIVATE(rd1, rd2, rd3)
1207 atoms: DO imm = 1, SIZE(pot%mm_atom_index)
1208 IF (par_scheme == do_par_atom) THEN
1209 myind = imm + (iradtyp - 1)*SIZE(pot%mm_atom_index)
1210 IF (mod(myind, para_env%num_pe) /= para_env%mepos) cycle atoms
1211 END IF
1212 lindmm = pot%mm_atom_index(imm)
1213 indmm = mm_atom_index(lindmm)
1214 ra(:) = pbc(mm_particles(indmm)%r - dommoqm, mm_cell) + dommoqm
1215 IF (shells) THEN
1216 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
1217 END IF
1218 qt = mm_charges(lindmm)
1219 ! Possible Spherical Cutoff
1220 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1221 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
1222 qt = qt*sph_chrg_factor
1223 END IF
1224 IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
1225 rt1 = ra(1)
1226 rt2 = ra(2)
1227 rt3 = ra(3)
1228 ft1 = 0.0_dp
1229 ft2 = 0.0_dp
1230 ft3 = 0.0_dp
1231 loopongrid: DO k = bo(1, 3), bo(2, 3)
1232 my_k = k - gbo(1, 3)
1233 xs3 = real(my_k, dp)*dr3
1234 my_j = bo(1, 2) - gbo(1, 2)
1235 xs2 = real(my_j, dp)*dr2
1236 rv3 = rt3 - xs3
1237 DO j = bo(1, 2), bo(2, 2)
1238 my_i = bo(1, 1) - gbo(1, 1)
1239 xs1 = real(my_i, dp)*dr1
1240 rv2 = rt2 - xs2
1241 DO i = bo(1, 1), bo(2, 1)
1242 rv1 = rt1 - xs1
1243 r2 = rv1*rv1 + rv2*rv2 + rv3*rv3
1244 r = sqrt(r2)
1245 ix = floor(r/dx) + 1
1246 rx = (r - real(ix - 1, dp)*dx)/dx
1247 rx2 = rx*rx
1248 term = pot0_2(1, ix)*(-6.0_dp*(rx - rx2)) &
1249 + pot0_2(2, ix)*(1.0_dp - 4.0_dp*rx + 3.0_dp*rx2) &
1250 + pot0_2(1, ix + 1)*(6.0_dp*(rx - rx2)) &
1251 + pot0_2(2, ix + 1)*(-2.0_dp*rx + 3.0_dp*rx2)
1252 fac = grid(i, j, k)*term
1253 IF (r == 0.0_dp) THEN
1254 rd1 = 1.0_dp
1255 rd2 = 1.0_dp
1256 rd3 = 1.0_dp
1257 ELSE
1258 rd1 = rv1/r
1259 rd2 = rv2/r
1260 rd3 = rv3/r
1261 END IF
1262 ft1 = ft1 + fac*rd1
1263 ft2 = ft2 + fac*rd2
1264 ft3 = ft3 + fac*rd3
1265 xs1 = xs1 + dr1
1266 END DO
1267 xs2 = xs2 + dr2
1268 END DO
1269 END DO loopongrid
1270 qt = -qt*dvol/dx
1271 lforces(1, lindmm) = ft1*qt
1272 lforces(2, lindmm) = ft2*qt
1273 lforces(3, lindmm) = ft3*qt
1274
1275 forces(1, lindmm) = forces(1, lindmm) + lforces(1, lindmm)
1276 forces(2, lindmm) = forces(2, lindmm) + lforces(2, lindmm)
1277 forces(3, lindmm) = forces(3, lindmm) + lforces(3, lindmm)
1278 END DO atoms
1279 !$OMP END PARALLEL DO
1280 END DO radius
1281 !
1282 ! Debug Statement
1283 !
1284 IF (debug_this_module) THEN
1285 CALL debug_qmmm_forces_with_gauss_lr(pgfs=pgfs, &
1286 aug_pools=aug_pools, &
1287 rho=cgrid, &
1288 num_mm_atoms=num_mm_atoms, &
1289 mm_charges=mm_charges, &
1290 mm_atom_index=mm_atom_index, &
1291 mm_particles=mm_particles, &
1292 coarser_grid_level=coarser_grid_level, &
1293 debug_force=lforces, &
1294 potentials=potentials, &
1295 para_env=para_env, &
1296 mm_cell=mm_cell, &
1297 dommoqm=dommoqm, &
1298 iw=iw, &
1299 par_scheme=par_scheme, &
1300 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1301 shells=shells)
1302 END IF
1303
1304 DEALLOCATE (lforces)
1305 CALL timestop(handle)
1306 END SUBROUTINE qmmm_forces_with_gaussian_lr
1307
1308! **************************************************************************************************
1309!> \brief Evaluates numerically QM/MM forces and compares them with
1310!> the analytically computed ones.
1311!> It is evaluated only when debug_this_module is set to .TRUE.
1312!> \param rho ...
1313!> \param qs_env ...
1314!> \param qmmm_env ...
1315!> \param Analytical_Forces ...
1316!> \param mm_particles ...
1317!> \param mm_atom_index ...
1318!> \param num_mm_atoms ...
1319!> \param interp_section ...
1320!> \param mm_cell ...
1321!> \par History
1322!> 08.2004 created [tlaino]
1323!> \author Teodoro Laino
1324! **************************************************************************************************
1325 SUBROUTINE qmmm_debug_forces(rho, qs_env, qmmm_env, Analytical_Forces, &
1326 mm_particles, mm_atom_index, num_mm_atoms, &
1327 interp_section, mm_cell)
1328 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1329 TYPE(qs_environment_type), POINTER :: qs_env
1330 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1331 REAL(kind=dp), DIMENSION(:, :), POINTER :: analytical_forces
1332 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1333 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1334 INTEGER, INTENT(IN) :: num_mm_atoms
1335 TYPE(section_vals_type), POINTER :: interp_section
1336 TYPE(cell_type), POINTER :: mm_cell
1337
1338 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_debug_forces'
1339
1340 INTEGER :: handle, i, indmm, iw, j, k
1341 REAL(kind=dp) :: coord_save
1342 REAL(kind=dp), DIMENSION(2) :: energy
1343 REAL(kind=dp), DIMENSION(3) :: err
1344 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1345 TYPE(cp_logger_type), POINTER :: logger
1346 TYPE(mp_para_env_type), POINTER :: para_env
1347 TYPE(pw_env_type), POINTER :: pw_env
1348 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1349 TYPE(pw_r3d_rs_type) :: v_qmmm_rspace
1350 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
1351 TYPE(section_vals_type), POINTER :: input_section, print_section
1352
1353 CALL timeset(routinen, handle)
1354 NULLIFY (num_forces)
1355 CALL get_qs_env(qs_env=qs_env, &
1356 pw_env=pw_env, &
1357 input=input_section, &
1358 para_env=para_env)
1359
1360 print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
1361 logger => cp_get_default_logger()
1362 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".qmmmLog")
1363 CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools)
1364 CALL pw_pools(1)%pool%create_pw(v_qmmm_rspace)
1365 ALLOCATE (num_forces(3, num_mm_atoms))
1366 ks_qmmm_env_loc => qs_env%ks_qmmm_env
1367 IF (iw > 0) WRITE (iw, '(/A)') "DEBUG SECTION:"
1368 atoms: DO i = 1, num_mm_atoms
1369 indmm = mm_atom_index(i)
1370 coords: DO j = 1, 3
1371 coord_save = mm_particles(indmm)%r(j)
1372 energy = 0.0_dp
1373 diff: DO k = 1, 2
1374 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1375 CALL pw_zero(v_qmmm_rspace)
1376 SELECT CASE (qmmm_env%qmmm_coupl_type)
1377 CASE (do_qmmm_coulomb)
1378 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1379 CASE (do_qmmm_pcharge)
1380 cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1381 CASE (do_qmmm_gauss, do_qmmm_swave)
1382 CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env, &
1383 v_qmmm=v_qmmm_rspace, &
1384 mm_particles=mm_particles, &
1385 aug_pools=qmmm_env%aug_pools, &
1386 para_env=para_env, &
1387 eps_mm_rspace=qmmm_env%eps_mm_rspace, &
1388 cube_info=ks_qmmm_env_loc%cube_info, &
1389 pw_pools=pw_pools, &
1390 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
1391 coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
1392 interp_section=interp_section, &
1393 mm_cell=mm_cell)
1394 CASE (do_qmmm_none)
1395 cycle diff
1396 CASE DEFAULT
1397 cpabort("Unknown QM/MM coupling")
1398 END SELECT
1399 energy(k) = pw_integral_ab(rho, v_qmmm_rspace)
1400 END DO diff
1401 IF (iw > 0) THEN
1402 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1403 "DEBUG :: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1404 END IF
1405 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1406 mm_particles(indmm)%r(j) = coord_save
1407 END DO coords
1408 END DO atoms
1409
1410 SELECT CASE (qmmm_env%qmmm_coupl_type)
1411 CASE (do_qmmm_coulomb)
1412 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1413 CASE (do_qmmm_pcharge)
1414 cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1415 CASE (do_qmmm_gauss, do_qmmm_swave)
1416 IF (iw > 0) WRITE (iw, '(/A/)') "CHECKING NUMERICAL Vs ANALYTICAL FORCES (Err%):"
1417 DO i = 1, num_mm_atoms
1418 indmm = mm_atom_index(i)
1419 err = 0.0_dp
1420 DO k = 1, 3
1421 IF (abs(num_forces(k, i)) >= 5.0e-5_dp) THEN
1422 err(k) = (analytical_forces(k, i) - num_forces(k, i))/num_forces(k, i)*100.0_dp
1423 END IF
1424 END DO
1425 IF (iw > 0) THEN
1426 WRITE (iw, 100) indmm, analytical_forces(1, i), num_forces(1, i), err(1), &
1427 analytical_forces(2, i), num_forces(2, i), err(2), &
1428 analytical_forces(3, i), num_forces(3, i), err(3)
1429 END IF
1430 cpassert(abs(err(1)) <= maxerr)
1431 cpassert(abs(err(2)) <= maxerr)
1432 cpassert(abs(err(3)) <= maxerr)
1433 END DO
1434 CASE (do_qmmm_none)
1435 IF (iw > 0) WRITE (iw, '(T3,A)') "No QM/MM Derivatives to debug. Just Mechanical Coupling!"
1436 CASE DEFAULT
1437 cpabort("Unknown QM/MM coupling")
1438 END SELECT
1439 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
1440
1441 CALL pw_pools(1)%pool%give_back_pw(v_qmmm_rspace)
1442 DEALLOCATE (num_forces)
1443 CALL timestop(handle)
1444100 FORMAT(i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1445 END SUBROUTINE qmmm_debug_forces
1446
1447! **************************************************************************************************
1448!> \brief Debugs the integrate_gf_rspace_NoPBC.. It may helps ;-P
1449!> \param ilevel ...
1450!> \param zetp ...
1451!> \param rp ...
1452!> \param W ...
1453!> \param pwgrid ...
1454!> \param cube_info ...
1455!> \param eps_mm_rspace ...
1456!> \param aug_pools ...
1457!> \param debug_force ...
1458!> \param mm_cell ...
1459!> \param auxbas_grid ...
1460!> \param n_rep_real ...
1461!> \param iw ...
1462!> \par History
1463!> 08.2004 created [tlaino]
1464!> \author Teodoro Laino
1465! **************************************************************************************************
1466 SUBROUTINE debug_integrate_gf_rspace_nopbc(ilevel, zetp, rp, W, pwgrid, cube_info, &
1467 eps_mm_rspace, aug_pools, debug_force, &
1468 mm_cell, auxbas_grid, n_rep_real, iw)
1469 INTEGER, INTENT(IN) :: ilevel
1470 REAL(kind=dp), INTENT(IN) :: zetp
1471 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rp
1472 REAL(kind=dp), INTENT(IN) :: w
1473 TYPE(pw_r3d_rs_type), INTENT(IN) :: pwgrid
1474 TYPE(cube_info_type), INTENT(IN) :: cube_info
1475 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
1476 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1477 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: debug_force
1478 TYPE(cell_type), POINTER :: mm_cell
1479 INTEGER, INTENT(IN) :: auxbas_grid
1480 INTEGER, DIMENSION(3), INTENT(IN) :: n_rep_real
1481 INTEGER, INTENT(IN) :: iw
1482
1483 CHARACTER(len=*), PARAMETER :: routinen = 'debug_integrate_gf_rspace_NoPBC'
1484
1485 INTEGER :: handle, i, igrid, k, ngrids
1486 INTEGER, DIMENSION(2, 3) :: bo2
1487 INTEGER, SAVE :: icount
1488 REAL(kind=dp), DIMENSION(2) :: energy
1489 REAL(kind=dp), DIMENSION(3) :: err, force, myrp
1490 REAL(kind=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat
1491 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1492
1493 DATA icount/0/
1494 ! Statements
1495 CALL timeset(routinen, handle)
1496 !Statements
1497 ngrids = SIZE(aug_pools)
1498 CALL pw_pools_create_pws(aug_pools, grids)
1499 DO igrid = 1, ngrids
1500 CALL pw_zero(grids(igrid))
1501 END DO
1502 bo2 = grids(auxbas_grid)%pw_grid%bounds
1503 ALLOCATE (xdat(bo2(1, 1):bo2(2, 1)))
1504 ALLOCATE (ydat(bo2(1, 2):bo2(2, 2)))
1505 ALLOCATE (zdat(bo2(1, 3):bo2(2, 3)))
1506
1507 icount = icount + 1
1508 DO i = 1, 3
1509 DO k = 1, 2
1510 myrp = rp
1511 myrp(i) = myrp(i) + (-1.0_dp)**k*dx
1512 CALL pw_zero(grids(ilevel))
1513 CALL collocate_gf_rspace_nopbc(zetp=zetp, &
1514 rp=myrp, &
1515 scale=-1.0_dp, &
1516 w=w, &
1517 pwgrid=grids(ilevel), &
1518 cube_info=cube_info, &
1519 eps_mm_rspace=eps_mm_rspace, &
1520 xdat=xdat, &
1521 ydat=ydat, &
1522 zdat=zdat, &
1523 bo2=bo2, &
1524 n_rep_real=n_rep_real, &
1525 mm_cell=mm_cell)
1526
1527 energy(k) = pw_integral_ab(pwgrid, grids(ilevel))
1528 END DO
1529 force(i) = (energy(2) - energy(1))/(2.0_dp*dx)
1530 END DO
1531 err = 0.0_dp
1532 IF (all(force /= 0.0_dp)) THEN
1533 err(1) = (debug_force(1) - force(1))/force(1)*100.0_dp
1534 err(2) = (debug_force(2) - force(2))/force(2)*100.0_dp
1535 err(3) = (debug_force(3) - force(3))/force(3)*100.0_dp
1536 END IF
1537 IF (iw > 0) THEN
1538 WRITE (iw, 100) icount, debug_force(1), force(1), err(1), &
1539 debug_force(2), force(2), err(2), &
1540 debug_force(3), force(3), err(3)
1541 END IF
1542 cpassert(abs(err(1)) <= maxerr)
1543 cpassert(abs(err(2)) <= maxerr)
1544 cpassert(abs(err(3)) <= maxerr)
1545
1546 IF (ASSOCIATED(xdat)) THEN
1547 DEALLOCATE (xdat)
1548 END IF
1549 IF (ASSOCIATED(ydat)) THEN
1550 DEALLOCATE (ydat)
1551 END IF
1552 IF (ASSOCIATED(zdat)) THEN
1553 DEALLOCATE (zdat)
1554 END IF
1555
1556 CALL pw_pools_give_back_pws(aug_pools, grids)
1557 CALL timestop(handle)
1558100 FORMAT("Collocation : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1559 END SUBROUTINE debug_integrate_gf_rspace_nopbc
1560
1561! **************************************************************************************************
1562!> \brief Debugs qmmm_forces_with_gaussian_LG.. It may helps too ... ;-]
1563!> \param pgfs ...
1564!> \param aug_pools ...
1565!> \param rho ...
1566!> \param mm_charges ...
1567!> \param mm_atom_index ...
1568!> \param mm_particles ...
1569!> \param num_mm_atoms ...
1570!> \param coarser_grid_level ...
1571!> \param per_potentials ...
1572!> \param debug_force ...
1573!> \param para_env ...
1574!> \param mm_cell ...
1575!> \param dOmmOqm ...
1576!> \param iw ...
1577!> \param par_scheme ...
1578!> \param qmmm_spherical_cutoff ...
1579!> \param shells ...
1580!> \par History
1581!> 08.2004 created [tlaino]
1582!> \author Teodoro Laino
1583! **************************************************************************************************
1584 SUBROUTINE debug_qmmm_forces_with_gauss_lg(pgfs, aug_pools, rho, mm_charges, mm_atom_index, &
1585 mm_particles, num_mm_atoms, coarser_grid_level, per_potentials, &
1586 debug_force, para_env, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1587
1588 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1589 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1590 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1591 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1592 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1593 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1594 INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level
1595 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
1596 REAL(kind=dp), DIMENSION(:, :) :: debug_force
1597 TYPE(mp_para_env_type), POINTER :: para_env
1598 TYPE(cell_type), POINTER :: mm_cell
1599 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1600 INTEGER, INTENT(IN) :: iw, par_scheme
1601 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1602 LOGICAL :: shells
1603
1604 CHARACTER(len=*), PARAMETER :: routinen = 'debug_qmmm_forces_with_gauss_LG'
1605
1606 INTEGER :: handle, i, igrid, indmm, j, k, ngrids
1607 REAL(kind=dp) :: coord_save
1608 REAL(kind=dp), DIMENSION(2) :: energy
1609 REAL(kind=dp), DIMENSION(3) :: err
1610 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1611 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1612
1613 ALLOCATE (num_forces(3, num_mm_atoms))
1614 CALL timeset(routinen, handle)
1615 ngrids = SIZE(aug_pools)
1616 CALL pw_pools_create_pws(aug_pools, grids)
1617 DO igrid = 1, ngrids
1618 CALL pw_zero(grids(igrid))
1619 END DO
1620 atoms: DO i = 1, num_mm_atoms
1621 indmm = mm_atom_index(i)
1622 coords: DO j = 1, 3
1623 coord_save = mm_particles(indmm)%r(j)
1624 energy = 0.0_dp
1625 diff: DO k = 1, 2
1626 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1627 CALL pw_zero(grids(coarser_grid_level))
1628
1629 CALL qmmm_elec_with_gaussian_lg(pgfs=pgfs, &
1630 cgrid=grids(coarser_grid_level), &
1631 mm_charges=mm_charges, &
1632 mm_atom_index=mm_atom_index, &
1633 mm_particles=mm_particles, &
1634 para_env=para_env, &
1635 per_potentials=per_potentials, &
1636 mm_cell=mm_cell, &
1637 dommoqm=dommoqm, &
1638 par_scheme=par_scheme, &
1639 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1640 shells=shells)
1641
1642 energy(k) = pw_integral_ab(rho, grids(coarser_grid_level))
1643 END DO diff
1644 IF (iw > 0) THEN
1645 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1646 "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1647 END IF
1648 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1649 mm_particles(indmm)%r(j) = coord_save
1650 END DO coords
1651 END DO atoms
1652
1653 DO i = 1, num_mm_atoms
1654 indmm = mm_atom_index(i)
1655 err = 0.0_dp
1656 IF (all(num_forces /= 0.0_dp)) THEN
1657 err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1658 err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1659 err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1660 END IF
1661 IF (iw > 0) THEN
1662 WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1663 debug_force(2, i), num_forces(2, i), err(2), &
1664 debug_force(3, i), num_forces(3, i), err(3)
1665 END IF
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) THEN
1761 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1762 "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1763 END IF
1764 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1765 mm_particles(indmm)%r(j) = coord_save
1766 END DO coords
1767 END DO atoms
1768
1769 DO i = 1, num_mm_atoms
1770 indmm = mm_atom_index(i)
1771 err = 0.0_dp
1772 IF (all(num_forces(:, i) /= 0.0_dp)) THEN
1773 err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1774 err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1775 err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1776 END IF
1777 IF (iw > 0) THEN
1778 WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1779 debug_force(2, i), num_forces(2, i), err(2), &
1780 debug_force(3, i), num_forces(3, i), err(3)
1781 END IF
1782 cpassert(abs(err(1)) <= maxerr)
1783 cpassert(abs(err(2)) <= maxerr)
1784 cpassert(abs(err(3)) <= maxerr)
1785 END DO
1786
1787 DEALLOCATE (num_forces)
1788 CALL pw_pools_give_back_pws(aug_pools, grids)
1789 CALL timestop(handle)
1790100 FORMAT("MM Atom LR : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1791 END SUBROUTINE debug_qmmm_forces_with_gauss_lr
1792
1793END 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:103
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:296
input_constants::do_qmmm_none
integer, parameter, public do_qmmm_none
Definition input_constants.F:396
input_constants::do_qmmm_pcharge
integer, parameter, public do_qmmm_pcharge
Definition input_constants.F:396
input_constants::do_qmmm_coulomb
integer, parameter, public do_qmmm_coulomb
Definition input_constants.F:396
input_constants::do_qmmm_swave
integer, parameter, public do_qmmm_swave
Definition input_constants.F:396
input_constants::do_qmmm_gauss
integer, parameter, public do_qmmm_gauss
Definition input_constants.F:396
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:160
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: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:242
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:756
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:531
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: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, 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:533
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:757
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:158
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
message_passing::mp_request_type
Definition message_passing.F:606
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