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