Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines 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 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 : !***
84 : CONTAINS
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 60 : 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 60 : CALL timeset(routineN, handle)
109 :
110 60 : IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
111 : !calculate preconditioner matrix for CG if necessary
112 12 : IF (qs_env%calc_image_preconditioner) THEN
113 2 : IF (qmmm_env%image_charge_pot%image_restart) THEN
114 : CALL restart_image_matrix(image_matrix=qs_env%image_matrix, &
115 0 : qs_env=qs_env, qmmm_env=qmmm_env)
116 : ELSE
117 : CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
118 2 : 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 12 : 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 48 : qs_env=qs_env)
129 : END IF
130 :
131 : ! calculate the image/metal potential with the optimized coefficients
132 60 : 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 60 : energy=energy, qs_env=qs_env)
137 :
138 60 : CALL timestop(handle)
139 :
140 60 : 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 48 : 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 48 : 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 48 : V0 = -qmmm_env%image_charge_pot%V0
171 48 : eta = qmmm_env%image_charge_pot%eta
172 48 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
173 :
174 144 : ALLOCATE (pot_const(natom))
175 48 : IF (.NOT. ASSOCIATED(coeff)) THEN
176 16 : ALLOCATE (coeff(natom))
177 : END IF
178 144 : coeff = 0.0_dp
179 :
180 : CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
181 48 : pot_const)
182 : !add integral V0*ga(r)
183 144 : 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 48 : pot_const, natom, info)
189 48 : CPASSERT(info == 0)
190 :
191 240 : coeff = pot_const
192 :
193 48 : DEALLOCATE (pot_const)
194 :
195 48 : CALL timestop(handle)
196 :
197 48 : 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 12 : 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 : REAL(KIND=dp), PARAMETER :: eps_rsnew = 1.0E-16_dp
217 :
218 : INTEGER :: handle, iter_steps, natom, output_unit
219 : REAL(KIND=dp) :: alpha, eta, rsnew, rsold, V0
220 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: Ad, d, pot_const, r, vmetal_const, z
221 : TYPE(cp_logger_type), POINTER :: logger
222 : TYPE(pw_r3d_rs_type) :: auxpot_Ad_rspace, v_metal_rspace_guess
223 : TYPE(section_vals_type), POINTER :: input
224 :
225 12 : CALL timeset(routineN, handle)
226 :
227 12 : NULLIFY (pot_const, vmetal_const, logger, input)
228 12 : logger => cp_get_default_logger()
229 :
230 : !minus sign V0: account for the fact that v_hartree has the opposite sign
231 12 : V0 = -qmmm_env%image_charge_pot%V0
232 12 : eta = qmmm_env%image_charge_pot%eta
233 12 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
234 :
235 36 : ALLOCATE (pot_const(natom))
236 24 : ALLOCATE (vmetal_const(natom))
237 24 : ALLOCATE (r(natom))
238 24 : ALLOCATE (d(natom))
239 24 : ALLOCATE (z(natom))
240 24 : ALLOCATE (Ad(natom))
241 12 : IF (.NOT. ASSOCIATED(coeff)) THEN
242 4 : ALLOCATE (coeff(natom))
243 : END IF
244 :
245 : CALL integrate_potential_ga_rspace(v_hartree_rspace, qmmm_env, qs_env, &
246 12 : pot_const)
247 :
248 : !add integral V0*ga(r)
249 36 : pot_const(:) = -pot_const(:) + V0*SQRT((pi/eta)**3)
250 :
251 : !initial guess for coeff
252 36 : coeff = 1.0_dp
253 36 : d = 0.0_dp
254 36 : z = 0.0_dp
255 36 : r = 0.0_dp
256 12 : rsold = 0.0_dp
257 12 : rsnew = 0.0_dp
258 12 : iter_steps = 0
259 :
260 : !calculate first guess of image/metal potential
261 : CALL calculate_potential_metal(v_metal_rspace=v_metal_rspace_guess, &
262 12 : coeff=coeff, qs_env=qs_env)
263 : CALL integrate_potential_ga_rspace(potential=v_metal_rspace_guess, &
264 12 : qmmm_env=qmmm_env, qs_env=qs_env, int_res=vmetal_const)
265 :
266 : ! modify coefficients iteratively
267 60 : r = pot_const - vmetal_const
268 276 : z = MATMUL(qs_env%image_matrix, r)
269 60 : d = z
270 36 : rsold = DOT_PRODUCT(r, z)
271 :
272 6 : DO
273 : !calculate A*d
274 54 : Ad = 0.0_dp
275 : CALL calculate_potential_metal(v_metal_rspace= &
276 18 : auxpot_Ad_rspace, coeff=d, qs_env=qs_env)
277 : CALL integrate_potential_ga_rspace(potential= &
278 : auxpot_Ad_rspace, qmmm_env=qmmm_env, &
279 18 : qs_env=qs_env, int_res=Ad)
280 :
281 54 : alpha = rsold/DOT_PRODUCT(d, Ad)
282 90 : coeff = coeff + alpha*d
283 :
284 90 : r = r - alpha*Ad
285 414 : z = MATMUL(qs_env%image_matrix, r)
286 54 : rsnew = DOT_PRODUCT(r, z)
287 18 : iter_steps = iter_steps + 1
288 : ! SQRT(rsnew) < 1.0E-08
289 18 : IF (rsnew < eps_rsnew) THEN
290 12 : CALL auxpot_Ad_rspace%release()
291 : EXIT
292 : END IF
293 30 : d = z + rsnew/rsold*d
294 6 : rsold = rsnew
295 24 : CALL auxpot_Ad_rspace%release()
296 : END DO
297 :
298 : ! print iteration info
299 : CALL get_qs_env(qs_env=qs_env, &
300 12 : input=input)
301 : output_unit = cp_print_key_unit_nr(logger, input, &
302 : "QMMM%PRINT%PROGRAM_RUN_INFO", &
303 12 : extension=".qmmmLog")
304 12 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A,T74,I7)") &
305 6 : "Number of iteration steps for determination of image coefficients:", iter_steps
306 : CALL cp_print_key_finished_output(output_unit, logger, input, &
307 12 : "QMMM%PRINT%PROGRAM_RUN_INFO")
308 :
309 12 : IF (iter_steps < 25) THEN
310 12 : qs_env%calc_image_preconditioner = .FALSE.
311 : ELSE
312 0 : qs_env%calc_image_preconditioner = .TRUE.
313 : END IF
314 :
315 12 : CALL v_metal_rspace_guess%release()
316 12 : DEALLOCATE (pot_const)
317 12 : DEALLOCATE (vmetal_const)
318 12 : DEALLOCATE (r)
319 12 : DEALLOCATE (d, z)
320 12 : DEALLOCATE (Ad)
321 :
322 12 : CALL timestop(handle)
323 :
324 24 : END SUBROUTINE calc_image_coeff_iterative
325 :
326 : ! ****************************************************************************
327 : !> \brief calculates the integral V(r)*ga(r)
328 : !> \param potential any potential
329 : !> \param qmmm_env qmmm environment
330 : !> \param qs_env qs environment
331 : !> \param int_res result of the integration
332 : !> \param atom_num atom index, needed when calculating image_matrix
333 : !> \param atom_num_ref index of reference atom, needed when calculating
334 : !> image_matrix
335 : ! **************************************************************************************************
336 98 : SUBROUTINE integrate_potential_ga_rspace(potential, qmmm_env, qs_env, int_res, &
337 : atom_num, atom_num_ref)
338 :
339 : TYPE(pw_r3d_rs_type), INTENT(IN) :: potential
340 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
341 : TYPE(qs_environment_type), POINTER :: qs_env
342 : REAL(KIND=dp), DIMENSION(:), POINTER :: int_res
343 : INTEGER, INTENT(IN), OPTIONAL :: atom_num, atom_num_ref
344 :
345 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential_ga_rspace'
346 :
347 : INTEGER :: atom_a, atom_b, atom_ref, handle, iatom, &
348 : j, k, natom, npme
349 98 : INTEGER, DIMENSION(:), POINTER :: cores
350 : REAL(KIND=dp) :: eps_rho_rspace, radius
351 : REAL(KIND=dp), DIMENSION(3) :: ra
352 98 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab
353 : TYPE(cell_type), POINTER :: cell
354 : TYPE(dft_control_type), POINTER :: dft_control
355 : TYPE(mp_para_env_type), POINTER :: para_env
356 : TYPE(pw_env_type), POINTER :: pw_env
357 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
358 : TYPE(realspace_grid_type), POINTER :: rs_v
359 :
360 98 : CALL timeset(routineN, handle)
361 :
362 98 : NULLIFY (cores, hab, cell, auxbas_rs_desc, pw_env, para_env, &
363 98 : dft_control, rs_v)
364 98 : ALLOCATE (hab(1, 1))
365 :
366 98 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
367 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
368 98 : auxbas_rs_grid=rs_v)
369 98 : CALL transfer_pw2rs(rs_v, potential)
370 :
371 : CALL get_qs_env(qs_env=qs_env, &
372 : cell=cell, &
373 : dft_control=dft_control, &
374 98 : para_env=para_env, pw_env=pw_env)
375 :
376 98 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
377 :
378 98 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
379 98 : k = 1
380 98 : IF (PRESENT(atom_num)) k = atom_num
381 :
382 98 : CALL reallocate(cores, 1, natom - k + 1)
383 294 : int_res = 0.0_dp
384 98 : npme = 0
385 290 : cores = 0
386 :
387 290 : DO iatom = k, natom
388 290 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
389 : ! replicated realspace grid, split the atoms up between procs
390 192 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
391 96 : npme = npme + 1
392 96 : cores(npme) = iatom
393 : END IF
394 : ELSE
395 0 : npme = npme + 1
396 0 : cores(npme) = iatom
397 : END IF
398 : END DO
399 :
400 194 : DO j = 1, npme
401 :
402 96 : iatom = cores(j)
403 96 : atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
404 :
405 96 : IF (PRESENT(atom_num) .AND. PRESENT(atom_num_ref)) THEN
406 : ! shift the function since potential only calculate for ref atom
407 6 : atom_b = qmmm_env%image_charge_pot%image_mm_list(k)
408 6 : atom_ref = qmmm_env%image_charge_pot%image_mm_list(atom_num_ref)
409 : ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell) &
410 : - pbc(qmmm_env%image_charge_pot%particles_all(atom_b)%r, cell) &
411 24 : + pbc(qmmm_env%image_charge_pot%particles_all(atom_ref)%r, cell)
412 :
413 : ELSE
414 90 : ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
415 : END IF
416 :
417 96 : hab(1, 1) = 0.0_dp
418 :
419 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
420 : ra=ra, rb=ra, rp=ra, &
421 : zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
422 96 : prefactor=1.0_dp, cutoff=1.0_dp)
423 :
424 : CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
425 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
426 : rs_v, hab, o1=0, o2=0, &
427 : radius=radius, calculate_forces=.FALSE., &
428 96 : use_subpatch=.TRUE., subpatch_pattern=0)
429 :
430 194 : int_res(iatom) = hab(1, 1)
431 :
432 : END DO
433 :
434 490 : CALL para_env%sum(int_res)
435 :
436 98 : DEALLOCATE (hab, cores)
437 :
438 98 : CALL timestop(handle)
439 :
440 98 : END SUBROUTINE integrate_potential_ga_rspace
441 :
442 : ! **************************************************************************************************
443 : !> \brief calculates the image forces on the MM atoms
444 : !> \param potential any potential, in this case: Hartree potential
445 : !> \param coeff expansion coefficients of the image charge density, i.e.
446 : !> rho_metal=sum_a c_a*g_a
447 : !> \param forces structure storing the force contribution of the image charges
448 : !> for the metal (MM) atoms
449 : !> \param qmmm_env qmmm environment
450 : !> \param qs_env qs environment
451 : ! **************************************************************************************************
452 20 : SUBROUTINE integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, &
453 : qs_env)
454 :
455 : TYPE(pw_r3d_rs_type), INTENT(IN) :: potential
456 : REAL(KIND=dp), DIMENSION(:), POINTER :: coeff
457 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
458 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
459 : TYPE(qs_environment_type), POINTER :: qs_env
460 :
461 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_potential_devga_rspace'
462 :
463 : INTEGER :: atom_a, handle, iatom, j, natom, npme
464 20 : INTEGER, DIMENSION(:), POINTER :: cores
465 : LOGICAL :: use_virial
466 : REAL(KIND=dp) :: eps_rho_rspace, radius
467 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
468 20 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
469 : TYPE(cell_type), POINTER :: cell
470 : TYPE(dft_control_type), POINTER :: dft_control
471 : TYPE(mp_para_env_type), POINTER :: para_env
472 : TYPE(pw_env_type), POINTER :: pw_env
473 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
474 : TYPE(realspace_grid_type), POINTER :: rs_v
475 : TYPE(virial_type), POINTER :: virial
476 :
477 20 : CALL timeset(routineN, handle)
478 :
479 20 : NULLIFY (cores, hab, pab, cell, auxbas_rs_desc, pw_env, para_env, &
480 20 : dft_control, rs_v, virial)
481 20 : use_virial = .FALSE.
482 :
483 20 : ALLOCATE (hab(1, 1))
484 20 : ALLOCATE (pab(1, 1))
485 :
486 20 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
487 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
488 20 : auxbas_rs_grid=rs_v)
489 20 : CALL transfer_pw2rs(rs_v, potential)
490 :
491 : CALL get_qs_env(qs_env=qs_env, &
492 : cell=cell, &
493 : dft_control=dft_control, &
494 : para_env=para_env, pw_env=pw_env, &
495 20 : virial=virial)
496 :
497 20 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
498 :
499 : IF (use_virial) THEN
500 0 : CPABORT("Virial not implemented for image charge method")
501 : END IF
502 :
503 20 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
504 :
505 20 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
506 :
507 20 : IF (.NOT. ASSOCIATED(forces)) THEN
508 30 : ALLOCATE (forces(3, natom))
509 : END IF
510 :
511 180 : forces(:, :) = 0.0_dp
512 :
513 20 : CALL reallocate(cores, 1, natom)
514 20 : npme = 0
515 60 : cores = 0
516 :
517 60 : DO iatom = 1, natom
518 60 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
519 : ! replicated realspace grid, split the atoms up between procs
520 40 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
521 20 : npme = npme + 1
522 20 : cores(npme) = iatom
523 : END IF
524 : ELSE
525 0 : npme = npme + 1
526 0 : cores(npme) = iatom
527 : END IF
528 : END DO
529 :
530 40 : DO j = 1, npme
531 :
532 20 : iatom = cores(j)
533 20 : atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
534 20 : ra(:) = pbc(qmmm_env%image_charge_pot%particles_all(atom_a)%r, cell)
535 20 : hab(1, 1) = 0.0_dp
536 20 : pab(1, 1) = 1.0_dp
537 20 : force_a(:) = 0.0_dp
538 20 : force_b(:) = 0.0_dp
539 :
540 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
541 : ra=ra, rb=ra, rp=ra, &
542 : zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, &
543 : pab=pab, o1=0, o2=0, & ! without map_consistent
544 20 : prefactor=1.0_dp, cutoff=1.0_dp)
545 :
546 : CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, &
547 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
548 : rs_v, hab, pab, o1=0, o2=0, &
549 : radius=radius, calculate_forces=.TRUE., &
550 : force_a=force_a, force_b=force_b, use_subpatch=.TRUE., &
551 20 : subpatch_pattern=0)
552 :
553 80 : force_a(:) = coeff(iatom)*force_a(:)
554 100 : forces(:, iatom) = force_a(:)
555 :
556 : END DO
557 :
558 340 : CALL para_env%sum(forces)
559 :
560 20 : DEALLOCATE (hab, pab, cores)
561 :
562 : ! print info on gradients if wanted
563 20 : CALL print_gradients_image_atoms(forces, qs_env)
564 :
565 20 : CALL timestop(handle)
566 :
567 20 : END SUBROUTINE integrate_potential_devga_rspace
568 :
569 : !****************************************************************************
570 : !> \brief calculate image matrix T depending on constraints on image atoms
571 : !> in case coefficients are estimated not iteratively
572 : !> \param qs_env qs environment
573 : !> \param qmmm_env qmmm environment
574 : ! **************************************************************************************************
575 20 : SUBROUTINE conditional_calc_image_matrix(qs_env, qmmm_env)
576 :
577 : TYPE(qs_environment_type), POINTER :: qs_env
578 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
579 :
580 20 : IF (.NOT. qmmm_env%image_charge_pot%coeff_iterative) THEN
581 28 : SELECT CASE (qmmm_env%image_charge_pot%state_image_matrix)
582 : CASE (calc_always)
583 : CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
584 12 : ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
585 : CASE (calc_once)
586 : !if all image atoms are fully constrained, calculate image matrix
587 : !only for the first MD or GEO_OPT step
588 : CALL calculate_image_matrix(image_matrix=qs_env%image_matrix, &
589 2 : ipiv=qs_env%ipiv, qs_env=qs_env, qmmm_env=qmmm_env)
590 2 : qmmm_env%image_charge_pot%state_image_matrix = calc_once_done
591 2 : IF (qmmm_env%center_qm_subsys0) THEN
592 : CALL cp_warn(__LOCATION__, &
593 : "The image atoms are fully "// &
594 : "constrained and the image matrix is only calculated once. "// &
595 0 : "To be safe, set CENTER to NEVER ")
596 : END IF
597 : CASE (calc_once_done)
598 : ! do nothing image matrix is stored
599 : CASE DEFAULT
600 16 : CPABORT("No initialization for image charges available?")
601 : END SELECT
602 : END IF
603 :
604 20 : END SUBROUTINE conditional_calc_image_matrix
605 :
606 : !****************************************************************************
607 : !> \brief calculate image matrix T
608 : !> \param image_matrix matrix T
609 : !> \param ipiv pivoting prior to DGETRS (for Gaussian elimination)
610 : !> \param qs_env qs environment
611 : !> \param qmmm_env qmmm environment
612 : ! **************************************************************************************************
613 16 : SUBROUTINE calculate_image_matrix(image_matrix, ipiv, qs_env, qmmm_env)
614 :
615 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: image_matrix
616 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: ipiv
617 : TYPE(qs_environment_type), POINTER :: qs_env
618 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
619 :
620 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_image_matrix'
621 :
622 : INTEGER :: handle, natom, output_unit, stat
623 : TYPE(cp_logger_type), POINTER :: logger
624 : TYPE(section_vals_type), POINTER :: input
625 :
626 16 : CALL timeset(routineN, handle)
627 16 : NULLIFY (input, logger)
628 :
629 16 : logger => cp_get_default_logger()
630 :
631 16 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
632 :
633 16 : IF (.NOT. ASSOCIATED(image_matrix)) THEN
634 40 : ALLOCATE (image_matrix(natom, natom))
635 : END IF
636 16 : IF (PRESENT(ipiv)) THEN
637 14 : IF (.NOT. ASSOCIATED(ipiv)) THEN
638 24 : ALLOCATE (ipiv(natom))
639 : END IF
640 42 : ipiv = 0
641 : END IF
642 :
643 16 : CALL get_qs_env(qs_env, input=input)
644 : !print info
645 : output_unit = cp_print_key_unit_nr(logger, input, &
646 : "QMMM%PRINT%PROGRAM_RUN_INFO", &
647 16 : extension=".qmmmLog")
648 16 : IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
649 2 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") &
650 1 : "Calculating image matrix"
651 : ELSE
652 14 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T2,A)") &
653 7 : "Calculating image matrix"
654 : END IF
655 : CALL cp_print_key_finished_output(output_unit, logger, input, &
656 16 : "QMMM%PRINT%PROGRAM_RUN_INFO")
657 :
658 : ! Calculate image matrix using either GPW or MME method
659 20 : SELECT CASE (qmmm_env%image_charge_pot%image_matrix_method)
660 : CASE (do_eri_gpw)
661 4 : CALL calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
662 : CASE (do_eri_mme)
663 12 : CALL calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
664 : CASE DEFAULT
665 16 : CPABORT("Unknown method for calculating image matrix")
666 : END SELECT
667 :
668 16 : IF (qmmm_env%image_charge_pot%coeff_iterative) THEN
669 : !inversion --> preconditioner matrix for CG
670 2 : CALL invmat_symm(qs_env%image_matrix)
671 2 : CALL write_image_matrix(qs_env%image_matrix, qs_env)
672 : ELSE
673 : !pivoting prior to DGETRS (Gaussian elimination)
674 14 : IF (PRESENT(ipiv)) THEN
675 14 : CALL dgetrf(natom, natom, image_matrix, natom, ipiv, stat)
676 14 : CPASSERT(stat == 0)
677 : END IF
678 : END IF
679 :
680 16 : CALL timestop(handle)
681 :
682 16 : END SUBROUTINE calculate_image_matrix
683 :
684 : ! **************************************************************************************************
685 : !> \brief calculate image matrix T using GPW method
686 : !> \param image_matrix matrix T
687 : !> \param qs_env qs environment
688 : !> \param qmmm_env qmmm environment
689 : ! **************************************************************************************************
690 4 : SUBROUTINE calculate_image_matrix_gpw(image_matrix, qs_env, qmmm_env)
691 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: image_matrix
692 : TYPE(qs_environment_type), POINTER :: qs_env
693 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
694 :
695 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_image_matrix_gpw'
696 :
697 : INTEGER :: handle, iatom, iatom_ref, natom
698 : REAL(KIND=dp), DIMENSION(:), POINTER :: int_res
699 : TYPE(mp_para_env_type), POINTER :: para_env
700 : TYPE(pw_c1d_gs_type) :: rho_gb, vb_gspace
701 : TYPE(pw_env_type), POINTER :: pw_env
702 : TYPE(pw_poisson_type), POINTER :: poisson_env
703 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
704 : TYPE(pw_r3d_rs_type) :: vb_rspace
705 :
706 4 : CALL timeset(routineN, handle)
707 4 : NULLIFY (pw_env, auxbas_pw_pool, poisson_env, para_env, int_res)
708 :
709 4 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
710 12 : ALLOCATE (int_res(natom))
711 :
712 28 : image_matrix = 0.0_dp
713 :
714 4 : CALL get_qs_env(qs_env, pw_env=pw_env, para_env=para_env)
715 :
716 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
717 4 : poisson_env=poisson_env)
718 4 : CALL auxbas_pw_pool%create_pw(rho_gb)
719 4 : CALL auxbas_pw_pool%create_pw(vb_gspace)
720 4 : CALL auxbas_pw_pool%create_pw(vb_rspace)
721 :
722 : ! calculate vb only once for one reference atom
723 4 : iatom_ref = 1 !
724 : !collocate gaussian of reference MM atom on grid
725 4 : CALL pw_zero(rho_gb)
726 4 : CALL calculate_rho_single_gaussian(rho_gb, qs_env, iatom_ref)
727 : !calculate potential vb like hartree potential
728 4 : CALL pw_zero(vb_gspace)
729 4 : CALL pw_poisson_solve(poisson_env, rho_gb, vhartree=vb_gspace)
730 4 : CALL pw_zero(vb_rspace)
731 4 : CALL pw_transfer(vb_gspace, vb_rspace)
732 4 : CALL pw_scale(vb_rspace, vb_rspace%pw_grid%dvol)
733 :
734 12 : DO iatom = 1, natom
735 : !calculate integral vb_rspace*ga
736 24 : int_res = 0.0_dp
737 : CALL integrate_potential_ga_rspace(vb_rspace, qs_env%qmmm_env_qm, &
738 : qs_env, int_res, atom_num=iatom, &
739 8 : atom_num_ref=iatom_ref)
740 40 : image_matrix(iatom, iatom:natom) = int_res(iatom:natom)
741 28 : image_matrix(iatom + 1:natom, iatom) = int_res(iatom + 1:natom)
742 : END DO
743 :
744 4 : CALL vb_gspace%release()
745 4 : CALL vb_rspace%release()
746 4 : CALL rho_gb%release()
747 :
748 4 : DEALLOCATE (int_res)
749 :
750 4 : CALL timestop(handle)
751 4 : END SUBROUTINE calculate_image_matrix_gpw
752 :
753 : ! **************************************************************************************************
754 : !> \brief calculate image matrix T using MME (MiniMax-Ewald) method
755 : !> \param image_matrix matrix T
756 : !> \param qs_env qs environment
757 : !> \param qmmm_env qmmm environment
758 : ! **************************************************************************************************
759 12 : SUBROUTINE calculate_image_matrix_mme(image_matrix, qs_env, qmmm_env)
760 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: image_matrix
761 : TYPE(qs_environment_type), POINTER :: qs_env
762 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
763 :
764 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_image_matrix_mme'
765 :
766 : INTEGER :: atom_a, handle, iatom, natom
767 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zeta
768 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ra
769 : TYPE(mp_para_env_type), POINTER :: para_env
770 :
771 12 : CALL timeset(routineN, handle)
772 12 : NULLIFY (para_env)
773 :
774 12 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
775 60 : ALLOCATE (zeta(natom), ra(3, natom))
776 :
777 36 : zeta(:) = qmmm_env%image_charge_pot%eta
778 :
779 36 : DO iatom = 1, natom
780 24 : atom_a = qmmm_env%image_charge_pot%image_mm_list(iatom)
781 108 : ra(:, iatom) = qmmm_env%image_charge_pot%particles_all(atom_a)%r(:)
782 : END DO
783 :
784 12 : CALL get_qs_env(qs_env, para_env=para_env)
785 :
786 : CALL integrate_s_mme(qmmm_env%image_charge_pot%eri_mme_param, &
787 12 : zeta, zeta, ra, ra, image_matrix, para_env)
788 :
789 12 : CALL timestop(handle)
790 24 : END SUBROUTINE calculate_image_matrix_mme
791 :
792 : ! **************************************************************************************************
793 : !> \brief high-level integration routine for 2c integrals over s-type functions.
794 : !> Parallelization over pairs of functions.
795 : !> \param param ...
796 : !> \param zeta ...
797 : !> \param zetb ...
798 : !> \param ra ...
799 : !> \param rb ...
800 : !> \param hab ...
801 : !> \param para_env ...
802 : ! **************************************************************************************************
803 12 : SUBROUTINE integrate_s_mme(param, zeta, zetb, ra, rb, hab, para_env)
804 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
805 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, zetb
806 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ra, rb
807 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: hab
808 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
809 :
810 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_s_mme'
811 :
812 : INTEGER :: G_count, handle, ipgf, ipgf_prod, jpgf, &
813 : npgf_prod, npgfa, npgfb, R_count
814 : INTEGER, DIMENSION(2) :: limits
815 : REAL(KIND=dp), DIMENSION(3) :: rab
816 :
817 12 : CALL timeset(routineN, handle)
818 12 : G_count = 0; R_count = 0
819 :
820 84 : hab(:, :) = 0.0_dp
821 :
822 12 : npgfa = SIZE(zeta)
823 12 : npgfb = SIZE(zetb)
824 12 : npgf_prod = npgfa*npgfb ! total number of integrals
825 :
826 12 : limits = get_limit(npgf_prod, para_env%num_pe, para_env%mepos)
827 :
828 36 : DO ipgf_prod = limits(1), limits(2)
829 24 : ipgf = (ipgf_prod - 1)/npgfb + 1
830 24 : jpgf = MOD(ipgf_prod - 1, npgfb) + 1
831 96 : rab(:) = ra(:, ipgf) - rb(:, jpgf)
832 : CALL eri_mme_2c_integrate(param%par, 0, 0, 0, 0, zeta(ipgf), &
833 36 : zetb(jpgf), rab, hab, ipgf - 1, jpgf - 1, G_count=G_count, R_count=R_count)
834 : END DO
835 :
836 12 : CALL cp_eri_mme_update_local_counts(param, para_env, G_count_2c=G_count, R_count_2c=R_count)
837 156 : CALL para_env%sum(hab)
838 12 : CALL timestop(handle)
839 :
840 12 : END SUBROUTINE integrate_s_mme
841 :
842 : ! **************************************************************************************************
843 : !> \brief calculates potential of the metal (image potential) given a set of
844 : !> coefficients coeff
845 : !> \param v_metal_rspace potential generated by rho_metal in real space
846 : !> \param coeff expansion coefficients of the image charge density, i.e.
847 : !> rho_metal=sum_a c_a*g_a
848 : !> \param rho_hartree_gspace Kohn Sham density in reciprocal space
849 : !> \param energy structure where energies are stored
850 : !> \param qs_env qs environment
851 : ! **************************************************************************************************
852 180 : SUBROUTINE calculate_potential_metal(v_metal_rspace, coeff, rho_hartree_gspace, energy, &
853 : qs_env)
854 :
855 : TYPE(pw_r3d_rs_type), INTENT(OUT) :: v_metal_rspace
856 : REAL(KIND=dp), DIMENSION(:), POINTER :: coeff
857 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_hartree_gspace
858 : TYPE(qs_energy_type), OPTIONAL, POINTER :: energy
859 : TYPE(qs_environment_type), POINTER :: qs_env
860 :
861 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_potential_metal'
862 :
863 : INTEGER :: handle
864 : REAL(KIND=dp) :: en_external, en_vmetal_rhohartree, &
865 : total_rho_metal
866 : TYPE(pw_c1d_gs_type) :: rho_metal, v_metal_gspace
867 : TYPE(pw_env_type), POINTER :: pw_env
868 : TYPE(pw_poisson_type), POINTER :: poisson_env
869 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
870 :
871 90 : CALL timeset(routineN, handle)
872 :
873 90 : NULLIFY (pw_env, auxbas_pw_pool, poisson_env)
874 : en_vmetal_rhohartree = 0.0_dp
875 : en_external = 0.0_dp
876 :
877 90 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
878 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
879 90 : poisson_env=poisson_env)
880 :
881 90 : CALL auxbas_pw_pool%create_pw(rho_metal)
882 :
883 90 : CALL auxbas_pw_pool%create_pw(v_metal_gspace)
884 :
885 90 : CALL auxbas_pw_pool%create_pw(v_metal_rspace)
886 :
887 90 : CALL pw_zero(rho_metal)
888 : CALL calculate_rho_metal(rho_metal, coeff, total_rho_metal=total_rho_metal, &
889 90 : qs_env=qs_env)
890 :
891 90 : CALL pw_zero(v_metal_gspace)
892 : CALL pw_poisson_solve(poisson_env, rho_metal, &
893 90 : vhartree=v_metal_gspace)
894 :
895 90 : IF (PRESENT(rho_hartree_gspace)) THEN
896 : en_vmetal_rhohartree = 0.5_dp*pw_integral_ab(v_metal_gspace, &
897 60 : rho_hartree_gspace)
898 60 : en_external = qs_env%qmmm_env_qm%image_charge_pot%V0*total_rho_metal
899 60 : energy%image_charge = en_vmetal_rhohartree - 0.5_dp*en_external
900 : CALL print_image_energy_terms(en_vmetal_rhohartree, en_external, &
901 60 : total_rho_metal, qs_env)
902 : END IF
903 :
904 90 : CALL pw_zero(v_metal_rspace)
905 90 : CALL pw_transfer(v_metal_gspace, v_metal_rspace)
906 90 : CALL pw_scale(v_metal_rspace, v_metal_rspace%pw_grid%dvol)
907 90 : CALL v_metal_gspace%release()
908 90 : CALL rho_metal%release()
909 :
910 90 : CALL timestop(handle)
911 :
912 90 : END SUBROUTINE calculate_potential_metal
913 :
914 : ! ****************************************************************************
915 : !> \brief Add potential of metal (image charge pot) to Hartree Potential
916 : !> \param v_hartree Hartree potential (in real space)
917 : !> \param v_metal potential generated by rho_metal (in real space)
918 : !> \param qs_env qs environment
919 : ! **************************************************************************************************
920 60 : SUBROUTINE add_image_pot_to_hartree_pot(v_hartree, v_metal, qs_env)
921 :
922 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_hartree
923 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_metal
924 : TYPE(qs_environment_type), POINTER :: qs_env
925 :
926 : CHARACTER(len=*), PARAMETER :: routineN = 'add_image_pot_to_hartree_pot'
927 :
928 : INTEGER :: handle, output_unit
929 : TYPE(cp_logger_type), POINTER :: logger
930 : TYPE(section_vals_type), POINTER :: input
931 :
932 60 : CALL timeset(routineN, handle)
933 :
934 60 : NULLIFY (input, logger)
935 60 : logger => cp_get_default_logger()
936 :
937 : !add image charge potential
938 60 : CALL pw_axpy(v_metal, v_hartree)
939 :
940 : ! print info
941 : CALL get_qs_env(qs_env=qs_env, &
942 60 : input=input)
943 : output_unit = cp_print_key_unit_nr(logger, input, &
944 : "QMMM%PRINT%PROGRAM_RUN_INFO", &
945 60 : extension=".qmmmLog")
946 60 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") &
947 30 : "Adding image charge potential to the Hartree potential."
948 : CALL cp_print_key_finished_output(output_unit, logger, input, &
949 60 : "QMMM%PRINT%PROGRAM_RUN_INFO")
950 :
951 60 : CALL timestop(handle)
952 :
953 60 : END SUBROUTINE add_image_pot_to_hartree_pot
954 :
955 : !****************************************************************************
956 : !> \brief writes image matrix T to file when used as preconditioner for
957 : !> calculating image coefficients iteratively
958 : !> \param image_matrix matrix T
959 : !> \param qs_env qs environment
960 : ! **************************************************************************************************
961 2 : SUBROUTINE write_image_matrix(image_matrix, qs_env)
962 :
963 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: image_matrix
964 : TYPE(qs_environment_type), POINTER :: qs_env
965 :
966 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_image_matrix'
967 :
968 : CHARACTER(LEN=default_path_length) :: filename
969 : INTEGER :: handle, rst_unit
970 : TYPE(cp_logger_type), POINTER :: logger
971 : TYPE(mp_para_env_type), POINTER :: para_env
972 : TYPE(section_vals_type), POINTER :: print_key, qmmm_section
973 :
974 2 : CALL timeset(routineN, handle)
975 :
976 2 : NULLIFY (qmmm_section, print_key, logger, para_env)
977 2 : logger => cp_get_default_logger()
978 : rst_unit = -1
979 :
980 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
981 2 : input=qmmm_section)
982 :
983 : print_key => section_vals_get_subs_vals(qmmm_section, &
984 2 : "QMMM%PRINT%IMAGE_CHARGE_RESTART")
985 :
986 2 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
987 : qmmm_section, "QMMM%PRINT%IMAGE_CHARGE_RESTART"), &
988 : cp_p_file)) THEN
989 :
990 : rst_unit = cp_print_key_unit_nr(logger, qmmm_section, &
991 : "QMMM%PRINT%IMAGE_CHARGE_RESTART", &
992 : extension=".Image", &
993 : file_status="REPLACE", &
994 : file_action="WRITE", &
995 2 : file_form="UNFORMATTED")
996 :
997 2 : IF (rst_unit > 0) filename = cp_print_key_generate_filename(logger, &
998 : print_key, extension=".IMAGE", &
999 1 : my_local=.FALSE.)
1000 :
1001 2 : IF (rst_unit > 0) THEN
1002 7 : WRITE (rst_unit) image_matrix
1003 : END IF
1004 :
1005 : CALL cp_print_key_finished_output(rst_unit, logger, qmmm_section, &
1006 2 : "QMMM%PRINT%IMAGE_CHARGE_RESTART")
1007 : END IF
1008 :
1009 2 : CALL timestop(handle)
1010 :
1011 2 : END SUBROUTINE write_image_matrix
1012 :
1013 : !****************************************************************************
1014 : !> \brief restarts image matrix T when used as preconditioner for calculating
1015 : !> image coefficients iteratively
1016 : !> \param image_matrix matrix T
1017 : !> \param qs_env qs environment
1018 : !> \param qmmm_env qmmm environment
1019 : ! **************************************************************************************************
1020 0 : SUBROUTINE restart_image_matrix(image_matrix, qs_env, qmmm_env)
1021 :
1022 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: image_matrix
1023 : TYPE(qs_environment_type), POINTER :: qs_env
1024 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
1025 :
1026 : CHARACTER(LEN=*), PARAMETER :: routineN = 'restart_image_matrix'
1027 :
1028 : CHARACTER(LEN=default_path_length) :: image_filename
1029 : INTEGER :: handle, natom, output_unit, rst_unit
1030 : LOGICAL :: exist
1031 : TYPE(cp_logger_type), POINTER :: logger
1032 : TYPE(mp_para_env_type), POINTER :: para_env
1033 : TYPE(section_vals_type), POINTER :: qmmm_section
1034 :
1035 0 : CALL timeset(routineN, handle)
1036 :
1037 0 : NULLIFY (qmmm_section, logger, para_env)
1038 0 : logger => cp_get_default_logger()
1039 0 : exist = .FALSE.
1040 0 : rst_unit = -1
1041 :
1042 0 : natom = SIZE(qmmm_env%image_charge_pot%image_mm_list)
1043 :
1044 0 : IF (.NOT. ASSOCIATED(image_matrix)) THEN
1045 0 : ALLOCATE (image_matrix(natom, natom))
1046 : END IF
1047 :
1048 0 : image_matrix = 0.0_dp
1049 :
1050 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, &
1051 0 : input=qmmm_section)
1052 :
1053 : CALL section_vals_val_get(qmmm_section, "QMMM%IMAGE_CHARGE%IMAGE_RESTART_FILE_NAME", &
1054 0 : c_val=image_filename)
1055 :
1056 0 : INQUIRE (FILE=image_filename, exist=exist)
1057 :
1058 0 : IF (exist) THEN
1059 0 : IF (para_env%is_source()) THEN
1060 : CALL open_file(file_name=image_filename, &
1061 : file_status="OLD", &
1062 : file_form="UNFORMATTED", &
1063 : file_action="READ", &
1064 0 : unit_number=rst_unit)
1065 :
1066 0 : READ (rst_unit) qs_env%image_matrix
1067 : END IF
1068 :
1069 0 : CALL para_env%bcast(qs_env%image_matrix)
1070 :
1071 0 : IF (para_env%is_source()) CALL close_file(unit_number=rst_unit)
1072 :
1073 : output_unit = cp_print_key_unit_nr(logger, qmmm_section, &
1074 : "QMMM%PRINT%PROGRAM_RUN_INFO", &
1075 0 : extension=".qmmmLog")
1076 0 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") &
1077 0 : "Restarted image matrix"
1078 : ELSE
1079 0 : CPABORT("Restart file for image matrix not found")
1080 : END IF
1081 :
1082 0 : qmmm_env%image_charge_pot%image_restart = .FALSE.
1083 :
1084 0 : CALL timestop(handle)
1085 :
1086 0 : END SUBROUTINE restart_image_matrix
1087 :
1088 : ! ****************************************************************************
1089 : !> \brief Print info on image gradients on image MM atoms
1090 : !> \param forces structure storing the force contribution of the image charges
1091 : !> for the metal (MM) atoms (actually these are only the gradients)
1092 : !> \param qs_env qs environment
1093 : ! **************************************************************************************************
1094 20 : SUBROUTINE print_gradients_image_atoms(forces, qs_env)
1095 :
1096 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
1097 : TYPE(qs_environment_type), POINTER :: qs_env
1098 :
1099 : INTEGER :: atom_a, iatom, natom, output_unit
1100 : REAL(KIND=dp), DIMENSION(3) :: sum_gradients
1101 : TYPE(cp_logger_type), POINTER :: logger
1102 : TYPE(section_vals_type), POINTER :: input
1103 :
1104 20 : NULLIFY (input, logger)
1105 20 : logger => cp_get_default_logger()
1106 :
1107 20 : sum_gradients = 0.0_dp
1108 20 : natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1109 :
1110 60 : DO iatom = 1, natom
1111 180 : sum_gradients(:) = sum_gradients(:) + forces(:, iatom)
1112 : END DO
1113 :
1114 20 : CALL get_qs_env(qs_env=qs_env, input=input)
1115 :
1116 : output_unit = cp_print_key_unit_nr(logger, input, &
1117 20 : "QMMM%PRINT%DERIVATIVES", extension=".Log")
1118 20 : IF (output_unit > 0) THEN
1119 : WRITE (unit=output_unit, fmt="(/1X,A)") &
1120 0 : "Image gradients [a.u.] on MM image charge atoms after QMMM calculation: "
1121 : WRITE (unit=output_unit, fmt="(T4,A4,T27,A1,T50,A1,T74,A1)") &
1122 0 : "Atom", "X", "Y", "Z"
1123 0 : DO iatom = 1, natom
1124 0 : atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1125 : WRITE (unit=output_unit, fmt="(T2,I6,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
1126 0 : atom_a, forces(:, iatom)
1127 : END DO
1128 :
1129 0 : WRITE (unit=output_unit, fmt="(T2,A)") REPEAT("-", 79)
1130 : WRITE (unit=output_unit, fmt="(T2,A,T22,ES12.5,T45,ES12.5,T69,ES12.5)") &
1131 0 : "sum gradients:", sum_gradients
1132 0 : WRITE (unit=output_unit, fmt="(/)")
1133 : END IF
1134 :
1135 : CALL cp_print_key_finished_output(output_unit, logger, input, &
1136 20 : "QMMM%PRINT%DERIVATIVES")
1137 :
1138 20 : END SUBROUTINE print_gradients_image_atoms
1139 :
1140 : ! ****************************************************************************
1141 : !> \brief Print image coefficients
1142 : !> \param image_coeff expansion coefficients of the image charge density
1143 : !> \param qs_env qs environment
1144 : ! **************************************************************************************************
1145 10 : SUBROUTINE print_image_coefficients(image_coeff, qs_env)
1146 :
1147 : REAL(KIND=dp), DIMENSION(:), POINTER :: image_coeff
1148 : TYPE(qs_environment_type), POINTER :: qs_env
1149 :
1150 : INTEGER :: atom_a, iatom, natom, output_unit
1151 : REAL(KIND=dp) :: normalize_factor, sum_coeff
1152 : TYPE(cp_logger_type), POINTER :: logger
1153 : TYPE(section_vals_type), POINTER :: input
1154 :
1155 10 : NULLIFY (input, logger)
1156 10 : logger => cp_get_default_logger()
1157 :
1158 10 : sum_coeff = 0.0_dp
1159 10 : natom = SIZE(qs_env%qmmm_env_qm%image_charge_pot%image_mm_list)
1160 10 : normalize_factor = SQRT((qs_env%qmmm_env_qm%image_charge_pot%eta/pi)**3)
1161 :
1162 30 : DO iatom = 1, natom
1163 30 : sum_coeff = sum_coeff + image_coeff(iatom)
1164 : END DO
1165 :
1166 10 : CALL get_qs_env(qs_env=qs_env, input=input)
1167 :
1168 : output_unit = cp_print_key_unit_nr(logger, input, &
1169 10 : "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=".Log")
1170 10 : IF (output_unit > 0) THEN
1171 2 : WRITE (unit=output_unit, fmt="(/)")
1172 : WRITE (unit=output_unit, fmt="(T2,A)") &
1173 2 : "Image charges [a.u.] after QMMM calculation: "
1174 2 : WRITE (unit=output_unit, fmt="(T4,A4,T67,A)") "Atom", "Image charge"
1175 2 : WRITE (unit=output_unit, fmt="(T4,A,T67,A)") REPEAT("-", 4), REPEAT("-", 12)
1176 :
1177 6 : DO iatom = 1, natom
1178 4 : atom_a = qs_env%qmmm_env_qm%image_charge_pot%image_mm_list(iatom)
1179 : !opposite sign for image_coeff; during the calculation they have
1180 : !the 'wrong' sign to ensure consistency with v_hartree which has
1181 : !the opposite sign
1182 : WRITE (unit=output_unit, fmt="(T2,I6,T65,ES16.9)") &
1183 6 : atom_a, -image_coeff(iatom)/normalize_factor
1184 : END DO
1185 :
1186 2 : WRITE (unit=output_unit, fmt="(T2,A)") REPEAT("-", 79)
1187 : WRITE (unit=output_unit, fmt="(T2,A,T65,ES16.9)") &
1188 2 : "sum image charges:", -sum_coeff/normalize_factor
1189 : END IF
1190 :
1191 : CALL cp_print_key_finished_output(output_unit, logger, input, &
1192 10 : "QMMM%PRINT%IMAGE_CHARGE_INFO")
1193 :
1194 10 : END SUBROUTINE print_image_coefficients
1195 :
1196 : ! ****************************************************************************
1197 : !> \brief Print detailed image charge energies
1198 : !> \param en_vmetal_rhohartree energy contribution of the image charges
1199 : !> without external potential, i.e. 0.5*integral(v_metal*rho_hartree)
1200 : !> \param en_external additional energy contribution of the image charges due
1201 : !> to an external potential, i.e. V0*total_rho_metal
1202 : !> \param total_rho_metal total induced image charge density
1203 : !> \param qs_env qs environment
1204 : ! **************************************************************************************************
1205 60 : SUBROUTINE print_image_energy_terms(en_vmetal_rhohartree, en_external, &
1206 : total_rho_metal, qs_env)
1207 :
1208 : REAL(KIND=dp), INTENT(IN) :: en_vmetal_rhohartree, en_external, &
1209 : total_rho_metal
1210 : TYPE(qs_environment_type), POINTER :: qs_env
1211 :
1212 : INTEGER :: output_unit
1213 : TYPE(cp_logger_type), POINTER :: logger
1214 : TYPE(section_vals_type), POINTER :: input
1215 :
1216 60 : NULLIFY (input, logger)
1217 60 : logger => cp_get_default_logger()
1218 :
1219 60 : CALL get_qs_env(qs_env=qs_env, input=input)
1220 :
1221 : output_unit = cp_print_key_unit_nr(logger, input, &
1222 60 : "QMMM%PRINT%IMAGE_CHARGE_INFO", extension=".Log")
1223 :
1224 60 : IF (output_unit > 0) THEN
1225 : WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1226 6 : "Total induced charge density [a.u.]:", total_rho_metal
1227 6 : WRITE (unit=output_unit, fmt="(T3,A)") "Image energy terms: "
1228 : WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1229 6 : "Coulomb energy of QM and image charge density [a.u.]:", en_vmetal_rhohartree
1230 : WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1231 6 : "External potential energy term [a.u.]:", -0.5_dp*en_external
1232 : WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
1233 6 : "Total image charge energy [a.u.]:", en_vmetal_rhohartree - 0.5_dp*en_external
1234 : END IF
1235 :
1236 : CALL cp_print_key_finished_output(output_unit, logger, input, &
1237 60 : "QMMM%PRINT%IMAGE_CHARGE_INFO")
1238 :
1239 60 : END SUBROUTINE print_image_energy_terms
1240 :
1241 30 : END MODULE qmmm_image_charge
|
