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