(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
15MODULE 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
21 USE cp_ddapc_forces, ONLY: ewald_ddapc_force,&
22 reset_ch_pulay,&
23 restraint_functional_force,&
24 solvation_ddapc_force
25 USE cp_ddapc_util, ONLY: get_ddapc,&
26 modify_hartree_pot,&
27 restraint_functional_potential
28 USE cp_log_handling, ONLY: cp_get_default_logger,&
29 cp_logger_type
30 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
31 cp_print_key_unit_nr
32 USE dbcsr_api, ONLY: dbcsr_copy,&
33 dbcsr_p_type,&
34 dbcsr_set
35 USE input_constants, ONLY: do_spin_density
36 USE input_section_types, ONLY: section_vals_get_subs_vals,&
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
60 qs_ks_ddapc
61
62CONTAINS
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
554END MODULE cp_ddapc
pw_methods::pw_integral_ab
Definition pw_methods.F:221
pw_methods::pw_scale
Definition pw_methods.F:93
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:28
bibliography::blochl1995
integer, save, public blochl1995
Definition bibliography.F:43
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_ddapc_forces
Density Derived atomic point charges from a QM calculation (see J. Chem. Phys. Vol....
Definition cp_ddapc_forces.F:15
cp_ddapc_forces::solvation_ddapc_force
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...
Definition cp_ddapc_forces.F:458
cp_ddapc_forces::ewald_ddapc_force
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.
Definition cp_ddapc_forces.F:76
cp_ddapc_forces::restraint_functional_force
subroutine, public restraint_functional_force(qs_env, ddapc_restraint_control, dq, charges, n_gauss, particle_set)
computes derivatives for DDAPC restraint
Definition cp_ddapc_forces.F:390
cp_ddapc_forces::reset_ch_pulay
subroutine, public reset_ch_pulay(qs_env)
Evaluation of the pulay forces due to the fitted charge density.
Definition cp_ddapc_forces.F:298
cp_ddapc_util
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition cp_ddapc_util.F:15
cp_ddapc_util::get_ddapc
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.
Definition cp_ddapc_util.F:203
cp_ddapc_util::modify_hartree_pot
subroutine, public modify_hartree_pot(v_hartree_gspace, density_fit_section, particle_set, m, ami, radii, charges)
Modify the Hartree potential.
Definition cp_ddapc_util.F:612
cp_ddapc_util::restraint_functional_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
Definition cp_ddapc_util.F:535
cp_ddapc
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition cp_ddapc.F:15
cp_ddapc::cp_ddapc_apply_cd
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
cp_ddapc::qs_ks_ddapc
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
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_spin_density
integer, parameter, public do_spin_density
Definition input_constants.F:42
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::ddapc_restraint_type
Definition cp_control_types.F:258
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
particle_types::particle_type
Definition particle_types.F:35
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215