(git:ed6f26b)
Loading...
Searching...
No Matches
qmmm_gpw_forces.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines to compute energy and forces in a QM/MM calculation
10!> \par History
11!> 05.2004 created [tlaino]
12!> \author Teodoro Laino
13! **************************************************************************************************
14MODULE qmmm_gpw_forces
15 USE cell_types, ONLY: cell_type,&
16 pbc
17 USE cp_control_types, ONLY: dft_control_type
18 USE cp_log_handling, ONLY: cp_get_default_logger,&
19 cp_logger_type
20 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
21 cp_print_key_unit_nr
22 USE cp_spline_utils, ONLY: pw_restrict_s3,&
23 spline3_nopbc_interp,&
24 spline3_pbc_interp
25 USE cube_utils, ONLY: cube_info_type
26 USE input_constants, ONLY: do_par_atom,&
27 do_qmmm_coulomb,&
28 do_qmmm_gauss,&
29 do_qmmm_none,&
30 do_qmmm_pcharge,&
31 do_qmmm_swave
32 USE input_section_types, ONLY: section_vals_get_subs_vals,&
33 section_vals_type,&
34 section_vals_val_get
35 USE kinds, ONLY: dp
36 USE message_passing, ONLY: mp_comm_type,&
37 mp_para_env_type,&
38 mp_request_type
39 USE mm_collocate_potential, ONLY: collocate_gf_rspace_nopbc,&
40 integrate_gf_rspace_nopbc
41 USE particle_types, ONLY: particle_type
42 USE pw_env_types, ONLY: pw_env_get,&
43 pw_env_type
44 USE pw_methods, ONLY: pw_axpy,&
45 pw_integral_ab,&
46 pw_transfer,&
47 pw_zero
48 USE pw_pool_types, ONLY: pw_pool_p_type,&
49 pw_pool_type,&
50 pw_pools_create_pws,&
51 pw_pools_give_back_pws
52 USE pw_types, ONLY: pw_c1d_gs_type,&
53 pw_r3d_rs_type
54 USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,&
55 qmmm_gaussian_type
56 USE qmmm_gpw_energy, ONLY: qmmm_elec_with_gaussian,&
57 qmmm_elec_with_gaussian_lg,&
58 qmmm_elec_with_gaussian_lr
59 USE qmmm_se_forces, ONLY: deriv_se_qmmm_matrix
60 USE qmmm_tb_methods, ONLY: deriv_tb_qmmm_matrix,&
61 deriv_tb_qmmm_matrix_pc
62 USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
63 qmmm_per_pot_p_type,&
64 qmmm_per_pot_type,&
65 qmmm_pot_p_type,&
66 qmmm_pot_type
67 USE qmmm_util, ONLY: spherical_cutoff_factor
68 USE qs_energy_types, ONLY: qs_energy_type
69 USE qs_environment_types, ONLY: get_qs_env,&
70 qs_environment_type
71 USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
72 USE qs_rho_types, ONLY: qs_rho_get,&
73 qs_rho_type
74#include "./base/base_uses.f90"
75
76 IMPLICIT NONE
77
78 PRIVATE
79 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
80 REAL(KIND=dp), PARAMETER, PRIVATE :: dx = 0.01_dp ! Debug Variables
81 REAL(KIND=dp), PARAMETER, PRIVATE :: maxerr = 10.0_dp ! Debug Variables
82 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gpw_forces'
83 PUBLIC :: qmmm_forces
84
85CONTAINS
86
87! **************************************************************************************************
88!> \brief General driver to Compute the contribution
89!> to the forces due to the QM/MM potential
90!> \param qs_env ...
91!> \param qmmm_env ...
92!> \param mm_particles ...
93!> \param calc_force ...
94!> \param mm_cell ...
95!> \par History
96!> 06.2004 created [tlaino]
97!> \author Teodoro Laino
98! **************************************************************************************************
99 SUBROUTINE qmmm_forces(qs_env, qmmm_env, mm_particles, calc_force, mm_cell)
100 TYPE(qs_environment_type), POINTER :: qs_env
101 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
102 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
103 LOGICAL, INTENT(in), OPTIONAL :: calc_force
104 TYPE(cell_type), POINTER :: mm_cell
105
106 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces'
107
108 INTEGER :: handle, iatom, image_indmm, imm, indmm, &
109 ispin, iw
110 LOGICAL :: gapw, need_f, periodic
111 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges, &
112 forces_added_shells
113 TYPE(cp_logger_type), POINTER :: logger
114 TYPE(dft_control_type), POINTER :: dft_control
115 TYPE(mp_para_env_type), POINTER :: para_env
116 TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
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), DIMENSION(:), POINTER :: rho_r
121 TYPE(pw_r3d_rs_type), POINTER :: rho_tot_r, rho_tot_r2
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, rho_tot_r, energy, forces, &
133 forces_added_charges, input_section, rho0_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 dft_control=dft_control)
143
144 CALL qs_rho_get(rho, rho_r=rho_r)
145
146 logger => cp_get_default_logger()
147 ks_qmmm_env_loc => qs_env%ks_qmmm_env
148 interp_section => section_vals_get_subs_vals(input_section, "QMMM%INTERPOLATOR")
149 print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
150 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
151 extension=".qmmmLog")
152 gapw = dft_control%qs_control%gapw
153 ! If forces are required allocate these temporary arrays
154 IF (need_f) THEN
155 ALLOCATE (forces(3, qmmm_env%num_mm_atoms))
156 ALLOCATE (forces_added_charges(3, qmmm_env%added_charges%num_mm_atoms))
157 ALLOCATE (forces_added_shells(3, qmmm_env%added_shells%num_mm_atoms))
158 forces(:, :) = 0.0_dp
159 forces_added_charges(:, :) = 0.0_dp
160 forces_added_shells(:, :) = 0.0_dp
161 END IF
162 IF (dft_control%qs_control%semi_empirical) THEN
163 ! SEMIEMPIRICAL
164 SELECT CASE (qmmm_env%qmmm_coupl_type)
165 CASE (do_qmmm_coulomb)
166 CALL deriv_se_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env, &
167 need_f, forces, forces_added_charges)
168 CASE (do_qmmm_pcharge)
169 cpabort("Point Charge QM/MM electrostatic coupling not yet implemented for SE.")
170 CASE (do_qmmm_gauss, do_qmmm_swave)
171 cpabort("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.")
172 CASE (do_qmmm_none)
173 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
174 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
175 CASE DEFAULT
176 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
177 cpabort("")
178 END SELECT
179 ELSEIF (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 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
195 cpabort("")
196 END SELECT
197 IF (need_f) THEN
198 forces(:, :) = forces(:, :)/real(para_env%num_pe, kind=dp)
199 forces_added_charges(:, :) = forces_added_charges(:, :)/real(para_env%num_pe, kind=dp)
200 END IF
201 ELSE
202 ! GPW/GAPW
203 CALL pw_env_get(pw_env=pw_env, &
204 pw_pools=pw_pools, &
205 auxbas_pw_pool=auxbas_pool)
206 NULLIFY (rho_tot_r)
207 ALLOCATE (rho_tot_r)
208 CALL auxbas_pool%create_pw(rho_tot_r)
209 ! IF GAPW the core charge is replaced by the compensation charge
210 IF (gapw) THEN
211 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
212 CALL pw_transfer(rho_core, rho_tot_r)
213 energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
214 NULLIFY (rho_tot_r2)
215 ALLOCATE (rho_tot_r2)
216 CALL auxbas_pool%create_pw(rho_tot_r2)
217 CALL pw_transfer(rho0_s_gs, rho_tot_r2)
218 CALL pw_axpy(rho_tot_r2, rho_tot_r)
219 CALL auxbas_pool%give_back_pw(rho_tot_r2)
220 DEALLOCATE (rho_tot_r2)
221 ELSE
222 CALL pw_transfer(rho0_s_gs, rho_tot_r)
223 !
224 ! QM/MM Nuclear Electrostatic Potential already included through rho0
225 !
226 energy%qmmm_nu = 0.0_dp
227 END IF
228 ELSE
229 CALL pw_transfer(rho_core, rho_tot_r)
230 !
231 ! Computes the QM/MM Nuclear Electrostatic Potential
232 !
233 energy%qmmm_nu = pw_integral_ab(rho_tot_r, ks_qmmm_env_loc%v_qmmm_rspace)
234 END IF
235 IF (need_f) THEN
236 !
237 DO ispin = 1, SIZE(rho_r)
238 CALL pw_axpy(rho_r(ispin), rho_tot_r)
239 END DO
240 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Evaluating forces on MM atoms due to the:"
241 ! Electrostatic Interaction type...
242 SELECT CASE (qmmm_env%qmmm_coupl_type)
243 CASE (do_qmmm_coulomb)
244 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
245 CASE (do_qmmm_pcharge)
246 cpabort("Point Charge QM/MM electrostatic coupling not yet implemented for GPW/GAPW.")
247 CASE (do_qmmm_gauss, do_qmmm_swave)
248 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
249 "- QM/MM Coupling computed collocating the Gaussian Potential Functions."
250 CALL qmmm_forces_with_gaussian(rho=rho_tot_r, &
251 qmmm_env=qmmm_env, &
252 mm_particles=mm_particles, &
253 aug_pools=qmmm_env%aug_pools, &
254 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
255 coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
256 para_env=para_env, &
257 pw_pools=pw_pools, &
258 eps_mm_rspace=qmmm_env%eps_mm_rspace, &
259 cube_info=ks_qmmm_env_loc%cube_info, &
260 forces=forces, &
261 forces_added_charges=forces_added_charges, &
262 forces_added_shells=forces_added_shells, &
263 interp_section=interp_section, &
264 iw=iw, &
265 mm_cell=mm_cell)
266 CASE (do_qmmm_none)
267 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
268 "- No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
269 CASE DEFAULT
270 IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "- Unknown Coupling..."
271 cpabort("")
272 END SELECT
273 END IF
274 END IF
275 ! Correct Total Energy adding the contribution of the QM/MM nuclear interaction
276 energy%total = energy%total + energy%qmmm_nu
277 ! Proceed if gradients are requested..
278 IF (need_f) THEN
279 !ikuo Temporary change to alleviate compiler problems on Intel with
280 !array dimension of 0
281 IF (qmmm_env%num_mm_atoms /= 0) CALL para_env%sum(forces)
282 IF (qmmm_env%added_charges%num_mm_atoms /= 0) CALL para_env%sum(forces_added_charges)
283 IF (qmmm_env%added_shells%num_mm_atoms /= 0) CALL para_env%sum(forces_added_shells)
284 ! Debug Forces
285 IF (debug_this_module) THEN
286 IF (dft_control%qs_control%semi_empirical .OR. &
287 dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
288 WRITE (iw, *) "NO DEBUG AVAILABLE in module"//trim(routinen)
289 ELSE
290 ! Print Out Forces
291 IF (iw > 0) THEN
292 DO imm = 1, SIZE(qmmm_env%mm_atom_index)
293 WRITE (iw, *) "ANALYTICAL FORCES:"
294 indmm = qmmm_env%mm_atom_index(imm)
295 WRITE (iw, '(I6,3F15.9)') indmm, forces(:, imm)
296 END DO
297 END IF
298 CALL qmmm_debug_forces(rho=rho_tot_r, &
299 qs_env=qs_env, &
300 qmmm_env=qmmm_env, &
301 analytical_forces=forces, &
302 mm_particles=mm_particles, &
303 mm_atom_index=qmmm_env%mm_atom_index, &
304 num_mm_atoms=qmmm_env%num_mm_atoms, &
305 interp_section=interp_section, &
306 mm_cell=mm_cell)
307 END IF
308 END IF
309 END IF
310 ! Give back rho_tot_t to auxbas_pool only for GPW/GAPW
311 IF ((.NOT. dft_control%qs_control%semi_empirical) .AND. &
312 (.NOT. dft_control%qs_control%dftb) .AND. (.NOT. dft_control%qs_control%xtb)) THEN
313 CALL auxbas_pool%give_back_pw(rho_tot_r)
314 DEALLOCATE (rho_tot_r)
315 END IF
316 IF (iw > 0) THEN
317 IF (.NOT. gapw) WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
318 "QM/MM Nuclear Electrostatic Potential :", energy%qmmm_nu
319 WRITE (iw, '(T2,"QMMM|",1X,A,T66,F15.9)') &
320 "QMMM Total Energy (QM + QMMM electronic + QMMM nuclear):", energy%total
321 WRITE (iw, '(T2,"QMMM|",1X,A)') "MM energy NOT included in the above term!"// &
322 " Check for: FORCE_EVAL ( QMMM )"
323 WRITE (iw, '(T2,"QMMM|",1X,A)') "that includes both QM, QMMM and MM energy terms!"
324 END IF
325 IF (need_f) THEN
326 ! Transfer Forces
327 DO imm = 1, qmmm_env%num_mm_atoms
328 indmm = qmmm_env%mm_atom_index(imm)
329
330 !add image forces to Forces
331 IF (qmmm_env%image_charge) THEN
332 DO iatom = 1, qmmm_env%num_image_mm_atoms
333 image_indmm = qmmm_env%image_charge_pot%image_mm_list(iatom)
334 IF (image_indmm .EQ. indmm) THEN
335 forces(:, imm) = forces(:, imm) &
336 + qmmm_env%image_charge_pot%image_forcesMM(:, iatom)
337 END IF
338 END DO
339 END IF
340
341 ! Hack: In Forces there the gradients indeed...
342 ! Minux sign to take care of this misunderstanding...
343 mm_particles(indmm)%f(:) = -forces(:, imm) + mm_particles(indmm)%f(:)
344 END DO
345 DEALLOCATE (forces)
346 IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
347 DO imm = 1, qmmm_env%added_charges%num_mm_atoms
348 indmm = qmmm_env%added_charges%mm_atom_index(imm)
349 ! Hack: In Forces there the gradients indeed...
350 ! Minux sign to take care of this misunderstanding...
351 qmmm_env%added_charges%added_particles(indmm)%f(:) = -forces_added_charges(:, imm)
352 END DO
353 END IF
354 DEALLOCATE (forces_added_charges)
355 IF (qmmm_env%added_shells%num_mm_atoms .GT. 0) THEN
356 DO imm = 1, qmmm_env%added_shells%num_mm_atoms
357 indmm = qmmm_env%added_shells%mm_core_index(imm)
358 ! Hack: In Forces there the gradients indeed...
359 ! Minux sign to take care of this misunderstanding...
360 qmmm_env%added_shells%added_particles(imm)%f(:) = qmmm_env%added_shells%added_particles(imm)%f(:) - &
361 forces_added_shells(:, imm)
362
363 END DO
364 END IF
365 DEALLOCATE (forces_added_shells)
366 END IF
367 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
368 CALL timestop(handle)
369
370 END SUBROUTINE qmmm_forces
371
372! **************************************************************************************************
373!> \brief Evaluates the contribution to the forces due to the
374!> QM/MM potential computed collocating the Electrostatic
375!> Gaussian Potential.
376!> \param rho ...
377!> \param qmmm_env ...
378!> \param mm_particles ...
379!> \param aug_pools ...
380!> \param auxbas_grid ...
381!> \param coarser_grid ...
382!> \param cube_info ...
383!> \param para_env ...
384!> \param eps_mm_rspace ...
385!> \param pw_pools ...
386!> \param Forces ...
387!> \param Forces_added_charges ...
388!> \param Forces_added_shells ...
389!> \param interp_section ...
390!> \param iw ...
391!> \param mm_cell ...
392!> \par History
393!> 06.2004 created [tlaino]
394!> \author Teodoro Laino
395! **************************************************************************************************
396 SUBROUTINE qmmm_forces_with_gaussian(rho, qmmm_env, mm_particles, &
397 aug_pools, auxbas_grid, coarser_grid, cube_info, para_env, &
398 eps_mm_rspace, pw_pools, Forces, Forces_added_charges, Forces_added_shells, &
399 interp_section, iw, mm_cell)
400 TYPE(pw_r3d_rs_type), POINTER :: rho
401 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
402 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
403 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
404 INTEGER, INTENT(IN) :: auxbas_grid, coarser_grid
405 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
406 TYPE(mp_para_env_type), POINTER :: para_env
407 REAL(kind=dp), INTENT(IN) :: eps_mm_rspace
408 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
409 REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges, &
410 forces_added_shells
411 TYPE(section_vals_type), POINTER :: interp_section
412 INTEGER, INTENT(IN) :: iw
413 TYPE(cell_type), POINTER :: mm_cell
414
415 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_forces_with_gaussian'
416
417 INTEGER :: handle, i, igrid, j, k, kind_interp, me, &
418 ngrids, stat
419 INTEGER, DIMENSION(3) :: glb, gub, lb, ub
420 INTEGER, DIMENSION(:), POINTER :: pos_of_x
421 LOGICAL :: shells
422 REAL(kind=dp), DIMENSION(:, :), POINTER :: tmp
423 TYPE(mp_comm_type) :: group
424 TYPE(mp_request_type) :: request
425 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
426
427! Statements
428
429 CALL timeset(routinen, handle)
430 NULLIFY (tmp)
431 cpassert(ASSOCIATED(mm_particles))
432 cpassert(ASSOCIATED(qmmm_env%mm_atom_chrg))
433 cpassert(ASSOCIATED(qmmm_env%mm_atom_index))
434 cpassert(ASSOCIATED(forces))
435 !Statements
436 ngrids = SIZE(pw_pools)
437 CALL pw_pools_create_pws(aug_pools, grids)
438 DO igrid = 1, ngrids
439 CALL pw_zero(grids(igrid))
440 END DO
441 ! Collocate Density on multigrids
442 lb = rho%pw_grid%bounds_local(1, :)
443 ub = rho%pw_grid%bounds_local(2, :)
444 grids(auxbas_grid)%array(lb(1):ub(1), &
445 lb(2):ub(2), &
446 lb(3):ub(3)) = rho%array
447 ! copy the boundaries
448 DO i = lb(1), ub(1)
449 grids(auxbas_grid)%array(i, ub(2) + 1, ub(3) + 1) = rho%array(i, lb(2), lb(3))
450 END DO
451 DO k = lb(3), ub(3)
452 DO i = lb(1), ub(1)
453 grids(auxbas_grid)%array(i, ub(2) + 1, k) = rho%array(i, lb(2), k)
454 END DO
455 END DO
456 DO j = lb(2), ub(2)
457 DO i = lb(1), ub(1)
458 grids(auxbas_grid)%array(i, j, ub(3) + 1) = rho%array(i, j, lb(3))
459 END DO
460 END DO
461 pos_of_x => grids(auxbas_grid)%pw_grid%para%pos_of_x
462 group = grids(auxbas_grid)%pw_grid%para%group
463 me = grids(auxbas_grid)%pw_grid%para%group%mepos
464 glb = rho%pw_grid%bounds(1, :)
465 gub = rho%pw_grid%bounds(2, :)
466 IF ((pos_of_x(glb(1)) .EQ. me) .AND. (pos_of_x(gub(1)) .EQ. me)) THEN
467 DO k = lb(3), ub(3)
468 DO j = lb(2), ub(2)
469 grids(auxbas_grid)%array(ub(1) + 1, j, k) = rho%array(lb(1), j, k)
470 END DO
471 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, k) = rho%array(lb(1), lb(2), k)
472 END DO
473 DO j = lb(2), ub(2)
474 grids(auxbas_grid)%array(ub(1) + 1, j, ub(3) + 1) = rho%array(lb(1), j, lb(3))
475 END DO
476 grids(auxbas_grid)%array(ub(1) + 1, ub(2) + 1, ub(3) + 1) = rho%array(lb(1), lb(2), lb(3))
477 ELSE IF (pos_of_x(glb(1)) .EQ. me) THEN
478 ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
479 rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)), &
480 stat=stat)
481 cpassert(stat == 0)
482 tmp = rho%array(lb(1), :, :)
483 CALL group%isend(msgin=tmp, dest=pos_of_x(rho%pw_grid%bounds(2, 1)), &
484 request=request, tag=112)
485 CALL request%wait()
486 ELSE IF (pos_of_x(gub(1)) .EQ. me) THEN
487 ALLOCATE (tmp(rho%pw_grid%bounds_local(1, 2):rho%pw_grid%bounds_local(2, 2), &
488 rho%pw_grid%bounds_local(1, 3):rho%pw_grid%bounds_local(2, 3)), &
489 stat=stat)
490 cpassert(stat == 0)
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("")
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 .GT. 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
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
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) &
1216 ra(:) = pbc(mm_particles(lindmm)%r - dommoqm, mm_cell) + dommoqm
1217 qt = mm_charges(lindmm)
1218 ! Possible Spherical Cutoff
1219 IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
1220 CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
1221 qt = qt*sph_chrg_factor
1222 END IF
1223 IF (abs(qt) <= epsilon(0.0_dp)) cycle atoms
1224 rt1 = ra(1)
1225 rt2 = ra(2)
1226 rt3 = ra(3)
1227 ft1 = 0.0_dp
1228 ft2 = 0.0_dp
1229 ft3 = 0.0_dp
1230 loopongrid: DO k = bo(1, 3), bo(2, 3)
1231 my_k = k - gbo(1, 3)
1232 xs3 = real(my_k, dp)*dr3
1233 my_j = bo(1, 2) - gbo(1, 2)
1234 xs2 = real(my_j, dp)*dr2
1235 rv3 = rt3 - xs3
1236 DO j = bo(1, 2), bo(2, 2)
1237 my_i = bo(1, 1) - gbo(1, 1)
1238 xs1 = real(my_i, dp)*dr1
1239 rv2 = rt2 - xs2
1240 DO i = bo(1, 1), bo(2, 1)
1241 rv1 = rt1 - xs1
1242 r2 = rv1*rv1 + rv2*rv2 + rv3*rv3
1243 r = sqrt(r2)
1244 ix = floor(r/dx) + 1
1245 rx = (r - real(ix - 1, dp)*dx)/dx
1246 rx2 = rx*rx
1247 term = pot0_2(1, ix)*(-6.0_dp*(rx - rx2)) &
1248 + pot0_2(2, ix)*(1.0_dp - 4.0_dp*rx + 3.0_dp*rx2) &
1249 + pot0_2(1, ix + 1)*(6.0_dp*(rx - rx2)) &
1250 + pot0_2(2, ix + 1)*(-2.0_dp*rx + 3.0_dp*rx2)
1251 fac = grid(i, j, k)*term
1252 IF (r == 0.0_dp) THEN
1253 rd1 = 1.0_dp
1254 rd2 = 1.0_dp
1255 rd3 = 1.0_dp
1256 ELSE
1257 rd1 = rv1/r
1258 rd2 = rv2/r
1259 rd3 = rv3/r
1260 END IF
1261 ft1 = ft1 + fac*rd1
1262 ft2 = ft2 + fac*rd2
1263 ft3 = ft3 + fac*rd3
1264 xs1 = xs1 + dr1
1265 END DO
1266 xs2 = xs2 + dr2
1267 END DO
1268 END DO loopongrid
1269 qt = -qt*dvol/dx
1270 lforces(1, lindmm) = ft1*qt
1271 lforces(2, lindmm) = ft2*qt
1272 lforces(3, lindmm) = ft3*qt
1273
1274 forces(1, lindmm) = forces(1, lindmm) + lforces(1, lindmm)
1275 forces(2, lindmm) = forces(2, lindmm) + lforces(2, lindmm)
1276 forces(3, lindmm) = forces(3, lindmm) + lforces(3, lindmm)
1277 END DO atoms
1278 !$OMP END PARALLEL DO
1279 END DO radius
1280 !
1281 ! Debug Statement
1282 !
1283 IF (debug_this_module) THEN
1284 CALL debug_qmmm_forces_with_gauss_lr(pgfs=pgfs, &
1285 aug_pools=aug_pools, &
1286 rho=cgrid, &
1287 num_mm_atoms=num_mm_atoms, &
1288 mm_charges=mm_charges, &
1289 mm_atom_index=mm_atom_index, &
1290 mm_particles=mm_particles, &
1291 coarser_grid_level=coarser_grid_level, &
1292 debug_force=lforces, &
1293 potentials=potentials, &
1294 para_env=para_env, &
1295 mm_cell=mm_cell, &
1296 dommoqm=dommoqm, &
1297 iw=iw, &
1298 par_scheme=par_scheme, &
1299 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1300 shells=shells)
1301 END IF
1302
1303 DEALLOCATE (lforces)
1304 CALL timestop(handle)
1305 END SUBROUTINE qmmm_forces_with_gaussian_lr
1306
1307! **************************************************************************************************
1308!> \brief Evaluates numerically QM/MM forces and compares them with
1309!> the analytically computed ones.
1310!> It is evaluated only when debug_this_module is set to .TRUE.
1311!> \param rho ...
1312!> \param qs_env ...
1313!> \param qmmm_env ...
1314!> \param Analytical_Forces ...
1315!> \param mm_particles ...
1316!> \param mm_atom_index ...
1317!> \param num_mm_atoms ...
1318!> \param interp_section ...
1319!> \param mm_cell ...
1320!> \par History
1321!> 08.2004 created [tlaino]
1322!> \author Teodoro Laino
1323! **************************************************************************************************
1324 SUBROUTINE qmmm_debug_forces(rho, qs_env, qmmm_env, Analytical_Forces, &
1325 mm_particles, mm_atom_index, num_mm_atoms, &
1326 interp_section, mm_cell)
1327 TYPE(pw_r3d_rs_type), POINTER :: rho
1328 TYPE(qs_environment_type), POINTER :: qs_env
1329 TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1330 REAL(kind=dp), DIMENSION(:, :), POINTER :: analytical_forces
1331 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1332 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1333 INTEGER, INTENT(IN) :: num_mm_atoms
1334 TYPE(section_vals_type), POINTER :: interp_section
1335 TYPE(cell_type), POINTER :: mm_cell
1336
1337 CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_debug_forces'
1338
1339 INTEGER :: handle, i, indmm, iw, j, k
1340 REAL(kind=dp) :: coord_save
1341 REAL(kind=dp), DIMENSION(2) :: energy
1342 REAL(kind=dp), DIMENSION(3) :: err
1343 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1344 TYPE(cp_logger_type), POINTER :: logger
1345 TYPE(mp_para_env_type), POINTER :: para_env
1346 TYPE(pw_env_type), POINTER :: pw_env
1347 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1348 TYPE(pw_r3d_rs_type) :: v_qmmm_rspace
1349 TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
1350 TYPE(section_vals_type), POINTER :: input_section, print_section
1351
1352 CALL timeset(routinen, handle)
1353 NULLIFY (num_forces)
1354 CALL get_qs_env(qs_env=qs_env, &
1355 pw_env=pw_env, &
1356 input=input_section, &
1357 para_env=para_env)
1358
1359 print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
1360 logger => cp_get_default_logger()
1361 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".qmmmLog")
1362 CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools)
1363 CALL pw_pools(1)%pool%create_pw(v_qmmm_rspace)
1364 ALLOCATE (num_forces(3, num_mm_atoms))
1365 ks_qmmm_env_loc => qs_env%ks_qmmm_env
1366 IF (iw > 0) WRITE (iw, '(/A)') "DEBUG SECTION:"
1367 atoms: DO i = 1, num_mm_atoms
1368 indmm = mm_atom_index(i)
1369 coords: DO j = 1, 3
1370 coord_save = mm_particles(indmm)%r(j)
1371 energy = 0.0_dp
1372 diff: DO k = 1, 2
1373 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1374 CALL pw_zero(v_qmmm_rspace)
1375 SELECT CASE (qmmm_env%qmmm_coupl_type)
1376 CASE (do_qmmm_coulomb)
1377 cpabort("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1378 CASE (do_qmmm_pcharge)
1379 cpabort("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
1380 CASE (do_qmmm_gauss, do_qmmm_swave)
1381 CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env, &
1382 v_qmmm=v_qmmm_rspace, &
1383 mm_particles=mm_particles, &
1384 aug_pools=qmmm_env%aug_pools, &
1385 para_env=para_env, &
1386 eps_mm_rspace=qmmm_env%eps_mm_rspace, &
1387 cube_info=ks_qmmm_env_loc%cube_info, &
1388 pw_pools=pw_pools, &
1389 auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
1390 coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
1391 interp_section=interp_section, &
1392 mm_cell=mm_cell)
1393 CASE (do_qmmm_none)
1394 cycle diff
1395 CASE DEFAULT
1396 IF (iw > 0) WRITE (iw, '(T3,A)') "Unknown Coupling..."
1397 cpabort("")
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) &
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 cpassert(abs(err(1)) <= maxerr)
1430 cpassert(abs(err(2)) <= maxerr)
1431 cpassert(abs(err(3)) <= maxerr)
1432 END DO
1433 CASE (do_qmmm_none)
1434 IF (iw > 0) WRITE (iw, '(T3,A)') "No QM/MM Derivatives to debug. Just Mechanical Coupling!"
1435 CASE DEFAULT
1436 IF (iw > 0) WRITE (iw, '(T3,A)') "Unknown Coupling..."
1437 cpabort("")
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) &
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 cpassert(abs(err(1)) <= maxerr)
1542 cpassert(abs(err(2)) <= maxerr)
1543 cpassert(abs(err(3)) <= maxerr)
1544
1545 IF (ASSOCIATED(xdat)) THEN
1546 DEALLOCATE (xdat)
1547 END IF
1548 IF (ASSOCIATED(ydat)) THEN
1549 DEALLOCATE (ydat)
1550 END IF
1551 IF (ASSOCIATED(zdat)) THEN
1552 DEALLOCATE (zdat)
1553 END IF
1554
1555 CALL pw_pools_give_back_pws(aug_pools, grids)
1556 CALL timestop(handle)
1557100 FORMAT("Collocation : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1558 END SUBROUTINE debug_integrate_gf_rspace_nopbc
1559
1560! **************************************************************************************************
1561!> \brief Debugs qmmm_forces_with_gaussian_LG.. It may helps too ... ;-]
1562!> \param pgfs ...
1563!> \param aug_pools ...
1564!> \param rho ...
1565!> \param mm_charges ...
1566!> \param mm_atom_index ...
1567!> \param mm_particles ...
1568!> \param num_mm_atoms ...
1569!> \param coarser_grid_level ...
1570!> \param per_potentials ...
1571!> \param debug_force ...
1572!> \param para_env ...
1573!> \param mm_cell ...
1574!> \param dOmmOqm ...
1575!> \param iw ...
1576!> \param par_scheme ...
1577!> \param qmmm_spherical_cutoff ...
1578!> \param shells ...
1579!> \par History
1580!> 08.2004 created [tlaino]
1581!> \author Teodoro Laino
1582! **************************************************************************************************
1583 SUBROUTINE debug_qmmm_forces_with_gauss_lg(pgfs, aug_pools, rho, mm_charges, mm_atom_index, &
1584 mm_particles, num_mm_atoms, coarser_grid_level, per_potentials, &
1585 debug_force, para_env, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1586
1587 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1588 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1589 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1590 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1591 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1592 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1593 INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level
1594 TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
1595 REAL(kind=dp), DIMENSION(:, :) :: debug_force
1596 TYPE(mp_para_env_type), POINTER :: para_env
1597 TYPE(cell_type), POINTER :: mm_cell
1598 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1599 INTEGER, INTENT(IN) :: iw, par_scheme
1600 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1601 LOGICAL :: shells
1602
1603 CHARACTER(len=*), PARAMETER :: routinen = 'debug_qmmm_forces_with_gauss_LG'
1604
1605 INTEGER :: handle, i, igrid, indmm, j, k, ngrids
1606 REAL(kind=dp) :: coord_save
1607 REAL(kind=dp), DIMENSION(2) :: energy
1608 REAL(kind=dp), DIMENSION(3) :: err
1609 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1610 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1611
1612 ALLOCATE (num_forces(3, num_mm_atoms))
1613 CALL timeset(routinen, handle)
1614 ngrids = SIZE(aug_pools)
1615 CALL pw_pools_create_pws(aug_pools, grids)
1616 DO igrid = 1, ngrids
1617 CALL pw_zero(grids(igrid))
1618 END DO
1619 atoms: DO i = 1, num_mm_atoms
1620 indmm = mm_atom_index(i)
1621 coords: DO j = 1, 3
1622 coord_save = mm_particles(indmm)%r(j)
1623 energy = 0.0_dp
1624 diff: DO k = 1, 2
1625 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1626 CALL pw_zero(grids(coarser_grid_level))
1627
1628 CALL qmmm_elec_with_gaussian_lg(pgfs=pgfs, &
1629 cgrid=grids(coarser_grid_level), &
1630 mm_charges=mm_charges, &
1631 mm_atom_index=mm_atom_index, &
1632 mm_particles=mm_particles, &
1633 para_env=para_env, &
1634 per_potentials=per_potentials, &
1635 mm_cell=mm_cell, &
1636 dommoqm=dommoqm, &
1637 par_scheme=par_scheme, &
1638 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1639 shells=shells)
1640
1641 energy(k) = pw_integral_ab(rho, grids(coarser_grid_level))
1642 END DO diff
1643 IF (iw > 0) &
1644 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1645 "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1646 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1647 mm_particles(indmm)%r(j) = coord_save
1648 END DO coords
1649 END DO atoms
1650
1651 DO i = 1, num_mm_atoms
1652 indmm = mm_atom_index(i)
1653 err = 0.0_dp
1654 IF (all(num_forces /= 0.0_dp)) THEN
1655 err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1656 err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1657 err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1658 END IF
1659 IF (iw > 0) &
1660 WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1661 debug_force(2, i), num_forces(2, i), err(2), &
1662 debug_force(3, i), num_forces(3, i), err(3)
1663 cpassert(abs(err(1)) <= maxerr)
1664 cpassert(abs(err(2)) <= maxerr)
1665 cpassert(abs(err(3)) <= maxerr)
1666 END DO
1667
1668 DEALLOCATE (num_forces)
1669 CALL pw_pools_give_back_pws(aug_pools, grids)
1670 CALL timestop(handle)
1671100 FORMAT("MM Atom LR : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1672 END SUBROUTINE debug_qmmm_forces_with_gauss_lg
1673
1674! **************************************************************************************************
1675!> \brief Debugs qmmm_forces_with_gaussian_LR.. It may helps too ... ;-]
1676!> \param pgfs ...
1677!> \param aug_pools ...
1678!> \param rho ...
1679!> \param mm_charges ...
1680!> \param mm_atom_index ...
1681!> \param mm_particles ...
1682!> \param num_mm_atoms ...
1683!> \param coarser_grid_level ...
1684!> \param potentials ...
1685!> \param debug_force ...
1686!> \param para_env ...
1687!> \param mm_cell ...
1688!> \param dOmmOqm ...
1689!> \param iw ...
1690!> \param par_scheme ...
1691!> \param qmmm_spherical_cutoff ...
1692!> \param shells ...
1693!> \par History
1694!> 08.2004 created [tlaino]
1695!> \author Teodoro Laino
1696! **************************************************************************************************
1697 SUBROUTINE debug_qmmm_forces_with_gauss_lr(pgfs, aug_pools, rho, mm_charges, mm_atom_index, &
1698 mm_particles, num_mm_atoms, coarser_grid_level, potentials, &
1699 debug_force, para_env, mm_cell, dOmmOqm, iw, par_scheme, qmmm_spherical_cutoff, shells)
1700
1701 TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
1702 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: aug_pools
1703 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
1704 REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
1705 INTEGER, DIMENSION(:), POINTER :: mm_atom_index
1706 TYPE(particle_type), DIMENSION(:), POINTER :: mm_particles
1707 INTEGER, INTENT(IN) :: num_mm_atoms, coarser_grid_level
1708 TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
1709 REAL(kind=dp), DIMENSION(:, :) :: debug_force
1710 TYPE(mp_para_env_type), POINTER :: para_env
1711 TYPE(cell_type), POINTER :: mm_cell
1712 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: dommoqm
1713 INTEGER, INTENT(IN) :: iw, par_scheme
1714 REAL(kind=dp), DIMENSION(2), INTENT(IN) :: qmmm_spherical_cutoff
1715 LOGICAL :: shells
1716
1717 CHARACTER(len=*), PARAMETER :: routinen = 'debug_qmmm_forces_with_gauss_LR'
1718
1719 INTEGER :: handle, i, igrid, indmm, j, k, ngrids
1720 REAL(kind=dp) :: coord_save
1721 REAL(kind=dp), DIMENSION(2) :: energy
1722 REAL(kind=dp), DIMENSION(3) :: err
1723 REAL(kind=dp), DIMENSION(:, :), POINTER :: num_forces
1724 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: grids
1725
1726 ALLOCATE (num_forces(3, num_mm_atoms))
1727 CALL timeset(routinen, handle)
1728 ngrids = SIZE(aug_pools)
1729 CALL pw_pools_create_pws(aug_pools, grids)
1730 DO igrid = 1, ngrids
1731 CALL pw_zero(grids(igrid))
1732 END DO
1733 atoms: DO i = 1, num_mm_atoms
1734 indmm = mm_atom_index(i)
1735 coords: DO j = 1, 3
1736 coord_save = mm_particles(indmm)%r(j)
1737 energy = 0.0_dp
1738 diff: DO k = 1, 2
1739 mm_particles(indmm)%r(j) = coord_save + (-1)**k*dx
1740 CALL pw_zero(grids(coarser_grid_level))
1741
1742 CALL qmmm_elec_with_gaussian_lr(pgfs=pgfs, &
1743 grid=grids(coarser_grid_level), &
1744 mm_charges=mm_charges, &
1745 mm_atom_index=mm_atom_index, &
1746 mm_particles=mm_particles, &
1747 para_env=para_env, &
1748 potentials=potentials, &
1749 mm_cell=mm_cell, &
1750 dommoqm=dommoqm, &
1751 par_scheme=par_scheme, &
1752 qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
1753 shells=shells)
1754
1755 energy(k) = pw_integral_ab(rho, grids(coarser_grid_level))
1756 END DO diff
1757 IF (iw > 0) &
1758 WRITE (iw, '(A,I6,A,I3,A,2F15.9)') &
1759 "DEBUG LR:: MM Atom = ", indmm, " Coord = ", j, " Energies (+/-) :: ", energy(2), energy(1)
1760 num_forces(j, i) = (energy(2) - energy(1))/(2.0_dp*dx)
1761 mm_particles(indmm)%r(j) = coord_save
1762 END DO coords
1763 END DO atoms
1764
1765 DO i = 1, num_mm_atoms
1766 indmm = mm_atom_index(i)
1767 err = 0.0_dp
1768 IF (all(num_forces(:, i) /= 0.0_dp)) THEN
1769 err(1) = (debug_force(1, i) - num_forces(1, i))/num_forces(1, i)*100.0_dp
1770 err(2) = (debug_force(2, i) - num_forces(2, i))/num_forces(2, i)*100.0_dp
1771 err(3) = (debug_force(3, i) - num_forces(3, i))/num_forces(3, i)*100.0_dp
1772 END IF
1773 IF (iw > 0) &
1774 WRITE (iw, 100) indmm, debug_force(1, i), num_forces(1, i), err(1), &
1775 debug_force(2, i), num_forces(2, i), err(2), &
1776 debug_force(3, i), num_forces(3, i), err(3)
1777 cpassert(abs(err(1)) <= maxerr)
1778 cpassert(abs(err(2)) <= maxerr)
1779 cpassert(abs(err(3)) <= maxerr)
1780 END DO
1781
1782 DEALLOCATE (num_forces)
1783 CALL pw_pools_give_back_pws(aug_pools, grids)
1784 CALL timestop(handle)
1785100 FORMAT("MM Atom LR : ", i5, 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ", 2f15.9, " ( ", f7.2, " ) ")
1786 END SUBROUTINE debug_qmmm_forces_with_gauss_lr
1787
1788END MODULE qmmm_gpw_forces
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:51
cell_types::pbc
Definition cell_types.F:92
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_integral_ab
Definition pw_methods.F:221
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
pw_pool_types::pw_pools_create_pws
Definition pw_pool_types.F:139
pw_pool_types::pw_pools_give_back_pws
Definition pw_pool_types.F:150
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_spline_utils
utils to manipulate splines on the regular grid of a pw
Definition cp_spline_utils.F:14
cp_spline_utils::spline3_nopbc_interp
integer, parameter, public spline3_nopbc_interp
Definition cp_spline_utils.F:38
cp_spline_utils::spline3_pbc_interp
integer, parameter, public spline3_pbc_interp
Definition cp_spline_utils.F:38
cp_spline_utils::pw_restrict_s3
subroutine, public pw_restrict_s3(pw_fine_in, pw_coarse_out, coarse_pool, param_section)
restricts the function from a fine grid to a coarse one
Definition cp_spline_utils.F:55
cube_utils
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
Definition cube_utils.F:18
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_par_atom
integer, parameter, public do_par_atom
Definition input_constants.F:287
input_constants::do_qmmm_none
integer, parameter, public do_qmmm_none
Definition input_constants.F:385
input_constants::do_qmmm_pcharge
integer, parameter, public do_qmmm_pcharge
Definition input_constants.F:385
input_constants::do_qmmm_coulomb
integer, parameter, public do_qmmm_coulomb
Definition input_constants.F:385
input_constants::do_qmmm_swave
integer, parameter, public do_qmmm_swave
Definition input_constants.F:385
input_constants::do_qmmm_gauss
integer, parameter, public do_qmmm_gauss
Definition input_constants.F:385
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:153
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mm_collocate_potential
Calculate the MM potential by collocating the primitive Gaussian functions (pgf)
Definition mm_collocate_potential.F:15
mm_collocate_potential::integrate_gf_rspace_nopbc
subroutine, public integrate_gf_rspace_nopbc(zetp, rp, scale, w, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo, force, n_rep_real, mm_cell)
Main driver to integrate gaussian functions on a grid function without using periodic boundary condit...
Definition mm_collocate_potential.F:296
mm_collocate_potential::collocate_gf_rspace_nopbc
subroutine, public collocate_gf_rspace_nopbc(zetp, rp, scale, w, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo2, n_rep_real, mm_cell)
Main driver to collocate gaussian functions on grid without using periodic boundary conditions (NoPBC...
Definition mm_collocate_potential.F:173
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qmmm_gaussian_types
Sets the typo for the gaussian treatment of the qm/mm interaction.
Definition qmmm_gaussian_types.F:14
qmmm_gpw_energy
A collection of methods to treat the QM/MM electrostatic coupling.
Definition qmmm_gpw_energy.F:14
qmmm_gpw_energy::qmmm_elec_with_gaussian
subroutine, public qmmm_elec_with_gaussian(qmmm_env, v_qmmm, mm_particles, aug_pools, cube_info, para_env, eps_mm_rspace, pw_pools, auxbas_grid, coarser_grid, interp_section, mm_cell)
Compute the QM/MM electrostatic Interaction collocating the gaussian Electrostatic Potential.
Definition qmmm_gpw_energy.F:243
qmmm_gpw_energy::qmmm_elec_with_gaussian_lr
subroutine, public qmmm_elec_with_gaussian_lr(pgfs, grid, mm_charges, mm_atom_index, mm_particles, para_env, potentials, mm_cell, dommoqm, par_scheme, qmmm_spherical_cutoff, shells)
Compute the QM/MM electrostatic Interaction collocating (1/R - Sum_NG Gaussians) on the coarser grid ...
Definition qmmm_gpw_energy.F:757
qmmm_gpw_energy::qmmm_elec_with_gaussian_lg
subroutine, public qmmm_elec_with_gaussian_lg(pgfs, cgrid, mm_charges, mm_atom_index, mm_particles, para_env, per_potentials, mm_cell, dommoqm, par_scheme, qmmm_spherical_cutoff, shells)
Compute the QM/MM electrostatic Interaction collocating (1/R - Sum_NG Gaussians) on the coarser grid ...
Definition qmmm_gpw_energy.F:532
qmmm_gpw_forces
Routines to compute energy and forces in a QM/MM calculation.
Definition qmmm_gpw_forces.F:14
qmmm_gpw_forces::qmmm_forces
subroutine, public qmmm_forces(qs_env, qmmm_env, mm_particles, calc_force, mm_cell)
General driver to Compute the contribution to the forces due to the QM/MM potential.
Definition qmmm_gpw_forces.F:100
qmmm_se_forces
Calculation of the derivative of the QMMM Hamiltonian integral matrix <a|\sum_i q_i|b> for semi-empir...
Definition qmmm_se_forces.F:13
qmmm_se_forces::deriv_se_qmmm_matrix
subroutine, public deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian QMMM terms.
Definition qmmm_se_forces.F:78
qmmm_tb_methods
TB methods used with QMMM.
Definition qmmm_tb_methods.F:12
qmmm_tb_methods::deriv_tb_qmmm_matrix
subroutine, public deriv_tb_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition qmmm_tb_methods.F:593
qmmm_tb_methods::deriv_tb_qmmm_matrix_pc
subroutine, public deriv_tb_qmmm_matrix_pc(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, forces, forces_added_charges)
Constructs the derivative w.r.t. 1-el DFTB hamiltonian QMMM terms.
Definition qmmm_tb_methods.F:830
qmmm_types_low
Definition qmmm_types_low.F:13
qmmm_util
Definition qmmm_util.F:13
qmmm_util::spherical_cutoff_factor
subroutine, public spherical_cutoff_factor(spherical_cutoff, rij, factor)
Computes a spherical cutoff factor for the QMMM interactions.
Definition qmmm_util.F:615
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_ks_qmmm_types
Definition qs_ks_qmmm_types.F:13
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cube_utils::cube_info_type
Definition cube_utils.F:38
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qmmm_gaussian_types::qmmm_gaussian_p_type
represent a pointer to a qmmm_gaussian_type, to be able to create arrays of pointers
Definition qmmm_gaussian_types.F:42
qmmm_gaussian_types::qmmm_gaussian_type
Definition qmmm_gaussian_types.F:27
qmmm_types_low::qmmm_env_qm_type
...
Definition qmmm_types_low.F:195
qmmm_types_low::qmmm_per_pot_p_type
Definition qmmm_types_low.F:100
qmmm_types_low::qmmm_per_pot_type
Periodic Potential.
Definition qmmm_types_low.F:90
qmmm_types_low::qmmm_pot_p_type
Definition qmmm_types_low.F:83
qmmm_types_low::qmmm_pot_type
Real Space Potential.
Definition qmmm_types_low.F:76
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
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