(git:b279b6b)
qmmm_image_charge.F
Go to the documentation of this file.
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 for image charge calculation within QM/MM
10 !> \par History
11 !> 12.2011 created
12 !> \author Dorothea Golze
13 ! **************************************************************************************************
14 MODULE qmmm_image_charge
15  USE ao_util, ONLY: exp_radius_very_extended
16  USE cell_types, ONLY: cell_type,&
17  pbc
18  USE cp_control_types, ONLY: dft_control_type
19  USE cp_eri_mme_interface, ONLY: cp_eri_mme_param,&
20  cp_eri_mme_update_local_counts
21  USE cp_files, ONLY: close_file,&
22  open_file
23  USE cp_log_handling, ONLY: cp_get_default_logger,&
24  cp_logger_type
25  USE cp_output_handling, ONLY: cp_p_file,&
26  cp_print_key_finished_output,&
27  cp_print_key_generate_filename,&
28  cp_print_key_should_output,&
29  cp_print_key_unit_nr
30  USE eri_mme_integrate, ONLY: eri_mme_2c_integrate
31  USE input_constants, ONLY: calc_always,&
32  calc_once,&
33  calc_once_done,&
34  do_eri_gpw,&
35  do_eri_mme
36  USE input_section_types, ONLY: section_vals_get_subs_vals,&
37  section_vals_type,&
38  section_vals_val_get
39  USE kinds, ONLY: default_path_length,&
40  dp
41  USE mathconstants, ONLY: pi
42  USE memory_utilities, ONLY: reallocate
43  USE message_passing, ONLY: mp_para_env_type
44  USE pw_env_types, ONLY: pw_env_get,&
45  pw_env_type
46  USE pw_methods, ONLY: pw_axpy,&
47  pw_integral_ab,&
48  pw_scale,&
49  pw_transfer,&
50  pw_zero
51  USE pw_poisson_methods, ONLY: pw_poisson_solve
52  USE pw_poisson_types, ONLY: pw_poisson_type
53  USE pw_pool_types, ONLY: pw_pool_type
54  USE pw_types, ONLY: pw_c1d_gs_type,&
55  pw_r3d_rs_type
56  USE qmmm_types_low, ONLY: qmmm_env_qm_type
57  USE qs_collocate_density, ONLY: calculate_rho_metal,&
58  calculate_rho_single_gaussian
59  USE qs_energy_types, ONLY: qs_energy_type
60  USE qs_environment_types, ONLY: get_qs_env,&
61  qs_environment_type
62  USE qs_integrate_potential, ONLY: integrate_pgf_product
63  USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
64  realspace_grid_type,&
65  transfer_pw2rs
66  USE util, ONLY: get_limit
67  USE virial_types, ONLY: virial_type
68 #include "./base/base_uses.f90"
69 
70  IMPLICIT NONE
71  PRIVATE
72 
73  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
74  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_image_charge'
75 
76  PUBLIC :: calculate_image_pot, &
77  integrate_potential_devga_rspace, &
78  conditional_calc_image_matrix, &
79  add_image_pot_to_hartree_pot, &
80  print_image_coefficients
81 
82 !***
83 CONTAINS
84 ! **************************************************************************************************
85 !> \brief determines coefficients by solving image_matrix*coeff=-pot_const by
86 !> Gaussian elimination or in an iterative fashion and calculates
87 !> image/metal potential with these coefficients
88 !> \param v_hartree_rspace Hartree potential in real space
89 !> \param rho_hartree_gspace Kohn Sham density in reciprocal space
90 !> \param energy structure where energies are stored
91 !> \param qmmm_env qmmm environment
92 !> \param qs_env qs environment
93 ! **************************************************************************************************
94  SUBROUTINE calculate_image_pot(v_hartree_rspace, rho_hartree_gspace, energy, &
95  qmmm_env, qs_env)
96 
97  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree_rspace
98  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_hartree_gspace
99  TYPE(qs_energy_type), POINTER :: energy
100  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
101  TYPE(qs_environment_type), POINTER :: qs_env
102 
103  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_image_pot'
104 
105  INTEGER :: handle
106 
107  CALL timeset(routinen, handle)
108 
109  IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
110  !calculate preconditioner matrix for CG if necessary
111  IF (qs_env%calc_image_preconditioner) THEN
112  IF (qmmm_env%image_charge_pot%image_restart) THEN
113  CALL restart_image_matrix(image_matrix=qs_env%image_matrix, &
114  qs_env=qs_env, qmmm_env=qmmm_env)
115  ELSE
116  CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
117  qs_env=qs_env, qmmm_env=qmmm_env)
118  END IF
119  END IF
120  CALL calc_image_coeff_iterative(v_hartree_rspace=v_hartree_rspace, &
121  coeff=qs_env%image_coeff, qmmm_env=qmmm_env, &
122  qs_env=qs_env)
123 
124  ELSE
125  CALL calc_image_coeff_gaussalgorithm(v_hartree_rspace=v_hartree_rspace, &
126  coeff=qs_env%image_coeff, qmmm_env=qmmm_env, &
127  qs_env=qs_env)
128  END IF
129 
130  ! calculate the image/metal potential with the optimized coefficients
131  ALLOCATE (qs_env%ks_qmmm_env%v_metal_rspace)
132  CALL calculate_potential_metal(v_metal_rspace= &
133  qs_env%ks_qmmm_env%v_metal_rspace, coeff=qs_env%image_coeff, &
134  rho_hartree_gspace=rho_hartree_gspace, &
135  energy=energy, qs_env=qs_env)
136 
137  CALL timestop(handle)
138 
139  END SUBROUTINE calculate_image_pot
140 
141 ! **************************************************************************************************
142 !> \brief determines coefficients by solving the linear set of equations
143 !> image_matrix*coeff=-pot_const using a Gaussian elimination scheme
144 !> \param v_hartree_rspace Hartree potential in real space
145 !> \param coeff expansion coefficients of the image charge density, i.e.
146 !> rho_metal=sum_a c_a*g_a
147 !> \param qmmm_env qmmm environment
148 !> \param qs_env qs environment
149 ! **************************************************************************************************
150  SUBROUTINE calc_image_coeff_gaussalgorithm(v_hartree_rspace, coeff, qmmm_env, &
151  qs_env)
152 
153  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree_rspace
154  REAL(kind=dp), DIMENSION(:), POINTER :: coeff
155  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
156  TYPE(qs_environment_type), POINTER :: qs_env
157 
158  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_image_coeff_gaussalgorithm'
159 
160  INTEGER :: handle, info, natom
161  REAL(kind=dp) :: eta, v0
162  REAL(kind=dp), DIMENSION(:), POINTER :: pot_const
163 
164  CALL timeset(routinen, handle)
165 
166  NULLIFY (pot_const)
167 
168  !minus sign V0: account for the fact that v_hartree has the opposite sign
169  v0 = -qmmm_env%image_charge_pot%V0
170  eta = qmmm_env%image_charge_pot%eta
171  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
172 
173  ALLOCATE (pot_const(natom))
174  IF (.NOT. ASSOCIATED(coeff)) THEN
175  ALLOCATE (coeff(natom))
176  END IF
177  coeff = 0.0_dp
178 
179  CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
180  pot_const)
181  !add integral V0*ga(r)
182  pot_const(:) = -pot_const(:) + v0*sqrt((pi/eta)**3)
183 
184  !solve linear system of equations T*coeff=-pot_const
185  !LU factorization of T by DGETRF done in calculate_image_matrix
186  CALL dgetrs('N', natom, 1, qs_env%image_matrix, natom, qs_env%ipiv, &
187  pot_const, natom, info)
188  cpassert(info == 0)
189 
190  coeff = pot_const
191 
192  DEALLOCATE (pot_const)
193 
194  CALL timestop(handle)
195 
196  END SUBROUTINE calc_image_coeff_gaussalgorithm
197 
198 ! **************************************************************************************************
199 !> \brief determines image coefficients iteratively
200 !> \param v_hartree_rspace Hartree potential in real space
201 !> \param coeff expansion coefficients of the image charge density, i.e.
202 !> rho_metal=sum_a c_a*g_a
203 !> \param qmmm_env qmmm environment
204 !> \param qs_env qs environment
205 ! **************************************************************************************************
206  SUBROUTINE calc_image_coeff_iterative(v_hartree_rspace, coeff, qmmm_env, &
207  qs_env)
208 
209  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree_rspace
210  REAL(kind=dp), DIMENSION(:), POINTER :: coeff
211  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
212  TYPE(qs_environment_type), POINTER :: qs_env
213 
214  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_image_coeff_iterative'
215 
216  INTEGER :: handle, iter_steps, natom, output_unit
217  REAL(kind=dp) :: alpha, eta, rsnew, rsold, v0
218  REAL(kind=dp), DIMENSION(:), POINTER :: ad, d, pot_const, r, vmetal_const, z
219  TYPE(cp_logger_type), POINTER :: logger
220  TYPE(pw_r3d_rs_type) :: auxpot_ad_rspace, v_metal_rspace_guess
221  TYPE(section_vals_type), POINTER :: input
222 
223  CALL timeset(routinen, handle)
224 
225  NULLIFY (pot_const, vmetal_const, logger, input)
226  logger => cp_get_default_logger()
227 
228  !minus sign V0: account for the fact that v_hartree has the opposite sign
229  v0 = -qmmm_env%image_charge_pot%V0
230  eta = qmmm_env%image_charge_pot%eta
231  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
232 
233  ALLOCATE (pot_const(natom))
234  ALLOCATE (vmetal_const(natom))
235  ALLOCATE (r(natom))
236  ALLOCATE (d(natom))
237  ALLOCATE (z(natom))
238  ALLOCATE (ad(natom))
239  IF (.NOT. ASSOCIATED(coeff)) THEN
240  ALLOCATE (coeff(natom))
241  END IF
242 
243  CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
244  pot_const)
245 
246  !add integral V0*ga(r)
247  pot_const(:) = -pot_const(:) + v0*sqrt((pi/eta)**3)
248 
249  !initial guess for coeff
250  coeff = 1.0_dp
251  d = 0.0_dp
252  z = 0.0_dp
253  r = 0.0_dp
254  rsold = 0.0_dp
255  rsnew = 0.0_dp
256  iter_steps = 0
257 
258  !calculate first guess of image/metal potential
259  CALL calculate_potential_metal(v_metal_rspace=v_metal_rspace_guess, &
260  coeff=coeff, qs_env=qs_env)
261  CALL integrate_potential_ga_rspace(potential=v_metal_rspace_guess, &
262  qmmm_env=qmmm_env, qs_env=qs_env, int_res=vmetal_const)
263 
264  ! modify coefficients iteratively
265  r = pot_const - vmetal_const
266  z = matmul(qs_env%image_matrix, r)
267  d = z
268  rsold = dot_product(r, z)
269 
270  DO
271  !calculate A*d
272  ad = 0.0_dp
273  CALL calculate_potential_metal(v_metal_rspace= &
274  auxpot_ad_rspace, coeff=d, qs_env=qs_env)
275  CALL integrate_potential_ga_rspace(potential= &
276  auxpot_ad_rspace, qmmm_env=qmmm_env, &
277  qs_env=qs_env, int_res=ad)
278 
279  alpha = rsold/dot_product(d, ad)
280  coeff = coeff + alpha*d
281 
282  r = r - alpha*ad
283  z = matmul(qs_env%image_matrix, r)
284  rsnew = dot_product(r, z)
285  iter_steps = iter_steps + 1
286  ! SQRT(rsnew) < 1.0E-08
287  IF (rsnew < 1.0e-16) THEN
288  CALL auxpot_ad_rspace%release()
289  EXIT
290  END IF
291  d = z + rsnew/rsold*d
292  rsold = rsnew
293  CALL auxpot_ad_rspace%release()
294  END DO
295 
296  ! print iteration info
297  CALL get_qs_env(qs_env=qs_env, &
298  input=input)
299  output_unit = cp_print_key_unit_nr(logger, input, &
300  "QMMM%PRINT%PROGRAM_RUN_INFO", &
301  extension=".qmmmLog")
302  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T3,A,T74,I7)") &
303  "Number of iteration steps for determination of image coefficients:", iter_steps
304  CALL cp_print_key_finished_output(output_unit, logger, input, &
305  "QMMM%PRINT%PROGRAM_RUN_INFO")
306 
307  IF (iter_steps .LT. 25) THEN
308  qs_env%calc_image_preconditioner = .false.
309  ELSE
310  qs_env%calc_image_preconditioner = .true.
311  END IF
312 
313  CALL v_metal_rspace_guess%release()
314  DEALLOCATE (pot_const)
315  DEALLOCATE (vmetal_const)
316  DEALLOCATE (r)
317  DEALLOCATE (d, z)
318  DEALLOCATE (ad)
319 
320  CALL timestop(handle)
321 
322  END SUBROUTINE calc_image_coeff_iterative
323 
324 ! ****************************************************************************
325 !> \brief calculates the integral V(r)*ga(r)
326 !> \param potential any potential
327 !> \param qmmm_env qmmm environment
328 !> \param qs_env qs environment
329 !> \param int_res result of the integration
330 !> \param atom_num atom index, needed when calculating image_matrix
331 !> \param atom_num_ref index of reference atom, needed when calculating
332 !> image_matrix
333 ! **************************************************************************************************
334  SUBROUTINE integrate_potential_ga_rspace(potential, qmmm_env, qs_env, int_res, &
335  atom_num, atom_num_ref)
336 
337  TYPE(pw_r3d_rs_type), INTENT(IN) :: potential
338  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
339  TYPE(qs_environment_type), POINTER :: qs_env
340  REAL(kind=dp), DIMENSION(:), POINTER :: int_res
341  INTEGER, INTENT(IN), OPTIONAL :: atom_num, atom_num_ref
342 
343  CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_potential_ga_rspace'
344 
345  INTEGER :: atom_a, atom_b, atom_ref, handle, iatom, &
346  j, k, natom, npme
347  INTEGER, DIMENSION(:), POINTER :: cores
348  REAL(kind=dp) :: eps_rho_rspace, radius
349  REAL(kind=dp), DIMENSION(3) :: ra
350  REAL(kind=dp), DIMENSION(:, :), POINTER :: hab
351  TYPE(cell_type), POINTER :: cell
352  TYPE(dft_control_type), POINTER :: dft_control
353  TYPE(mp_para_env_type), POINTER :: para_env
354  TYPE(pw_env_type), POINTER :: pw_env
355  TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
356  TYPE(realspace_grid_type), POINTER :: rs_v
357 
358  CALL timeset(routinen, handle)
359 
360  NULLIFY (cores, hab, cell, auxbas_rs_desc, pw_env, para_env, &
361  dft_control, rs_v)
362  ALLOCATE (hab(1, 1))
363 
364  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
365  CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
366  auxbas_rs_grid=rs_v)
367  CALL transfer_pw2rs(rs_v, potential)
368 
369  CALL get_qs_env(qs_env=qs_env, &
370  cell=cell, &
371  dft_control=dft_control, &
372  para_env=para_env, pw_env=pw_env)
373 
374  eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
375 
376  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
377  k = 1
378  IF (PRESENT(atom_num)) k = atom_num
379 
380  CALL reallocate(cores, 1, natom - k + 1)
381  int_res = 0.0_dp
382  npme = 0
383  cores = 0
384 
385  DO iatom = k, natom
386  IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
387  ! replicated realspace grid, split the atoms up between procs
388  IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
389  npme = npme + 1
390  cores(npme) = iatom
391  END IF
392  ELSE
393  npme = npme + 1
394  cores(npme) = iatom
395  END IF
396  END DO
397 
398  DO j = 1, npme
399 
400  iatom = cores(j)
401  atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
402 
403  IF (PRESENT(atom_num) .AND. PRESENT(atom_num_ref)) THEN
404  ! shift the function since potential only calculate for ref atom
405  atom_b = qmmm_env%image_charge_pot%image_mm_list(k)
406  atom_ref = qmmm_env%image_charge_pot%image_mm_list(atom_num_ref)
407  ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell) &
408  - pbc(qmmm_env%image_charge_pot%particles_all(atom_b)%r, cell) &
409  + pbc(qmmm_env%image_charge_pot%particles_all(atom_ref)%r, cell)
410 
411  ELSE
412  ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
413  END IF
414 
415  hab(1, 1) = 0.0_dp
416 
417  radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
418  ra=ra, rb=ra, rp=ra, &
419  zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
420  prefactor=1.0_dp, cutoff=1.0_dp)
421 
422  CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
423  0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
424  rs_v, hab, o1=0, o2=0, &
425  radius=radius, calculate_forces=.false., &
426  use_subpatch=.true., subpatch_pattern=0)
427 
428  int_res(iatom) = hab(1, 1)
429 
430  END DO
431 
432  CALL para_env%sum(int_res)
433 
434  DEALLOCATE (hab, cores)
435 
436  CALL timestop(handle)
437 
438  END SUBROUTINE integrate_potential_ga_rspace
439 
440 ! **************************************************************************************************
441 !> \brief calculates the image forces on the MM atoms
442 !> \param potential any potential, in this case: Hartree potential
443 !> \param coeff expansion coefficients of the image charge density, i.e.
444 !> rho_metal=sum_a c_a*g_a
445 !> \param forces structure storing the force contribution of the image charges
446 !> for the metal (MM) atoms
447 !> \param qmmm_env qmmm environment
448 !> \param qs_env qs environment
449 ! **************************************************************************************************
450  SUBROUTINE integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, &
451  qs_env)
452 
453  TYPE(pw_r3d_rs_type), INTENT(IN) :: potential
454  REAL(kind=dp), DIMENSION(:), POINTER :: coeff
455  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
456  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
457  TYPE(qs_environment_type), POINTER :: qs_env
458 
459  CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_potential_devga_rspace'
460 
461  INTEGER :: atom_a, handle, iatom, j, natom, npme
462  INTEGER, DIMENSION(:), POINTER :: cores
463  LOGICAL :: use_virial
464  REAL(kind=dp) :: eps_rho_rspace, radius
465  REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
466  REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab
467  TYPE(cell_type), POINTER :: cell
468  TYPE(dft_control_type), POINTER :: dft_control
469  TYPE(mp_para_env_type), POINTER :: para_env
470  TYPE(pw_env_type), POINTER :: pw_env
471  TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
472  TYPE(realspace_grid_type), POINTER :: rs_v
473  TYPE(virial_type), POINTER :: virial
474 
475  CALL timeset(routinen, handle)
476 
477  NULLIFY (cores, hab, pab, cell, auxbas_rs_desc, pw_env, para_env, &
478  dft_control, rs_v, virial)
479  use_virial = .false.
480 
481  ALLOCATE (hab(1, 1))
482  ALLOCATE (pab(1, 1))
483 
484  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
485  CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
486  auxbas_rs_grid=rs_v)
487  CALL transfer_pw2rs(rs_v, potential)
488 
489  CALL get_qs_env(qs_env=qs_env, &
490  cell=cell, &
491  dft_control=dft_control, &
492  para_env=para_env, pw_env=pw_env, &
493  virial=virial)
494 
495  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
496 
497  IF (use_virial) THEN
498  cpabort("Virial not implemented for image charge method")
499  END IF
500 
501  eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
502 
503  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
504 
505  IF (.NOT. ASSOCIATED(forces)) THEN
506  ALLOCATE (forces(3, natom))
507  END IF
508 
509  forces(:, :) = 0.0_dp
510 
511  CALL reallocate(cores, 1, natom)
512  npme = 0
513  cores = 0
514 
515  DO iatom = 1, natom
516  IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
517  ! replicated realspace grid, split the atoms up between procs
518  IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
519  npme = npme + 1
520  cores(npme) = iatom
521  END IF
522  ELSE
523  npme = npme + 1
524  cores(npme) = iatom
525  END IF
526  END DO
527 
528  DO j = 1, npme
529 
530  iatom = cores(j)
531  atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
532  ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
533  hab(1, 1) = 0.0_dp
534  pab(1, 1) = 1.0_dp
535  force_a(:) = 0.0_dp
536  force_b(:) = 0.0_dp
537 
538  radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
539  ra=ra, rb=ra, rp=ra, &
540  zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
541  pab=pab, o1=0, o2=0, & ! without map_consistent
542  prefactor=1.0_dp, cutoff=1.0_dp)
543 
544  CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
545  0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
546  rs_v, hab, pab, o1=0, o2=0, &
547  radius=radius, calculate_forces=.true., &
548  force_a=force_a, force_b=force_b, use_subpatch=.true., &
549  subpatch_pattern=0)
550 
551  force_a(:) = coeff(iatom)*force_a(:)
552  forces(:, iatom) = force_a(:)
553 
554  END DO
555 
556  CALL para_env%sum(forces)
557 
558  DEALLOCATE (hab, pab, cores)
559 
560  ! print info on gradients if wanted
561  CALL print_gradients_image_atoms(forces, qs_env)
562 
563  CALL timestop(handle)
564 
565  END SUBROUTINE integrate_potential_devga_rspace
566 
567 !****************************************************************************
568 !> \brief calculate image matrix T depending on constraints on image atoms
569 !> in case coefficients are estimated not iteratively
570 !> \param qs_env qs environment
571 !> \param qmmm_env qmmm environment
572 ! **************************************************************************************************
573  SUBROUTINE conditional_calc_image_matrix(qs_env, qmmm_env)
574 
575  TYPE(qs_environment_type), POINTER :: qs_env
576  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
577 
578  IF (.NOT. qmmm_env%image_charge_pot%coeff_iterative) THEN
579  SELECT CASE (qmmm_env%image_charge_pot%state_image_matrix)
580  CASE (calc_always)
581  CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
582  ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
583  CASE (calc_once)
584  !if all image atoms are fully constrained, calculate image matrix
585  !only for the first MD or GEO_OPT step
586  CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
587  ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
588  qmmm_env%image_charge_pot%state_image_matrix = calc_once_done
589  IF (qmmm_env%center_qm_subsys0) &
590  CALL cp_warn(__location__, &
591  "The image atoms are fully "// &
592  "constrained and the image matrix is only calculated once. "// &
593  "To be safe, set CENTER to NEVER ")
594  CASE (calc_once_done)
595  ! do nothing image matrix is stored
596  CASE DEFAULT
597  cpabort("No initialization for image charges available?")
598  END SELECT
599  END IF
600 
601  END SUBROUTINE conditional_calc_image_matrix
602 
603 !****************************************************************************
604 !> \brief calculate image matrix T
605 !> \param image_matrix matrix T
606 !> \param ipiv pivoting prior to DGETRS (for Gaussian elimination)
607 !> \param qs_env qs environment
608 !> \param qmmm_env qmmm environment
609 ! **************************************************************************************************
610  SUBROUTINE calculate_image_matrix(image_matrix, ipiv, qs_env, qmmm_env)
611 
612  REAL(kind=dp), DIMENSION(:, :), POINTER :: image_matrix
613  INTEGER, DIMENSION(:), OPTIONAL, POINTER :: ipiv
614  TYPE(qs_environment_type), POINTER :: qs_env
615  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
616 
617  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_image_matrix'
618 
619  INTEGER :: handle, j, k, natom, output_unit, stat
620  TYPE(cp_logger_type), POINTER :: logger
621  TYPE(section_vals_type), POINTER :: input
622 
623  CALL timeset(routinen, handle)
624  NULLIFY (input, logger)
625 
626  logger => cp_get_default_logger()
627 
628  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
629 
630  IF (.NOT. ASSOCIATED(image_matrix)) THEN
631  ALLOCATE (image_matrix(natom, natom))
632  END IF
633  IF (PRESENT(ipiv)) THEN
634  IF (.NOT. ASSOCIATED(ipiv)) THEN
635  ALLOCATE (ipiv(natom))
636  END IF
637  ipiv = 0
638  END IF
639 
640  CALL get_qs_env(qs_env, input=input)
641  !print info
642  output_unit = cp_print_key_unit_nr(logger, input, &
643  "QMMM%PRINT%PROGRAM_RUN_INFO", &
644  extension=".qmmmLog")
645  IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
646  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T3,A)") &
647  "Calculating image matrix"
648  ELSE
649  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T2,A)") &
650  "Calculating image matrix"
651  END IF
652  CALL cp_print_key_finished_output(output_unit, logger, input, &
653  "QMMM%PRINT%PROGRAM_RUN_INFO")
654 
655  ! Calculate image matrix using either GPW or MME method
656  SELECT CASE (qmmm_env%image_charge_pot%image_matrix_method)
657  CASE (do_eri_gpw)
658  CALL calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
659  CASE (do_eri_mme)
660  CALL calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
661  CASE DEFAULT
662  cpabort("Unknown method for calculating image matrix")
663  END SELECT
664 
665  IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
666  !inversion --> preconditioner matrix for CG
667  CALL dpotrf('L', natom, qs_env%image_matrix, natom, stat)
668  cpassert(stat == 0)
669  CALL dpotri('L', natom, qs_env%image_matrix, natom, stat)
670  cpassert(stat == 0)
671  DO j = 1, natom
672  DO k = j + 1, natom
673  qs_env%image_matrix(j, k) = qs_env%image_matrix(k, j)
674  END DO
675  END DO
676  CALL write_image_matrix(qs_env%image_matrix, qs_env)
677  ELSE
678  !pivoting prior to DGETRS (Gaussian elimination)
679  IF (PRESENT(ipiv)) THEN
680  CALL dgetrf(natom, natom, image_matrix, natom, ipiv, stat)
681  cpassert(stat == 0)
682  END IF
683  END IF
684 
685  CALL timestop(handle)
686 
687  END SUBROUTINE calculate_image_matrix
688 
689 ! **************************************************************************************************
690 !> \brief calculate image matrix T using GPW method
691 !> \param image_matrix matrix T
692 !> \param qs_env qs environment
693 !> \param qmmm_env qmmm environment
694 ! **************************************************************************************************
695  SUBROUTINE calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
696  REAL(kind=dp), DIMENSION(:, :), POINTER :: image_matrix
697  TYPE(qs_environment_type), POINTER :: qs_env
698  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
699 
700  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_image_matrix_gpw'
701 
702  INTEGER :: handle, iatom, iatom_ref, natom
703  REAL(kind=dp), DIMENSION(:), POINTER :: int_res
704  TYPE(mp_para_env_type), POINTER :: para_env
705  TYPE(pw_c1d_gs_type) :: rho_gb, vb_gspace
706  TYPE(pw_env_type), POINTER :: pw_env
707  TYPE(pw_poisson_type), POINTER :: poisson_env
708  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
709  TYPE(pw_r3d_rs_type) :: vb_rspace
710 
711  CALL timeset(routinen, handle)
712  NULLIFY (pw_env, auxbas_pw_pool, poisson_env, para_env, int_res)
713 
714  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
715  ALLOCATE (int_res(natom))
716 
717  image_matrix = 0.0_dp
718 
719  CALL get_qs_env(qs_env, pw_env=pw_env, para_env=para_env)
720 
721  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
722  poisson_env=poisson_env)
723  CALL auxbas_pw_pool%create_pw(rho_gb)
724  CALL auxbas_pw_pool%create_pw(vb_gspace)
725  CALL auxbas_pw_pool%create_pw(vb_rspace)
726 
727  ! calculate vb only once for one reference atom
728  iatom_ref = 1 !
729  !collocate gaussian of reference MM atom on grid
730  CALL pw_zero(rho_gb)
731  CALL calculate_rho_single_gaussian(rho_gb, qs_env, iatom_ref)
732  !calculate potential vb like hartree potential
733  CALL pw_zero(vb_gspace)
734  CALL pw_poisson_solve(poisson_env, rho_gb, vhartree=vb_gspace)
735  CALL pw_zero(vb_rspace)
736  CALL pw_transfer(vb_gspace, vb_rspace)
737  CALL pw_scale(vb_rspace, vb_rspace%pw_grid%dvol)
738 
739  DO iatom = 1, natom
740  !calculate integral vb_rspace*ga
741  int_res = 0.0_dp
742  CALL integrate_potential_ga_rspace(vb_rspace, qs_env%qmmm_env_qm, &
743  qs_env, int_res, atom_num=iatom, &
744  atom_num_ref=iatom_ref)
745  image_matrix(iatom, iatom:natom) = int_res(iatom:natom)
746  image_matrix(iatom + 1:natom, iatom) = int_res(iatom + 1:natom)
747  END DO
748 
749  CALL vb_gspace%release()
750  CALL vb_rspace%release()
751  CALL rho_gb%release()
752 
753  DEALLOCATE (int_res)
754 
755  CALL timestop(handle)
756  END SUBROUTINE calculate_image_matrix_gpw
757 
758 ! **************************************************************************************************
759 !> \brief calculate image matrix T using MME (MiniMax-Ewald) method
760 !> \param image_matrix matrix T
761 !> \param qs_env qs environment
762 !> \param qmmm_env qmmm environment
763 ! **************************************************************************************************
764  SUBROUTINE calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
765  REAL(kind=dp), DIMENSION(:, :), POINTER :: image_matrix
766  TYPE(qs_environment_type), POINTER :: qs_env
767  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
768 
769  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_image_matrix_mme'
770 
771  INTEGER :: atom_a, handle, iatom, natom
772  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: zeta
773  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ra
774  TYPE(mp_para_env_type), POINTER :: para_env
775 
776  CALL timeset(routinen, handle)
777  NULLIFY (para_env)
778 
779  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
780  ALLOCATE (zeta(natom), ra(3, natom))
781 
782  zeta(:) = qmmm_env%image_charge_pot%eta
783 
784  DO iatom = 1, natom
785  atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
786  ra(:, iatom) = qmmm_env%image_charge_pot%particles_all(atom_a)%r(:)
787  END DO
788 
789  CALL get_qs_env(qs_env, para_env=para_env)
790 
791  CALL integrate_s_mme(qmmm_env%image_charge_pot%eri_mme_param, &
792  zeta, zeta, ra, ra, image_matrix, para_env)
793 
794  CALL timestop(handle)
795  END SUBROUTINE calculate_image_matrix_mme
796 
797 ! **************************************************************************************************
798 !> \brief high-level integration routine for 2c integrals over s-type functions.
799 !> Parallelization over pairs of functions.
800 !> \param param ...
801 !> \param zeta ...
802 !> \param zetb ...
803 !> \param ra ...
804 !> \param rb ...
805 !> \param hab ...
806 !> \param para_env ...
807 ! **************************************************************************************************
808  SUBROUTINE integrate_s_mme(param, zeta, zetb, ra, rb, hab, para_env)
809  TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
810  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: zeta, zetb
811  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
812  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: hab
813  TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
814 
815  CHARACTER(len=*), PARAMETER :: routinen = 'integrate_s_mme'
816 
817  INTEGER :: g_count, handle, ipgf, ipgf_prod, jpgf, &
818  npgf_prod, npgfa, npgfb, r_count
819  INTEGER, DIMENSION(2) :: limits
820  REAL(kind=dp), DIMENSION(3) :: rab
821 
822  CALL timeset(routinen, handle)
823  g_count = 0; r_count = 0
824 
825  hab(:, :) = 0.0_dp
826 
827  npgfa = SIZE(zeta)
828  npgfb = SIZE(zetb)
829  npgf_prod = npgfa*npgfb ! total number of integrals
830 
831  limits = get_limit(npgf_prod, para_env%num_pe, para_env%mepos)
832 
833  DO ipgf_prod = limits(1), limits(2)
834  ipgf = (ipgf_prod - 1)/npgfb + 1
835  jpgf = mod(ipgf_prod - 1, npgfb) + 1
836  rab(:) = ra(:, ipgf) - rb(:, jpgf)
837  CALL eri_mme_2c_integrate(param%par, 0, 0, 0, 0, zeta(ipgf), &
838  zetb(jpgf), rab, hab, ipgf - 1, jpgf - 1, g_count=g_count, r_count=r_count)
839  END DO
840 
841  CALL cp_eri_mme_update_local_counts(param, para_env, g_count_2c=g_count, r_count_2c=r_count)
842  CALL para_env%sum(hab)
843  CALL timestop(handle)
844 
845  END SUBROUTINE integrate_s_mme
846 
847 ! **************************************************************************************************
848 !> \brief calculates potential of the metal (image potential) given a set of
849 !> coefficients coeff
850 !> \param v_metal_rspace potential generated by rho_metal in real space
851 !> \param coeff expansion coefficients of the image charge density, i.e.
852 !> rho_metal=sum_a c_a*g_a
853 !> \param rho_hartree_gspace Kohn Sham density in reciprocal space
854 !> \param energy structure where energies are stored
855 !> \param qs_env qs environment
856 ! **************************************************************************************************
857  SUBROUTINE calculate_potential_metal(v_metal_rspace, coeff, rho_hartree_gspace, energy, &
858  qs_env)
859 
860  TYPE(pw_r3d_rs_type), INTENT(OUT) :: v_metal_rspace
861  REAL(kind=dp), DIMENSION(:), POINTER :: coeff
862  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_hartree_gspace
863  TYPE(qs_energy_type), OPTIONAL, POINTER :: energy
864  TYPE(qs_environment_type), POINTER :: qs_env
865 
866  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_potential_metal'
867 
868  INTEGER :: handle
869  REAL(kind=dp) :: en_external, en_vmetal_rhohartree, &
870  total_rho_metal
871  TYPE(pw_c1d_gs_type) :: rho_metal, v_metal_gspace
872  TYPE(pw_env_type), POINTER :: pw_env
873  TYPE(pw_poisson_type), POINTER :: poisson_env
874  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
875 
876  CALL timeset(routinen, handle)
877 
878  NULLIFY (pw_env, auxbas_pw_pool, poisson_env)
879  en_vmetal_rhohartree = 0.0_dp
880  en_external = 0.0_dp
881 
882  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
883  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
884  poisson_env=poisson_env)
885 
886  CALL auxbas_pw_pool%create_pw(rho_metal)
887 
888  CALL auxbas_pw_pool%create_pw(v_metal_gspace)
889 
890  CALL auxbas_pw_pool%create_pw(v_metal_rspace)
891 
892  CALL pw_zero(rho_metal)
893  CALL calculate_rho_metal(rho_metal, coeff, total_rho_metal=total_rho_metal, &
894  qs_env=qs_env)
895 
896  CALL pw_zero(v_metal_gspace)
897  CALL pw_poisson_solve(poisson_env, rho_metal, &
898  vhartree=v_metal_gspace)
899 
900  IF (PRESENT(rho_hartree_gspace)) THEN
901  en_vmetal_rhohartree = 0.5_dp*pw_integral_ab(v_metal_gspace, &
902  rho_hartree_gspace)
903  en_external = qs_env%qmmm_env_qm%image_charge_pot%V0*total_rho_metal
904  energy%image_charge = en_vmetal_rhohartree - 0.5_dp*en_external
905  CALL print_image_energy_terms(en_vmetal_rhohartree, en_external, &
906  total_rho_metal, qs_env)
907  END IF
908 
909  CALL pw_zero(v_metal_rspace)
910  CALL pw_transfer(v_metal_gspace, v_metal_rspace)
911  CALL pw_scale(v_metal_rspace, v_metal_rspace%pw_grid%dvol)
912  CALL v_metal_gspace%release()
913  CALL rho_metal%release()
914 
915  CALL timestop(handle)
916 
917  END SUBROUTINE calculate_potential_metal
918 
919 ! ****************************************************************************
920 !> \brief Add potential of metal (image charge pot) to Hartree Potential
921 !> \param v_hartree Hartree potential (in real space)
922 !> \param v_metal potential generated by rho_metal (in real space)
923 !> \param qs_env qs environment
924 ! **************************************************************************************************
925  SUBROUTINE add_image_pot_to_hartree_pot(v_hartree, v_metal, qs_env)
926 
927  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_hartree
928  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_metal
929  TYPE(qs_environment_type), POINTER :: qs_env
930 
931  CHARACTER(len=*), PARAMETER :: routinen = 'add_image_pot_to_hartree_pot'
932 
933  INTEGER :: handle, output_unit
934  TYPE(cp_logger_type), POINTER :: logger
935  TYPE(section_vals_type), POINTER :: input
936 
937  CALL timeset(routinen, handle)
938 
939  NULLIFY (input, logger)
940  logger => cp_get_default_logger()
941 
942  !add image charge potential
943  CALL pw_axpy(v_metal, v_hartree)
944 
945  ! print info
946  CALL get_qs_env(qs_env=qs_env, &
947  input=input)
948  output_unit = cp_print_key_unit_nr(logger, input, &
949  "QMMM%PRINT%PROGRAM_RUN_INFO", &
950  extension=".qmmmLog")
951  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T3,A)") &
952  "Adding image charge potential to the Hartree potential."
953  CALL cp_print_key_finished_output(output_unit, logger, input, &
954  "QMMM%PRINT%PROGRAM_RUN_INFO")
955 
956  CALL timestop(handle)
957 
958  END SUBROUTINE add_image_pot_to_hartree_pot
959 
960 !****************************************************************************
961 !> \brief writes image matrix T to file when used as preconditioner for
962 !> calculating image coefficients iteratively
963 !> \param image_matrix matrix T
964 !> \param qs_env qs environment
965 ! **************************************************************************************************
966  SUBROUTINE write_image_matrix(image_matrix, qs_env)
967 
968  REAL(kind=dp), DIMENSION(:, :), POINTER :: image_matrix
969  TYPE(qs_environment_type), POINTER :: qs_env
970 
971  CHARACTER(LEN=*), PARAMETER :: routinen = 'write_image_matrix'
972 
973  CHARACTER(LEN=default_path_length) :: filename
974  INTEGER :: handle, rst_unit
975  TYPE(cp_logger_type), POINTER :: logger
976  TYPE(mp_para_env_type), POINTER :: para_env
977  TYPE(section_vals_type), POINTER :: print_key, qmmm_section
978 
979  CALL timeset(routinen, handle)
980 
981  NULLIFY (qmmm_section, print_key, logger, para_env)
982  logger => cp_get_default_logger()
983  rst_unit = -1
984 
985  CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
986  input=qmmm_section)
987 
988  print_key => section_vals_get_subs_vals(qmmm_section, &
989  "QMMM%PRINT%IMAGE_CHARGE_RESTART")
990 
991  IF (btest(cp_print_key_should_output(logger%iter_info, &
992  qmmm_section, "QMMM%PRINT%IMAGE_CHARGE_RESTART"), &
993  cp_p_file)) THEN
994 
995  rst_unit = cp_print_key_unit_nr(logger, qmmm_section, &
996  "QMMM%PRINT%IMAGE_CHARGE_RESTART", &
997  extension=".Image", &
998  file_status="REPLACE", &
999  file_action="WRITE", &
1000  file_form="UNFORMATTED")
1001 
1002  IF (rst_unit > 0) filename = cp_print_key_generate_filename(logger, &
1003  print_key, extension=".IMAGE", &
1004  my_local=.false.)
1005 
1006  IF (rst_unit > 0) THEN
1007  WRITE (rst_unit) image_matrix
1008  END IF
1009 
1010  CALL cp_print_key_finished_output(rst_unit, logger, qmmm_section, &
1011  "QMMM%PRINT%IMAGE_CHARGE_RESTART")
1012  END IF
1013 
1014  CALL timestop(handle)
1015 
1016  END SUBROUTINE write_image_matrix
1017 
1018 !****************************************************************************
1019 !> \brief restarts image matrix T when used as preconditioner for calculating
1020 !> image coefficients iteratively
1021 !> \param image_matrix matrix T
1022 !> \param qs_env qs environment
1023 !> \param qmmm_env qmmm environment
1024 ! **************************************************************************************************
1025  SUBROUTINE restart_image_matrix(image_matrix, qs_env, qmmm_env)
1026 
1027  REAL(kind=dp), DIMENSION(:, :), POINTER :: image_matrix
1028  TYPE(qs_environment_type), POINTER :: qs_env
1029  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1030 
1031  CHARACTER(LEN=*), PARAMETER :: routinen = 'restart_image_matrix'
1032 
1033  CHARACTER(LEN=default_path_length) :: image_filename
1034  INTEGER :: handle, natom, output_unit, rst_unit
1035  LOGICAL :: exist
1036  TYPE(cp_logger_type), POINTER :: logger
1037  TYPE(mp_para_env_type), POINTER :: para_env
1038  TYPE(section_vals_type), POINTER :: qmmm_section
1039 
1040  CALL timeset(routinen, handle)
1041 
1042  NULLIFY (qmmm_section, logger, para_env)
1043  logger => cp_get_default_logger()
1044  exist = .false.
1045  rst_unit = -1
1046 
1047  natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
1048 
1049  IF (.NOT. ASSOCIATED(image_matrix)) THEN
1050  ALLOCATE (image_matrix(natom, natom))
1051  END IF
1052 
1053  image_matrix = 0.0_dp
1054 
1055  CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
1056  input=qmmm_section)
1057 
1058  CALL section_vals_val_get(qmmm_section, "QMMM%IMAGE_CHARGE%IMAGE_RESTART_FILE_NAME", &
1059  c_val=image_filename)
1060 
1061  INQUIRE (file=image_filename, exist=exist)
1062 
1063  IF (exist) THEN
1064  IF (para_env%is_source()) THEN
1065  CALL open_file(file_name=image_filename, &
1066  file_status="OLD", &
1067  file_form="UNFORMATTED", &
1068  file_action="READ", &
1069  unit_number=rst_unit)
1070 
1071  READ (rst_unit) qs_env%image_matrix
1072  END IF
1073 
1074  CALL para_env%bcast(qs_env%image_matrix)
1075 
1076  IF (para_env%is_source()) CALL close_file(unit_number=rst_unit)
1077 
1078  output_unit = cp_print_key_unit_nr(logger, qmmm_section, &
1079  "QMMM%PRINT%PROGRAM_RUN_INFO", &
1080  extension=".qmmmLog")
1081  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T3,A)") &
1082  "Restarted image matrix"
1083  ELSE
1084  cpabort("Restart file for image matrix not found")
1085  END IF
1086 
1087  qmmm_env%image_charge_pot%image_restart = .false.
1088 
1089  CALL timestop(handle)
1090 
1091  END SUBROUTINE restart_image_matrix
1092 
1093 ! ****************************************************************************
1094 !> \brief Print info on image gradients on image MM atoms
1095 !> \param forces structure storing the force contribution of the image charges
1096 !> for the metal (MM) atoms (actually these are only the gradients)
1097 !> \param qs_env qs environment
1098 ! **************************************************************************************************
1099  SUBROUTINE print_gradients_image_atoms(forces, qs_env)
1100 
1101  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
1102  TYPE(qs_environment_type), POINTER :: qs_env
1103 
1104  INTEGER :: atom_a, iatom, natom, output_unit
1105  REAL(kind=dp), DIMENSION(3) :: sum_gradients
1106  TYPE(cp_logger_type), POINTER :: logger
1107  TYPE(section_vals_type), POINTER :: input
1108 
1109  NULLIFY (input, logger)
1110  logger => cp_get_default_logger()
1111 
1112  sum_gradients = 0.0_dp
1113  natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1114 
1115  DO iatom = 1, natom
1116  sum_gradients(:) = sum_gradients(:) + forces(:, iatom)
1117  END DO
1118 
1119  CALL get_qs_env(qs_env=qs_env, input=input)
1120 
1121  output_unit = cp_print_key_unit_nr(logger, input, &
1122  "QMMM%PRINT%DERIVATIVES", extension=".Log")
1123  IF (output_unit > 0) THEN
1124  WRITE (unit=output_unit, fmt="(/1X,A)") &
1125  "Image gradients [a.u.] on MM image charge atoms after QMMM calculation: "
1126  WRITE (unit=output_unit, fmt="(T4,A4,T27,A1,T50,A1,T74,A1)") &
1127  "Atom", "X", "Y", "Z"
1128  DO iatom = 1, natom
1129  atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1130  WRITE (unit=output_unit, fmt="(T2,I6,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
1131  atom_a, forces(:, iatom)
1132  END DO
1133 
1134  WRITE (unit=output_unit, fmt="(T2,A)") repeat("-", 79)
1135  WRITE (unit=output_unit, fmt="(T2,A,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
1136  "sum gradients:", sum_gradients
1137  WRITE (unit=output_unit, fmt="(/)")
1138  END IF
1139 
1140  CALL cp_print_key_finished_output(output_unit, logger, input, &
1141  "QMMM%PRINT%DERIVATIVES")
1142 
1143  END SUBROUTINE print_gradients_image_atoms
1144 
1145 ! ****************************************************************************
1146 !> \brief Print image coefficients
1147 !> \param image_coeff expansion coefficients of the image charge density
1148 !> \param qs_env qs environment
1149 ! **************************************************************************************************
1150  SUBROUTINE print_image_coefficients(image_coeff, qs_env)
1151 
1152  REAL(kind=dp), DIMENSION(:), POINTER :: image_coeff
1153  TYPE(qs_environment_type), POINTER :: qs_env
1154 
1155  INTEGER :: atom_a, iatom, natom, output_unit
1156  REAL(kind=dp) :: normalize_factor, sum_coeff
1157  TYPE(cp_logger_type), POINTER :: logger
1158  TYPE(section_vals_type), POINTER :: input
1159 
1160  NULLIFY (input, logger)
1161  logger => cp_get_default_logger()
1162 
1163  sum_coeff = 0.0_dp
1164  natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1165  normalize_factor = sqrt((qs_env%qmmm_env_qm%image_charge_pot%eta/pi)**3)
1166 
1167  DO iatom = 1, natom
1168  sum_coeff = sum_coeff + image_coeff(iatom)
1169  END DO
1170 
1171  CALL get_qs_env(qs_env=qs_env, input=input)
1172 
1173  output_unit = cp_print_key_unit_nr(logger, input, &
1174  "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=".Log")
1175  IF (output_unit > 0) THEN
1176  WRITE (unit=output_unit, fmt="(/)")
1177  WRITE (unit=output_unit, fmt="(T2,A)") &
1178  "Image charges [a.u.] after QMMM calculation: "
1179  WRITE (unit=output_unit, fmt="(T4,A4,T67,A)") "Atom", "Image charge"
1180  WRITE (unit=output_unit, fmt="(T4,A,T67,A)") repeat("-", 4), repeat("-", 12)
1181 
1182  DO iatom = 1, natom
1183  atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1184  !opposite sign for image_coeff; during the calculation they have
1185  !the 'wrong' sign to ensure consistency with v_hartree which has
1186  !the opposite sign
1187  WRITE (unit=output_unit, fmt="(T2,I6,T65,ES16.9)") &
1188  atom_a, -image_coeff(iatom)/normalize_factor
1189  END DO
1190 
1191  WRITE (unit=output_unit, fmt="(T2,A)") repeat("-", 79)
1192  WRITE (unit=output_unit, fmt="(T2,A,T65,ES16.9)") &
1193  "sum image charges:", -sum_coeff/normalize_factor
1194  END IF
1195 
1196  CALL cp_print_key_finished_output(output_unit, logger, input, &
1197  "QMMM%PRINT%IMAGE_CHARGE_INFO")
1198 
1199  END SUBROUTINE print_image_coefficients
1200 
1201 ! ****************************************************************************
1202 !> \brief Print detailed image charge energies
1203 !> \param en_vmetal_rhohartree energy contribution of the image charges
1204 !> without external potential, i.e. 0.5*integral(v_metal*rho_hartree)
1205 !> \param en_external additional energy contribution of the image charges due
1206 !> to an external potential, i.e. V0*total_rho_metal
1207 !> \param total_rho_metal total induced image charge density
1208 !> \param qs_env qs environment
1209 ! **************************************************************************************************
1210  SUBROUTINE print_image_energy_terms(en_vmetal_rhohartree, en_external, &
1211  total_rho_metal, qs_env)
1212 
1213  REAL(kind=dp), INTENT(IN) :: en_vmetal_rhohartree, en_external, &
1214  total_rho_metal
1215  TYPE(qs_environment_type), POINTER :: qs_env
1216 
1217  INTEGER :: output_unit
1218  TYPE(cp_logger_type), POINTER :: logger
1219  TYPE(section_vals_type), POINTER :: input
1220 
1221  NULLIFY (input, logger)
1222  logger => cp_get_default_logger()
1223 
1224  CALL get_qs_env(qs_env=qs_env, input=input)
1225 
1226  output_unit = cp_print_key_unit_nr(logger, input, &
1227  "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=".Log")
1228 
1229  IF (output_unit > 0) THEN
1230  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1231  "Total induced charge density [a.u.]:", total_rho_metal
1232  WRITE (unit=output_unit, fmt="(T3,A)") "Image energy terms: "
1233  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1234  "Coulomb energy of QM and image charge density [a.u.]:", en_vmetal_rhohartree
1235  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1236  "External potential energy term [a.u.]:", -0.5_dp*en_external
1237  WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1238  "Total image charge energy [a.u.]:", en_vmetal_rhohartree - 0.5_dp*en_external
1239  END IF
1240 
1241  CALL cp_print_key_finished_output(output_unit, logger, input, &
1242  "QMMM%PRINT%IMAGE_CHARGE_INFO")
1243 
1244  END SUBROUTINE print_image_energy_terms
1245 
1246 END MODULE qmmm_image_charge
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
ao_util
All kind of helpful little routines.
Definition: ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition: ao_util.F:209
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_eri_mme_interface
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
Definition: cp_eri_mme_interface.F:12
cp_eri_mme_interface::cp_eri_mme_update_local_counts
subroutine, public cp_eri_mme_update_local_counts(param, para_env, G_count_2c, R_count_2c, GG_count_3c, GR_count_3c, RR_count_3c)
Update local counters to gather statistics on different paths taken in MME algorithm (each Ewald sum ...
Definition: cp_eri_mme_interface.F:714
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition: cp_output_handling.F:704
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
eri_mme_integrate
Minimax-Ewald (MME) method for calculating 2-center and 3-center electron repulsion integrals (ERI) o...
Definition: eri_mme_integrate.F:23
eri_mme_integrate::eri_mme_2c_integrate
subroutine, public eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, hab, o1, o2, G_count, R_count, normalize, potential, pot_par)
Low-level integration routine for 2-center ERIs.
Definition: eri_mme_integrate.F:72
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_eri_mme
integer, parameter, public do_eri_mme
Definition: input_constants.F:1226
input_constants::calc_once_done
integer, parameter, public calc_once_done
Definition: input_constants.F:400
input_constants::calc_once
integer, parameter, public calc_once
Definition: input_constants.F:400
input_constants::calc_always
integer, parameter, public calc_always
Definition: input_constants.F:400
input_constants::do_eri_gpw
integer, parameter, public do_eri_gpw
Definition: input_constants.F:1226
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_poisson_methods
Definition: pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qmmm_image_charge
Routines for image charge calculation within QM/MM.
Definition: qmmm_image_charge.F:14
qmmm_image_charge::calculate_image_pot
subroutine, public calculate_image_pot(v_hartree_rspace, rho_hartree_gspace, energy, qmmm_env, qs_env)
determines coefficients by solving image_matrix*coeff=-pot_const by Gaussian elimination or in an ite...
Definition: qmmm_image_charge.F:96
qmmm_image_charge::print_image_coefficients
subroutine, public print_image_coefficients(image_coeff, qs_env)
Print image coefficients.
Definition: qmmm_image_charge.F:1151
qmmm_image_charge::integrate_potential_devga_rspace
subroutine, public integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, qs_env)
calculates the image forces on the MM atoms
Definition: qmmm_image_charge.F:452
qmmm_image_charge::conditional_calc_image_matrix
subroutine, public conditional_calc_image_matrix(qs_env, qmmm_env)
calculate image matrix T depending on constraints on image atoms in case coefficients are estimated n...
Definition: qmmm_image_charge.F:574
qmmm_image_charge::add_image_pot_to_hartree_pot
subroutine, public add_image_pot_to_hartree_pot(v_hartree, v_metal, qs_env)
Add potential of metal (image charge pot) to Hartree Potential.
Definition: qmmm_image_charge.F:926
qmmm_types_low
Definition: qmmm_types_low.F:13
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_rho_metal
subroutine, public calculate_rho_metal(rho_metal, coeff, total_rho_metal, qs_env)
computes the image charge density on the grid (including coeffcients)
Definition: qs_collocate_density.F:1292
qs_collocate_density::calculate_rho_single_gaussian
subroutine, public calculate_rho_single_gaussian(rho_gb, qs_env, iatom_in)
collocate a single Gaussian on the grid
Definition: qs_collocate_density.F:1202
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition: qs_integrate_potential.F:22
realspace_grid_types
Definition: realspace_grid_types.F:18
realspace_grid_types::transfer_pw2rs
subroutine, public transfer_pw2rs(rs, pw)
...
Definition: realspace_grid_types.F:698
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333
virial_types
Definition: virial_types.F:13