Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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
21 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
22 : dbcsr_p_type,&
23 : dbcsr_set
24 : USE cp_ddapc_forces, ONLY: ewald_ddapc_force,&
25 : reset_ch_pulay,&
26 : restraint_functional_force,&
27 : solvation_ddapc_force
28 : USE cp_ddapc_util, ONLY: get_ddapc,&
29 : modify_hartree_pot,&
30 : restraint_functional_potential
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
34 : cp_print_key_unit_nr
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 :
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 1378 : 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 1378 : NULLIFY (v_hartree_rspace, dft_control)
103 :
104 1378 : CALL timeset(routineN, handle)
105 1378 : CALL cite_reference(Blochl1995)
106 : ! In case decouple periodic images and/or apply restraints to charges
107 1378 : logger => cp_get_default_logger()
108 1378 : ddapc_restraint_is_spin = .FALSE.
109 1378 : et_coupling_calc = .FALSE.
110 1378 : ddapc_size = 0
111 :
112 : ! no k-points
113 1378 : CPASSERT(SIZE(ks_matrix, 2) == 1)
114 :
115 : CALL get_qs_env(qs_env, &
116 : v_hartree_rspace=v_hartree_rspace, &
117 1378 : dft_control=dft_control)
118 :
119 1378 : IF (dft_control%qs_control%ddapc_restraint) THEN
120 732 : ddapc_size = SIZE(dft_control%qs_control%ddapc_restraint_control)
121 732 : IF (SIZE(energy%ddapc_restraint) .NE. ddapc_size) THEN
122 2 : DEALLOCATE (energy%ddapc_restraint)
123 6 : ALLOCATE (energy%ddapc_restraint(ddapc_size))
124 : END IF
125 :
126 1568 : DO i = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
127 836 : my_id = dft_control%qs_control%ddapc_restraint_control(i)%density_type
128 1568 : IF (my_id == do_spin_density .OR. ddapc_restraint_is_spin) ddapc_restraint_is_spin = .TRUE.
129 : END DO
130 732 : et_coupling_calc = dft_control%qs_control%et_coupling_calc
131 : END IF
132 :
133 1378 : explicit_potential = ddapc_restraint_is_spin .OR. et_coupling_calc
134 1378 : dft_control%qs_control%ddapc_explicit_potential = explicit_potential
135 1378 : dft_control%qs_control%ddapc_restraint_is_spin = ddapc_restraint_is_spin
136 1378 : IF (explicit_potential) THEN
137 164 : CALL auxbas_pw_pool%create_pw(v_spin_ddapc_rest_g)
138 164 : CALL pw_zero(v_spin_ddapc_rest_g)
139 : NULLIFY (v_spin_ddapc_rest_r)
140 164 : ALLOCATE (v_spin_ddapc_rest_r)
141 164 : CALL auxbas_pw_pool%create_pw(v_spin_ddapc_rest_r)
142 : END IF
143 :
144 1378 : 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 1378 : calculate_forces, Itype_of_density="FULL DENSITY")
149 1378 : IF (dft_control%qs_control%ddapc_restraint) THEN
150 : ! Restraints/Constraints
151 1568 : DO i = 1, ddapc_size
152 : NULLIFY (ddapc_restraint_control)
153 836 : 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 1568 : 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 1378 : calculate_forces, Itype_of_density="FULL DENSITY")
161 :
162 : ! CJM Copying the real-space Hartree potential to KS_ENV
163 1378 : IF ((.NOT. just_energy) .OR. et_coupling_calc) THEN
164 988 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
165 988 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
166 988 : IF (explicit_potential) THEN
167 92 : CALL pw_transfer(v_spin_ddapc_rest_g, v_spin_ddapc_rest_r)
168 92 : CALL pw_scale(v_spin_ddapc_rest_r, v_spin_ddapc_rest_r%pw_grid%dvol)
169 92 : IF (et_coupling_calc) THEN
170 0 : IF (qs_env%et_coupling%keep_matrix) THEN
171 0 : IF (qs_env%et_coupling%first_run) THEN
172 0 : NULLIFY (qs_env%et_coupling%rest_mat(1)%matrix)
173 0 : 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 0 : name="ET_RESTRAINT_MATRIX_B")
176 0 : 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 0 : qs_env=qs_env, calculate_forces=.FALSE.)
180 : qs_env%et_coupling%order_p = &
181 0 : dft_control%qs_control%ddapc_restraint_control(1)%ddapc_order_p
182 0 : qs_env%et_coupling%e1 = dft_control%qs_control%ddapc_restraint_control(1)%strength
183 0 : qs_env%et_coupling%keep_matrix = .FALSE.
184 : ELSE
185 0 : NULLIFY (qs_env%et_coupling%rest_mat(2)%matrix)
186 0 : 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 0 : name="ET_RESTRAINT_MATRIX_B")
189 0 : 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 0 : qs_env=qs_env, calculate_forces=.FALSE.)
193 : END IF
194 : END IF
195 : END IF
196 : END IF
197 : END IF
198 :
199 482 : IF (explicit_potential) THEN
200 164 : CALL auxbas_pw_pool%give_back_pw(v_spin_ddapc_rest_g)
201 : END IF
202 1378 : CALL timestop(handle)
203 :
204 1378 : 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 1674 : 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 1674 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges, radii
239 1674 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dq
240 : TYPE(cell_type), POINTER :: cell, super_cell
241 : TYPE(cp_logger_type), POINTER :: logger
242 1674 : 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 1674 : CALL timeset(routineN, handle)
248 1674 : need_f = .FALSE.
249 1674 : IF (PRESENT(calculate_forces)) need_f = calculate_forces
250 1674 : logger => cp_get_default_logger()
251 1674 : 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 682 : poisson_section, &
256 682 : force_env_section, &
257 682 : particle_set, &
258 682 : qmmm_periodic_section, &
259 682 : density_fit_section, &
260 682 : charges, &
261 682 : radii, &
262 682 : dq, &
263 682 : cell, &
264 682 : 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 682 : super_cell=super_cell)
271 682 : CPASSERT(ASSOCIATED(qs_env%cp_ddapc_ewald))
272 682 : poisson_section => section_vals_get_subs_vals(force_env_section, "DFT%POISSON")
273 :
274 682 : density_fit_section => section_vals_get_subs_vals(force_env_section, "DFT%DENSITY_FITTING")
275 :
276 682 : IF (qs_env%cp_ddapc_ewald%do_decoupling) THEN
277 346 : multipole_section => section_vals_get_subs_vals(poisson_section, "MULTIPOLE")
278 : END IF
279 682 : IF (qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
280 336 : qmmm_periodic_section => section_vals_get_subs_vals(force_env_section, "QMMM%PERIODIC")
281 336 : 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 682 : extension=".fitChargeLog")
286 : ! First we evaluate the charges at the corresponding SCF STEP
287 682 : 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 90 : 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 592 : Itype_of_density=Itype_of_density)
304 : END IF
305 : ! Evaluate the Ewald contribution to the decoupling/coupling E2 and E3
306 682 : IF (iw > 0) THEN
307 237300 : e_decpl = 0.5_dp*DOT_PRODUCT(charges, MATMUL(qs_env%cp_ddapc_env%Md, charges))
308 357 : WRITE (iw, FMT="(T3,A,T60,F20.10)") "Decoupling Energy: ", e_decpl
309 : END IF
310 682 : IF (qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .AND. (iw > 0)) THEN
311 169450 : e_recpl = 0.5_dp*DOT_PRODUCT(charges, MATMUL(qs_env%cp_ddapc_env%Mr, charges))
312 184 : 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 682 : charges)
321 : ! Modify the Hartree potential due to the decoupling/recoupling
322 682 : energy = 0.5_dp*pw_integral_ab(rho_tot_gspace, v_hartree_gspace)
323 682 : 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 90 : radii, dq, charges)
327 90 : 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 46 : radii, dq, charges)
331 : END IF
332 : END IF
333 : ! Clean the allocated arrays
334 682 : DEALLOCATE (charges)
335 682 : DEALLOCATE (radii)
336 682 : IF (ASSOCIATED(dq)) THEN
337 90 : DEALLOCATE (dq)
338 : END IF
339 : CALL cp_print_key_finished_output(iw, logger, multipole_section, &
340 682 : "PROGRAM_RUN_INFO")
341 : END IF
342 1674 : CALL timestop(handle)
343 1674 : 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 836 : 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 836 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges, radii
371 836 : 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 836 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
376 : TYPE(section_vals_type), POINTER :: density_fit_section, force_env_section, &
377 : restraint_section
378 :
379 836 : CALL timeset(routineN, handle)
380 836 : NULLIFY (dft_control, restraint_section, force_env_section, particle_set, &
381 836 : charges, radii, dq, cell, density_fit_section, super_cell)
382 836 : 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 836 : dft_control=dft_control)
390 :
391 836 : IF (PRESENT(calculate_forces)) need_f = calculate_forces
392 836 : apply_restrain = dft_control%qs_control%ddapc_restraint
393 836 : logger => cp_get_default_logger()
394 836 : IF (apply_restrain) THEN
395 : ! Initialize
396 836 : density_fit_section => section_vals_get_subs_vals(force_env_section, "DFT%DENSITY_FITTING")
397 836 : 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 836 : extension=".fitChargeLog")
400 : ! First we evaluate the charges at the corresponding SCF STEP
401 836 : my_id = ddapc_restraint_control%density_type
402 836 : 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 36 : 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 800 : 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 836 : 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 164 : 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 672 : ddapc_restraint_control, energy_res)
428 : END IF
429 :
430 836 : IF (need_f) THEN
431 : CALL restraint_functional_force(qs_env, &
432 : ddapc_restraint_control, &
433 : dq, &
434 : charges, &
435 : SIZE(radii), &
436 36 : particle_set)
437 : END IF
438 : ! Clean the allocated arrays
439 836 : DEALLOCATE (charges)
440 836 : DEALLOCATE (radii)
441 836 : IF (ASSOCIATED(dq)) THEN
442 36 : DEALLOCATE (dq)
443 : END IF
444 : CALL cp_print_key_finished_output(iw, logger, restraint_section, &
445 836 : "PROGRAM_RUN_INFO")
446 : END IF
447 836 : CALL timestop(handle)
448 836 : 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 1378 : 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 1378 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges, radii
477 1378 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dq
478 : TYPE(cell_type), POINTER :: cell, super_cell
479 : TYPE(cp_logger_type), POINTER :: logger
480 1378 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
481 : TYPE(section_vals_type), POINTER :: density_fit_section, force_env_section, &
482 : solvation_section
483 :
484 1378 : CALL timeset(routineN, handle)
485 1378 : need_f = .FALSE.
486 1378 : IF (PRESENT(calculate_forces)) need_f = calculate_forces
487 1378 : logger => cp_get_default_logger()
488 1378 : apply_solvation = qs_env%cp_ddapc_ewald%do_solvation
489 1378 : IF (apply_solvation) THEN
490 : ! Initialize
491 144 : NULLIFY (force_env_section, particle_set, charges, &
492 144 : 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 144 : super_cell=super_cell)
499 :
500 144 : 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 144 : extension=".fitChargeLog")
504 144 : 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 144 : 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 14 : 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 130 : Itype_of_density=Itype_of_density)
523 : END IF
524 : ! Evaluate the Ewald contribution to the decoupling/coupling E2 and E3
525 144 : IF (iw > 0) THEN
526 16452 : e_recpl = 0.5_dp*DOT_PRODUCT(charges, MATMUL(qs_env%cp_ddapc_env%Ms, charges))
527 72 : 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 144 : charges)
536 : ! Modify the Hartree potential due to the reaction field
537 144 : energy = 0.5_dp*pw_integral_ab(rho_tot_gspace, v_hartree_gspace)
538 144 : IF (need_f) THEN
539 : CALL solvation_ddapc_force(qs_env, solvation_section, particle_set, &
540 14 : radii, dq, charges)
541 : END IF
542 : ! Clean the allocated arrays
543 144 : DEALLOCATE (charges)
544 144 : DEALLOCATE (radii)
545 144 : IF (ASSOCIATED(dq)) THEN
546 14 : DEALLOCATE (dq)
547 : END IF
548 : CALL cp_print_key_finished_output(iw, logger, solvation_section, &
549 144 : "PROGRAM_RUN_INFO")
550 : END IF
551 1378 : CALL timestop(handle)
552 1378 : END SUBROUTINE cp_ddapc_apply_RF
553 :
554 613 : END MODULE cp_ddapc
|
