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