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