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