(git:0de0cc2)
cp_ddapc.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Density Derived atomic point charges from a QM calculation
10 !> (see Bloechl, J. Chem. Phys. Vol. 103 pp. 7422-7428)
11 !> \par History
12 !> 08.2005 created [tlaino]
13 !> \author Teodoro Laino
14 ! **************************************************************************************************
15 MODULE cp_ddapc
16  USE bibliography, ONLY: blochl1995,&
17  cite_reference
18  USE cell_types, ONLY: cell_type
19  USE cp_control_types, ONLY: ddapc_restraint_type,&
20  dft_control_type
25  USE cp_ddapc_util, ONLY: get_ddapc,&
29  cp_logger_type
32  USE dbcsr_api, ONLY: dbcsr_copy,&
33  dbcsr_p_type,&
34  dbcsr_set
37  section_vals_type
38  USE kinds, ONLY: dp
39  USE particle_types, ONLY: particle_type
40  USE pw_methods, ONLY: pw_integral_ab,&
41  pw_scale,&
42  pw_transfer,&
43  pw_zero
44  USE pw_pool_types, ONLY: pw_pool_type
45  USE pw_types, ONLY: pw_c1d_gs_type,&
46  pw_r3d_rs_type
47  USE qs_energy_types, ONLY: qs_energy_type
48  USE qs_environment_types, ONLY: get_qs_env,&
49  qs_environment_type
50  USE qs_integrate_potential, ONLY: integrate_v_rspace
51 #include "./base/base_uses.f90"
52 
53  IMPLICIT NONE
54  PRIVATE
55 
56  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
57  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc'
58 
59  PUBLIC :: cp_ddapc_apply_cd, & ! Apply Coupling/Decoupling to Periodic Images
61 
62 CONTAINS
63 
64 ! **************************************************************************************************
65 !> \brief Set of methods using DDAPC charges
66 !> \param qs_env ...
67 !> \param auxbas_pw_pool ...
68 !> \param rho_tot_gspace ...
69 !> \param v_hartree_gspace ...
70 !> \param v_spin_ddapc_rest_r ...
71 !> \param energy ...
72 !> \param calculate_forces ...
73 !> \param ks_matrix ...
74 !> \param just_energy ...
75 !> \par History
76 !> 08.2005 created [tlaino]
77 !> 08.2008 extended to restraint/constraint DDAPC charges [fschiff]
78 ! **************************************************************************************************
79  SUBROUTINE qs_ks_ddapc(qs_env, auxbas_pw_pool, rho_tot_gspace, v_hartree_gspace, &
80  v_spin_ddapc_rest_r, energy, calculate_forces, ks_matrix, just_energy)
81 
82  TYPE(qs_environment_type), POINTER :: qs_env
83  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
84  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_gspace, v_hartree_gspace
85  TYPE(pw_r3d_rs_type), POINTER :: v_spin_ddapc_rest_r
86  TYPE(qs_energy_type), POINTER :: energy
87  LOGICAL, INTENT(in) :: calculate_forces
88  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
89  LOGICAL, INTENT(in) :: just_energy
90 
91  CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_ks_ddapc'
92 
93  INTEGER :: ddapc_size, handle, i, my_id
94  LOGICAL :: ddapc_restraint_is_spin, &
95  et_coupling_calc, explicit_potential
96  TYPE(cp_logger_type), POINTER :: logger
97  TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control
98  TYPE(dft_control_type), POINTER :: dft_control
99  TYPE(pw_c1d_gs_type) :: v_spin_ddapc_rest_g
100  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
101 
102  NULLIFY (v_hartree_rspace, dft_control)
103 
104  CALL timeset(routinen, handle)
105  CALL cite_reference(blochl1995)
106  ! In case decouple periodic images and/or apply restraints to charges
107  logger => cp_get_default_logger()
108  ddapc_restraint_is_spin = .false.
109  et_coupling_calc = .false.
110  ddapc_size = 0
111 
112  ! no k-points
113  cpassert(SIZE(ks_matrix, 2) == 1)
114 
115  CALL get_qs_env(qs_env, &
116  v_hartree_rspace=v_hartree_rspace, &
117  dft_control=dft_control)
118 
119  IF (dft_control%qs_control%ddapc_restraint) THEN
120  ddapc_size = SIZE(dft_control%qs_control%ddapc_restraint_control)
121  IF (SIZE(energy%ddapc_restraint) .NE. ddapc_size) THEN
122  DEALLOCATE (energy%ddapc_restraint)
123  ALLOCATE (energy%ddapc_restraint(ddapc_size))
124  END IF
125 
126  DO i = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
127  my_id = dft_control%qs_control%ddapc_restraint_control(i)%density_type
128  IF (my_id == do_spin_density .OR. ddapc_restraint_is_spin) ddapc_restraint_is_spin = .true.
129  END DO
130  et_coupling_calc = dft_control%qs_control%et_coupling_calc
131  END IF
132 
133  explicit_potential = ddapc_restraint_is_spin .OR. et_coupling_calc
134  dft_control%qs_control%ddapc_explicit_potential = explicit_potential
135  dft_control%qs_control%ddapc_restraint_is_spin = ddapc_restraint_is_spin
136  IF (explicit_potential) THEN
137  CALL auxbas_pw_pool%create_pw(v_spin_ddapc_rest_g)
138  CALL pw_zero(v_spin_ddapc_rest_g)
139  NULLIFY (v_spin_ddapc_rest_r)
140  ALLOCATE (v_spin_ddapc_rest_r)
141  CALL auxbas_pw_pool%create_pw(v_spin_ddapc_rest_r)
142  END IF
143 
144  IF (calculate_forces) CALL reset_ch_pulay(qs_env)
145 
146  ! Decoupling/Recoupling
147  CALL cp_ddapc_apply_cd(qs_env, rho_tot_gspace, energy%hartree, v_hartree_gspace, &
148  calculate_forces, itype_of_density="FULL DENSITY")
149  IF (dft_control%qs_control%ddapc_restraint) THEN
150  ! Restraints/Constraints
151  DO i = 1, ddapc_size
152  NULLIFY (ddapc_restraint_control)
153  ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(i)
154 
155  CALL cp_ddapc_apply_rs(qs_env, energy%ddapc_restraint(i), v_hartree_gspace, &
156  v_spin_ddapc_rest_g, ddapc_restraint_control, calculate_forces)
157  END DO
158  END IF
159  CALL cp_ddapc_apply_rf(qs_env, rho_tot_gspace, energy%hartree, v_hartree_gspace, &
160  calculate_forces, itype_of_density="FULL DENSITY")
161 
162  ! CJM Copying the real-space Hartree potential to KS_ENV
163  IF ((.NOT. just_energy) .OR. et_coupling_calc) THEN
164  CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
165  CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
166  IF (explicit_potential) THEN
167  CALL pw_transfer(v_spin_ddapc_rest_g, v_spin_ddapc_rest_r)
168  CALL pw_scale(v_spin_ddapc_rest_r, v_spin_ddapc_rest_r%pw_grid%dvol)
169  IF (et_coupling_calc) THEN
170  IF (qs_env%et_coupling%keep_matrix) THEN
171  IF (qs_env%et_coupling%first_run) THEN
172  NULLIFY (qs_env%et_coupling%rest_mat(1)%matrix)
173  ALLOCATE (qs_env%et_coupling%rest_mat(1)%matrix)
174  CALL dbcsr_copy(qs_env%et_coupling%rest_mat(1)%matrix, ks_matrix(1, 1)%matrix, &
175  name="ET_RESTRAINT_MATRIX_B")
176  CALL dbcsr_set(qs_env%et_coupling%rest_mat(1)%matrix, 0.0_dp)
177  CALL integrate_v_rspace(v_spin_ddapc_rest_r, &
178  hmat=qs_env%et_coupling%rest_mat(1), &
179  qs_env=qs_env, calculate_forces=.false.)
180  qs_env%et_coupling%order_p = &
181  dft_control%qs_control%ddapc_restraint_control(1)%ddapc_order_p
182  qs_env%et_coupling%e1 = dft_control%qs_control%ddapc_restraint_control(1)%strength
183  qs_env%et_coupling%keep_matrix = .false.
184  ELSE
185  NULLIFY (qs_env%et_coupling%rest_mat(2)%matrix)
186  ALLOCATE (qs_env%et_coupling%rest_mat(2)%matrix)
187  CALL dbcsr_copy(qs_env%et_coupling%rest_mat(2)%matrix, ks_matrix(1, 1)%matrix, &
188  name="ET_RESTRAINT_MATRIX_B")
189  CALL dbcsr_set(qs_env%et_coupling%rest_mat(2)%matrix, 0.0_dp)
190  CALL integrate_v_rspace(v_spin_ddapc_rest_r, &
191  hmat=qs_env%et_coupling%rest_mat(2), &
192  qs_env=qs_env, calculate_forces=.false.)
193  END IF
194  END IF
195  END IF
196  END IF
197  END IF
198 
199  IF (explicit_potential) THEN
200  CALL auxbas_pw_pool%give_back_pw(v_spin_ddapc_rest_g)
201  END IF
202  CALL timestop(handle)
203 
204  END SUBROUTINE qs_ks_ddapc
205 
206 ! **************************************************************************************************
207 !> \brief Routine to couple/decouple periodic images with the Bloechl scheme
208 !>
209 !> The coupling/decoupling is obtaines evaluating terms E2 and E3 in
210 !> J. Chem. Phys. Vol. 103 pp. 7422-7428.. The E2 terms is just a
211 !> Ewald summation, and for performance reason I'm writing a specific
212 !> driver instead of using and setting-up the environment of the already
213 !> available routines
214 !> \param qs_env ...
215 !> \param rho_tot_gspace ...
216 !> \param energy ...
217 !> \param v_hartree_gspace ...
218 !> \param calculate_forces ...
219 !> \param Itype_of_density ...
220 !> \par History
221 !> 08.2005 created [tlaino]
222 !> \author Teodoro Laino
223 ! **************************************************************************************************
224  SUBROUTINE cp_ddapc_apply_cd(qs_env, rho_tot_gspace, energy, v_hartree_gspace, &
225  calculate_forces, Itype_of_density)
226  TYPE(qs_environment_type), POINTER :: qs_env
227  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_gspace
228  REAL(kind=dp), INTENT(INOUT) :: energy
229  TYPE(pw_c1d_gs_type), INTENT(IN) :: v_hartree_gspace
230  LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
231  CHARACTER(LEN=*) :: itype_of_density
232 
233  CHARACTER(len=*), PARAMETER :: routinen = 'cp_ddapc_apply_CD'
234 
235  INTEGER :: handle, iw
236  LOGICAL :: apply_decpl, need_f
237  REAL(kind=dp) :: e_decpl, e_recpl
238  REAL(kind=dp), DIMENSION(:), POINTER :: charges, radii
239  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dq
240  TYPE(cell_type), POINTER :: cell, super_cell
241  TYPE(cp_logger_type), POINTER :: logger
242  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
243  TYPE(section_vals_type), POINTER :: density_fit_section, force_env_section, &
244  multipole_section, poisson_section, &
245  qmmm_periodic_section
246 
247  CALL timeset(routinen, handle)
248  need_f = .false.
249  IF (PRESENT(calculate_forces)) need_f = calculate_forces
250  logger => cp_get_default_logger()
251  apply_decpl = qs_env%cp_ddapc_ewald%do_decoupling .OR. qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl
252  IF (apply_decpl) THEN
253  ! Initialize
254  NULLIFY (multipole_section, &
255  poisson_section, &
256  force_env_section, &
257  particle_set, &
258  qmmm_periodic_section, &
259  density_fit_section, &
260  charges, &
261  radii, &
262  dq, &
263  cell, &
264  super_cell)
265 
266  CALL get_qs_env(qs_env=qs_env, &
267  input=force_env_section, &
268  particle_set=particle_set, &
269  cell=cell, &
270  super_cell=super_cell)
271  cpassert(ASSOCIATED(qs_env%cp_ddapc_ewald))
272  poisson_section => section_vals_get_subs_vals(force_env_section, "DFT%POISSON")
273 
274  density_fit_section => section_vals_get_subs_vals(force_env_section, "DFT%DENSITY_FITTING")
275 
276  IF (qs_env%cp_ddapc_ewald%do_decoupling) THEN
277  multipole_section => section_vals_get_subs_vals(poisson_section, "MULTIPOLE")
278  END IF
279  IF (qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
280  qmmm_periodic_section => section_vals_get_subs_vals(force_env_section, "QMMM%PERIODIC")
281  multipole_section => section_vals_get_subs_vals(qmmm_periodic_section, "MULTIPOLE")
282  END IF
283  ! Start the real calculation
284  iw = cp_print_key_unit_nr(logger, multipole_section, "PROGRAM_RUN_INFO", &
285  extension=".fitChargeLog")
286  ! First we evaluate the charges at the corresponding SCF STEP
287  IF (need_f) THEN
288  CALL get_ddapc(qs_env, &
289  need_f, &
290  density_fit_section, &
291  qout1=charges, &
292  out_radii=radii, &
293  dq_out=dq, &
294  ext_rho_tot_g=rho_tot_gspace, &
295  itype_of_density=itype_of_density)
296  ELSE
297  CALL get_ddapc(qs_env, &
298  need_f, &
299  density_fit_section, &
300  qout1=charges, &
301  out_radii=radii, &
302  ext_rho_tot_g=rho_tot_gspace, &
303  itype_of_density=itype_of_density)
304  END IF
305  ! Evaluate the Ewald contribution to the decoupling/coupling E2 and E3
306  IF (iw > 0) THEN
307  e_decpl = 0.5_dp*dot_product(charges, matmul(qs_env%cp_ddapc_env%Md, charges))
308  WRITE (iw, fmt="(T3,A,T60,F20.10)") "Decoupling Energy: ", e_decpl
309  END IF
310  IF (qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .AND. (iw > 0)) THEN
311  e_recpl = 0.5_dp*dot_product(charges, matmul(qs_env%cp_ddapc_env%Mr, charges))
312  WRITE (iw, fmt="(T3,A,T60,F20.10)") "Recoupling Energy: ", e_recpl
313  END IF
314  CALL modify_hartree_pot(v_hartree_gspace, &
315  density_fit_section, &
316  particle_set, &
317  qs_env%cp_ddapc_env%Mt, &
318  qs_env%cp_ddapc_env%AmI, &
319  radii, &
320  charges)
321  ! Modify the Hartree potential due to the decoupling/recoupling
322  energy = 0.5_dp*pw_integral_ab(rho_tot_gspace, v_hartree_gspace)
323  IF (need_f) THEN
324  CALL ewald_ddapc_force(qs_env, qs_env%cp_ddapc_ewald%coeff_qm, &
325  .false., 1.0_dp, multipole_section, cell, particle_set, &
326  radii, dq, charges)
327  IF (qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
328  CALL ewald_ddapc_force(qs_env, qs_env%cp_ddapc_ewald%coeff_mm, &
329  .true., -1.0_dp, multipole_section, super_cell, particle_set, &
330  radii, dq, charges)
331  END IF
332  END IF
333  ! Clean the allocated arrays
334  DEALLOCATE (charges)
335  DEALLOCATE (radii)
336  IF (ASSOCIATED(dq)) THEN
337  DEALLOCATE (dq)
338  END IF
339  CALL cp_print_key_finished_output(iw, logger, multipole_section, &
340  "PROGRAM_RUN_INFO")
341  END IF
342  CALL timestop(handle)
343  END SUBROUTINE cp_ddapc_apply_cd
344 
345 ! **************************************************************************************************
346 !> \brief Routine to apply RESTRAINT/CONSTRAINTS to the density
347 !> with the Bloechl scheme
348 !> \param qs_env ...
349 !> \param energy_res ...
350 !> \param v_hartree_gspace ...
351 !> \param v_spin_ddapc_rest_g ...
352 !> \param ddapc_restraint_control ...
353 !> \param calculate_forces ...
354 !> \par History
355 !> 08.2005 created [tlaino]
356 !> \author Teodoro Laino
357 ! **************************************************************************************************
358  SUBROUTINE cp_ddapc_apply_rs(qs_env, energy_res, v_hartree_gspace, &
359  v_spin_ddapc_rest_g, ddapc_restraint_control, calculate_forces)
360  TYPE(qs_environment_type), POINTER :: qs_env
361  REAL(kind=dp), INTENT(INOUT), OPTIONAL :: energy_res
362  TYPE(pw_c1d_gs_type), INTENT(IN) :: v_hartree_gspace, v_spin_ddapc_rest_g
363  TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control
364  LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
365 
366  CHARACTER(len=*), PARAMETER :: routinen = 'cp_ddapc_apply_RS'
367 
368  INTEGER :: handle, iw, my_id
369  LOGICAL :: apply_restrain, need_f
370  REAL(kind=dp), DIMENSION(:), POINTER :: charges, radii
371  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dq
372  TYPE(cell_type), POINTER :: cell, super_cell
373  TYPE(cp_logger_type), POINTER :: logger
374  TYPE(dft_control_type), POINTER :: dft_control
375  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
376  TYPE(section_vals_type), POINTER :: density_fit_section, force_env_section, &
377  restraint_section
378 
379  CALL timeset(routinen, handle)
380  NULLIFY (dft_control, restraint_section, force_env_section, particle_set, &
381  charges, radii, dq, cell, density_fit_section, super_cell)
382  need_f = .false.
383 
384  CALL get_qs_env(qs_env=qs_env, &
385  input=force_env_section, &
386  particle_set=particle_set, &
387  cell=cell, &
388  super_cell=super_cell, &
389  dft_control=dft_control)
390 
391  IF (PRESENT(calculate_forces)) need_f = calculate_forces
392  apply_restrain = dft_control%qs_control%ddapc_restraint
393  logger => cp_get_default_logger()
394  IF (apply_restrain) THEN
395  ! Initialize
396  density_fit_section => section_vals_get_subs_vals(force_env_section, "DFT%DENSITY_FITTING")
397  restraint_section => section_vals_get_subs_vals(force_env_section, "DFT%QS%DDAPC_RESTRAINT")
398  iw = cp_print_key_unit_nr(logger, restraint_section, "PROGRAM_RUN_INFO", &
399  extension=".fitChargeLog")
400  ! First we evaluate the charges at the corresponding SCF STEP
401  my_id = ddapc_restraint_control%density_type
402  IF (need_f) THEN
403  CALL get_ddapc(qs_env, &
404  need_f, &
405  density_fit_section, &
406  density_type=my_id, &
407  qout1=charges, &
408  out_radii=radii, &
409  dq_out=dq, iwc=iw)
410  ELSE
411  CALL get_ddapc(qs_env, &
412  need_f, &
413  density_fit_section, &
414  density_type=my_id, &
415  qout1=charges, &
416  out_radii=radii, iwc=iw)
417  END IF
418 
419  ! Modify the Hartree potential due to the restrain or the v_spin_ddapc_rest_g
420  IF ((my_id == do_spin_density) .OR. dft_control%qs_control%et_coupling_calc) THEN
421  CALL restraint_functional_potential(v_spin_ddapc_rest_g, density_fit_section, &
422  particle_set, qs_env%cp_ddapc_env%AmI, radii, charges, &
423  ddapc_restraint_control, energy_res)
424  ELSE
425  CALL restraint_functional_potential(v_hartree_gspace, density_fit_section, &
426  particle_set, qs_env%cp_ddapc_env%AmI, radii, charges, &
427  ddapc_restraint_control, energy_res)
428  END IF
429 
430  IF (need_f) THEN
431  CALL restraint_functional_force(qs_env, &
432  ddapc_restraint_control, &
433  dq, &
434  charges, &
435  SIZE(radii), &
436  particle_set)
437  END IF
438  ! Clean the allocated arrays
439  DEALLOCATE (charges)
440  DEALLOCATE (radii)
441  IF (ASSOCIATED(dq)) THEN
442  DEALLOCATE (dq)
443  END IF
444  CALL cp_print_key_finished_output(iw, logger, restraint_section, &
445  "PROGRAM_RUN_INFO")
446  END IF
447  CALL timestop(handle)
448  END SUBROUTINE cp_ddapc_apply_rs
449 
450 ! **************************************************************************************************
451 !> \brief Routine to apply a reaction field during SCF (SCRF) with the Bloechl scheme
452 !> \param qs_env ...
453 !> \param rho_tot_gspace ...
454 !> \param energy ...
455 !> \param v_hartree_gspace ...
456 !> \param calculate_forces ...
457 !> \param Itype_of_density ...
458 !> \par History
459 !> 08.2005 created [tlaino]
460 !> \author Teodoro Laino
461 ! **************************************************************************************************
462  SUBROUTINE cp_ddapc_apply_rf(qs_env, rho_tot_gspace, energy, &
463  v_hartree_gspace, calculate_forces, Itype_of_density)
464  TYPE(qs_environment_type), POINTER :: qs_env
465  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_gspace
466  REAL(kind=dp), INTENT(INOUT) :: energy
467  TYPE(pw_c1d_gs_type), INTENT(IN) :: v_hartree_gspace
468  LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
469  CHARACTER(LEN=*) :: itype_of_density
470 
471  CHARACTER(len=*), PARAMETER :: routinen = 'cp_ddapc_apply_RF'
472 
473  INTEGER :: handle, iw
474  LOGICAL :: apply_solvation, need_f
475  REAL(kind=dp) :: e_recpl
476  REAL(kind=dp), DIMENSION(:), POINTER :: charges, radii
477  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dq
478  TYPE(cell_type), POINTER :: cell, super_cell
479  TYPE(cp_logger_type), POINTER :: logger
480  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
481  TYPE(section_vals_type), POINTER :: density_fit_section, force_env_section, &
482  solvation_section
483 
484  CALL timeset(routinen, handle)
485  need_f = .false.
486  IF (PRESENT(calculate_forces)) need_f = calculate_forces
487  logger => cp_get_default_logger()
488  apply_solvation = qs_env%cp_ddapc_ewald%do_solvation
489  IF (apply_solvation) THEN
490  ! Initialize
491  NULLIFY (force_env_section, particle_set, charges, &
492  radii, dq, cell, super_cell)
493 
494  CALL get_qs_env(qs_env=qs_env, &
495  input=force_env_section, &
496  particle_set=particle_set, &
497  cell=cell, &
498  super_cell=super_cell)
499 
500  solvation_section => section_vals_get_subs_vals(force_env_section, "DFT%SCRF")
501  ! Start the real calculation
502  iw = cp_print_key_unit_nr(logger, solvation_section, "PROGRAM_RUN_INFO", &
503  extension=".fitChargeLog")
504  density_fit_section => section_vals_get_subs_vals(force_env_section, "DFT%DENSITY_FITTING")
505  ! First we evaluate the charges at the corresponding SCF STEP
506  IF (need_f) THEN
507  CALL get_ddapc(qs_env, &
508  need_f, &
509  density_fit_section, &
510  qout1=charges, &
511  out_radii=radii, &
512  dq_out=dq, &
513  ext_rho_tot_g=rho_tot_gspace, &
514  itype_of_density=itype_of_density)
515  ELSE
516  CALL get_ddapc(qs_env, &
517  need_f, &
518  density_fit_section, &
519  qout1=charges, &
520  out_radii=radii, &
521  ext_rho_tot_g=rho_tot_gspace, &
522  itype_of_density=itype_of_density)
523  END IF
524  ! Evaluate the Ewald contribution to the decoupling/coupling E2 and E3
525  IF (iw > 0) THEN
526  e_recpl = 0.5_dp*dot_product(charges, matmul(qs_env%cp_ddapc_env%Ms, charges))
527  WRITE (iw, fmt="(T3,A,T60,F20.10)") "Solvation Energy: ", e_recpl
528  END IF
529  CALL modify_hartree_pot(v_hartree_gspace, &
530  density_fit_section, &
531  particle_set, &
532  qs_env%cp_ddapc_env%Ms, &
533  qs_env%cp_ddapc_env%AmI, &
534  radii, &
535  charges)
536  ! Modify the Hartree potential due to the reaction field
537  energy = 0.5_dp*pw_integral_ab(rho_tot_gspace, v_hartree_gspace)
538  IF (need_f) THEN
539  CALL solvation_ddapc_force(qs_env, solvation_section, particle_set, &
540  radii, dq, charges)
541  END IF
542  ! Clean the allocated arrays
543  DEALLOCATE (charges)
544  DEALLOCATE (radii)
545  IF (ASSOCIATED(dq)) THEN
546  DEALLOCATE (dq)
547  END IF
548  CALL cp_print_key_finished_output(iw, logger, solvation_section, &
549  "PROGRAM_RUN_INFO")
550  END IF
551  CALL timestop(handle)
552  END SUBROUTINE cp_ddapc_apply_rf
553 
554 END MODULE cp_ddapc
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public blochl1995
Definition: bibliography.F:43
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Density Derived atomic point charges from a QM calculation (see J. Chem. Phys. Vol....
subroutine, public solvation_ddapc_force(qs_env, solvation_section, particle_set, radii, dq, charges)
Evaluates the electrostatic potential due to a simple solvation model Spherical cavity in a dieletric...
recursive subroutine, public ewald_ddapc_force(qs_env, coeff, apply_qmmm_periodic, factor, multipole_section, cell, particle_set, radii, dq, charges)
Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling of periodic images.
subroutine, public restraint_functional_force(qs_env, ddapc_restraint_control, dq, charges, n_gauss, particle_set)
computes derivatives for DDAPC restraint
subroutine, public reset_ch_pulay(qs_env)
Evaluation of the pulay forces due to the fitted charge density.
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition: cp_ddapc_util.F:15
recursive subroutine, public get_ddapc(qs_env, calc_force, density_fit_section, density_type, qout1, qout2, out_radii, dq_out, ext_rho_tot_g, Itype_of_density, iwc)
Computes the Density Derived Atomic Point Charges.
subroutine, public modify_hartree_pot(v_hartree_gspace, density_fit_section, particle_set, M, AmI, radii, charges)
Modify the Hartree potential.
subroutine, public restraint_functional_potential(v_hartree_gspace, density_fit_section, particle_set, AmI, radii, charges, ddapc_restraint_control, energy_res)
modify hartree potential to handle restraints in DDAPC scheme
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition: cp_ddapc.F:15
subroutine, public cp_ddapc_apply_cd(qs_env, rho_tot_gspace, energy, v_hartree_gspace, calculate_forces, Itype_of_density)
Routine to couple/decouple periodic images with the Bloechl scheme.
Definition: cp_ddapc.F:226
subroutine, public qs_ks_ddapc(qs_env, auxbas_pw_pool, rho_tot_gspace, v_hartree_gspace, v_spin_ddapc_rest_r, energy, calculate_forces, ks_matrix, just_energy)
Set of methods using DDAPC charges.
Definition: cp_ddapc.F:81
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_spin_density
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Define the data structure for the particle information.
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Integrate single or product functions over a potential on a RS grid.