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_util
16 :
17 : USE atomic_charges, ONLY: print_atomic_charges
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: evaluate_restraint_functional
22 : USE cp_ddapc_methods, ONLY: build_A_matrix,&
23 : build_b_vector,&
24 : build_der_A_matrix_rows,&
25 : build_der_b_vector,&
26 : cleanup_g_dot_rvec_sin_cos,&
27 : prep_g_dot_rvec_sin_cos
28 : USE cp_ddapc_types, ONLY: cp_ddapc_create,&
29 : cp_ddapc_type
30 : USE cp_log_handling, ONLY: cp_get_default_logger,&
31 : cp_logger_type
32 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
33 : cp_print_key_unit_nr
34 : USE input_constants, ONLY: do_full_density,&
35 : do_spin_density
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_string_length,&
40 : dp
41 : USE mathconstants, ONLY: pi
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE particle_types, ONLY: particle_type
44 : USE pw_env_types, ONLY: pw_env_get,&
45 : pw_env_type
46 : USE pw_methods, ONLY: pw_axpy,&
47 : pw_copy,&
48 : pw_transfer
49 : USE pw_pool_types, ONLY: pw_pool_type
50 : USE pw_types, ONLY: pw_c1d_gs_type,&
51 : pw_r3d_rs_type
52 : USE qs_charges_types, ONLY: qs_charges_type
53 : USE qs_environment_types, ONLY: get_qs_env,&
54 : qs_environment_type
55 : USE qs_kind_types, ONLY: qs_kind_type
56 : USE qs_rho_types, ONLY: qs_rho_get,&
57 : qs_rho_type
58 : #include "./base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 : PRIVATE
62 :
63 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_util'
65 : PUBLIC :: get_ddapc, &
66 : restraint_functional_potential, &
67 : modify_hartree_pot, &
68 : cp_ddapc_init
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief Initialize the cp_ddapc_environment
74 : !> \param qs_env ...
75 : !> \par History
76 : !> 08.2005 created [tlaino]
77 : !> \author Teodoro Laino
78 : ! **************************************************************************************************
79 21868 : SUBROUTINE cp_ddapc_init(qs_env)
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 :
82 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_ddapc_init'
83 :
84 : INTEGER :: handle, i, iw, iw2, n_rep_val, num_gauss
85 : LOGICAL :: allocate_ddapc_env, unimplemented
86 : REAL(KIND=dp) :: gcut, pfact, rcmin, Vol
87 21868 : REAL(KIND=dp), DIMENSION(:), POINTER :: inp_radii, radii
88 : TYPE(cell_type), POINTER :: cell, super_cell
89 : TYPE(cp_logger_type), POINTER :: logger
90 : TYPE(dft_control_type), POINTER :: dft_control
91 : TYPE(mp_para_env_type), POINTER :: para_env
92 21868 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
93 : TYPE(pw_c1d_gs_type) :: rho_tot_g
94 : TYPE(pw_env_type), POINTER :: pw_env
95 : TYPE(pw_pool_type), POINTER :: auxbas_pool
96 : TYPE(qs_charges_type), POINTER :: qs_charges
97 : TYPE(qs_rho_type), POINTER :: rho
98 : TYPE(section_vals_type), POINTER :: density_fit_section
99 :
100 21868 : CALL timeset(routineN, handle)
101 21868 : logger => cp_get_default_logger()
102 21868 : NULLIFY (dft_control, rho, pw_env, &
103 21868 : radii, inp_radii, particle_set, qs_charges, para_env)
104 :
105 21868 : CALL get_qs_env(qs_env, dft_control=dft_control)
106 : allocate_ddapc_env = qs_env%cp_ddapc_ewald%do_solvation .OR. &
107 : qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
108 : qs_env%cp_ddapc_ewald%do_decoupling .OR. &
109 21868 : qs_env%cp_ddapc_ewald%do_restraint
110 : unimplemented = dft_control%qs_control%semi_empirical .OR. &
111 : dft_control%qs_control%dftb .OR. &
112 21868 : dft_control%qs_control%xtb
113 21868 : IF (allocate_ddapc_env .AND. unimplemented) THEN
114 0 : CPABORT("DDAP charges work only with GPW/GAPW code.")
115 : END IF
116 : allocate_ddapc_env = allocate_ddapc_env .OR. &
117 21868 : qs_env%cp_ddapc_ewald%do_property
118 21868 : allocate_ddapc_env = allocate_ddapc_env .AND. (.NOT. unimplemented)
119 21868 : IF (allocate_ddapc_env) THEN
120 : CALL get_qs_env(qs_env=qs_env, &
121 : dft_control=dft_control, &
122 : rho=rho, &
123 : pw_env=pw_env, &
124 : qs_charges=qs_charges, &
125 : particle_set=particle_set, &
126 : cell=cell, &
127 : super_cell=super_cell, &
128 248 : para_env=para_env)
129 248 : density_fit_section => section_vals_get_subs_vals(qs_env%input, "DFT%DENSITY_FITTING")
130 : iw = cp_print_key_unit_nr(logger, density_fit_section, &
131 248 : "PROGRAM_RUN_INFO", ".FitCharge")
132 248 : IF (iw > 0) THEN
133 41 : WRITE (iw, '(/,A)') " Initializing the DDAPC Environment"
134 : END IF
135 248 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pool)
136 248 : CALL auxbas_pool%create_pw(rho_tot_g)
137 248 : Vol = rho_tot_g%pw_grid%vol
138 : !
139 : ! Get Input Parameters
140 : !
141 248 : CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val)
142 248 : IF (n_rep_val /= 0) THEN
143 0 : CALL section_vals_val_get(density_fit_section, "RADII", r_vals=inp_radii)
144 0 : num_gauss = SIZE(inp_radii)
145 0 : ALLOCATE (radii(num_gauss))
146 0 : radii = inp_radii
147 : ELSE
148 248 : CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss)
149 248 : CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin)
150 248 : CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact)
151 744 : ALLOCATE (radii(num_gauss))
152 1192 : DO i = 1, num_gauss
153 944 : radii(i) = rcmin*pfact**(i - 1)
154 : END DO
155 : END IF
156 248 : CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
157 : ! Create DDAPC environment
158 : iw2 = cp_print_key_unit_nr(logger, density_fit_section, &
159 248 : "PROGRAM_RUN_INFO/CONDITION_NUMBER", ".FitCharge")
160 : ! Initialization of the cp_ddapc_env and of the cp_ddapc_ewald environment
161 : !NB pass qs_env%para_env for parallelization of ewald_ddapc_pot()
162 248 : ALLOCATE (qs_env%cp_ddapc_env)
163 : CALL cp_ddapc_create(para_env, &
164 : qs_env%cp_ddapc_env, &
165 : qs_env%cp_ddapc_ewald, &
166 : particle_set, &
167 : radii, &
168 : cell, &
169 : super_cell, &
170 : rho_tot_g, &
171 : gcut, &
172 : iw2, &
173 : Vol, &
174 248 : qs_env%input)
175 : CALL cp_print_key_finished_output(iw2, logger, density_fit_section, &
176 248 : "PROGRAM_RUN_INFO/CONDITION_NUMBER")
177 248 : DEALLOCATE (radii)
178 248 : CALL auxbas_pool%give_back_pw(rho_tot_g)
179 : END IF
180 21868 : CALL timestop(handle)
181 21868 : END SUBROUTINE cp_ddapc_init
182 :
183 : ! **************************************************************************************************
184 : !> \brief Computes the Density Derived Atomic Point Charges
185 : !> \param qs_env ...
186 : !> \param calc_force ...
187 : !> \param density_fit_section ...
188 : !> \param density_type ...
189 : !> \param qout1 ...
190 : !> \param qout2 ...
191 : !> \param out_radii ...
192 : !> \param dq_out ...
193 : !> \param ext_rho_tot_g ...
194 : !> \param Itype_of_density ...
195 : !> \param iwc ...
196 : !> \par History
197 : !> 08.2005 created [tlaino]
198 : !> \author Teodoro Laino
199 : ! **************************************************************************************************
200 1764 : RECURSIVE SUBROUTINE get_ddapc(qs_env, calc_force, density_fit_section, &
201 : density_type, qout1, qout2, out_radii, dq_out, ext_rho_tot_g, &
202 : Itype_of_density, iwc)
203 : TYPE(qs_environment_type), POINTER :: qs_env
204 : LOGICAL, INTENT(in), OPTIONAL :: calc_force
205 : TYPE(section_vals_type), POINTER :: density_fit_section
206 : INTEGER, OPTIONAL :: density_type
207 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: qout1, qout2, out_radii
208 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
209 : POINTER :: dq_out
210 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: ext_rho_tot_g
211 : CHARACTER(LEN=*), OPTIONAL :: Itype_of_density
212 : INTEGER, INTENT(IN), OPTIONAL :: iwc
213 :
214 : CHARACTER(len=*), PARAMETER :: routineN = 'get_ddapc'
215 :
216 : CHARACTER(LEN=default_string_length) :: type_of_density
217 : INTEGER :: handle, handle2, handle3, i, ii, &
218 : iparticle, iparticle0, ispin, iw, j, &
219 : myid, n_rep_val, ndim, nparticles, &
220 : num_gauss, pmax, pmin
221 : LOGICAL :: need_f
222 : REAL(KIND=dp) :: c1, c3, ch_dens, gcut, pfact, rcmin, Vol
223 1764 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: AmI_bv, AmI_cv, bv, cv, cvT_AmI, &
224 1764 : cvT_AmI_dAmj, dAmj_qv, qtot, qv
225 1764 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dbv, g_dot_rvec_cos, g_dot_rvec_sin
226 1764 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dAm, dqv, tv
227 1764 : REAL(KIND=dp), DIMENSION(:), POINTER :: inp_radii, radii
228 : TYPE(cell_type), POINTER :: cell, super_cell
229 : TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env
230 : TYPE(cp_logger_type), POINTER :: logger
231 : TYPE(dft_control_type), POINTER :: dft_control
232 1764 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
233 : TYPE(pw_c1d_gs_type) :: rho_tot_g
234 1764 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
235 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
236 : TYPE(pw_env_type), POINTER :: pw_env
237 : TYPE(pw_pool_type), POINTER :: auxbas_pool
238 1764 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
239 : TYPE(qs_charges_type), POINTER :: qs_charges
240 1764 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
241 : TYPE(qs_rho_type), POINTER :: rho
242 :
243 : !NB variables for doing build_der_A_matrix_rows in blocks
244 : !NB refactor math in inner loop - no need for dqv0
245 : !!NB refactor math in inner loop - new temporaries
246 :
247 : EXTERNAL dgemv, dgemm
248 :
249 1764 : CALL timeset(routineN, handle)
250 1764 : need_f = .FALSE.
251 1764 : IF (PRESENT(calc_force)) need_f = calc_force
252 1764 : logger => cp_get_default_logger()
253 1764 : NULLIFY (dft_control, rho, rho_core, rho0_s_gs, rhoz_cneo_s_gs, pw_env, rho_g, rho_r, &
254 1764 : radii, inp_radii, particle_set, qs_kind_set, qs_charges, cp_ddapc_env)
255 : CALL get_qs_env(qs_env=qs_env, &
256 : dft_control=dft_control, &
257 : rho=rho, &
258 : rho_core=rho_core, &
259 : rho0_s_gs=rho0_s_gs, &
260 : rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
261 : pw_env=pw_env, &
262 : qs_charges=qs_charges, &
263 : particle_set=particle_set, &
264 : qs_kind_set=qs_kind_set, &
265 : cell=cell, &
266 1764 : super_cell=super_cell)
267 :
268 1764 : CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
269 :
270 1764 : IF (PRESENT(iwc)) THEN
271 938 : iw = iwc
272 : ELSE
273 : iw = cp_print_key_unit_nr(logger, density_fit_section, &
274 826 : "PROGRAM_RUN_INFO", ".FitCharge")
275 : END IF
276 : CALL pw_env_get(pw_env=pw_env, &
277 1764 : auxbas_pw_pool=auxbas_pool)
278 1764 : CALL auxbas_pool%create_pw(rho_tot_g)
279 1764 : IF (PRESENT(ext_rho_tot_g)) THEN
280 : ! If provided use the input density in g-space
281 826 : CALL pw_transfer(ext_rho_tot_g, rho_tot_g)
282 826 : type_of_density = Itype_of_density
283 : ELSE
284 938 : IF (PRESENT(density_type)) THEN
285 836 : myid = density_type
286 : ELSE
287 : CALL section_vals_val_get(qs_env%input, &
288 102 : "PROPERTIES%FIT_CHARGE%TYPE_OF_DENSITY", i_val=myid)
289 : END IF
290 766 : SELECT CASE (myid)
291 : CASE (do_full_density)
292 : ! Otherwise build the total QS density (electron+nuclei) in G-space
293 766 : IF (dft_control%qs_control%gapw) THEN
294 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
295 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
296 : END IF
297 0 : CALL pw_transfer(rho0_s_gs, rho_tot_g)
298 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
299 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
300 : END IF
301 : ELSE
302 766 : CALL pw_transfer(rho_core, rho_tot_g)
303 : END IF
304 2172 : DO ispin = 1, SIZE(rho_g)
305 2172 : CALL pw_axpy(rho_g(ispin), rho_tot_g)
306 : END DO
307 766 : type_of_density = "FULL DENSITY"
308 : CASE (do_spin_density)
309 172 : CALL pw_copy(rho_g(1), rho_tot_g)
310 172 : CALL pw_axpy(rho_g(2), rho_tot_g, alpha=-1._dp)
311 1110 : type_of_density = "SPIN DENSITY"
312 : END SELECT
313 : END IF
314 1764 : Vol = rho_r(1)%pw_grid%vol
315 1764 : ch_dens = 0.0_dp
316 : ! should use pw_integrate
317 1764 : IF (rho_tot_g%pw_grid%have_g0) ch_dens = REAL(rho_tot_g%array(1), KIND=dp)
318 1764 : CALL logger%para_env%sum(ch_dens)
319 : !
320 : ! Get Input Parameters
321 : !
322 1764 : CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val)
323 1764 : IF (n_rep_val /= 0) THEN
324 0 : CALL section_vals_val_get(density_fit_section, "RADII", r_vals=inp_radii)
325 0 : num_gauss = SIZE(inp_radii)
326 0 : ALLOCATE (radii(num_gauss))
327 0 : radii = inp_radii
328 : ELSE
329 1764 : CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss)
330 1764 : CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin)
331 1764 : CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact)
332 5292 : ALLOCATE (radii(num_gauss))
333 7764 : DO i = 1, num_gauss
334 6000 : radii(i) = rcmin*pfact**(i - 1)
335 : END DO
336 : END IF
337 1764 : IF (PRESENT(out_radii)) THEN
338 1662 : IF (ASSOCIATED(out_radii)) THEN
339 0 : DEALLOCATE (out_radii)
340 : END IF
341 4986 : ALLOCATE (out_radii(SIZE(radii)))
342 9522 : out_radii = radii
343 : END IF
344 1764 : CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
345 1764 : cp_ddapc_env => qs_env%cp_ddapc_env
346 : !
347 : ! Start with the linear system
348 : !
349 1764 : ndim = SIZE(particle_set)*SIZE(radii)
350 5292 : ALLOCATE (bv(ndim))
351 3528 : ALLOCATE (qv(ndim))
352 5292 : ALLOCATE (qtot(SIZE(particle_set)))
353 3528 : ALLOCATE (cv(ndim))
354 1764 : CALL timeset(routineN//"-charges", handle2)
355 14424 : bv(:) = 0.0_dp
356 14424 : cv(:) = 1.0_dp/Vol
357 : CALL build_b_vector(bv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
358 14424 : particle_set, radii, rho_tot_g, gcut); bv(:) = bv(:)/Vol
359 1764 : CALL rho_tot_g%pw_grid%para%group%sum(bv)
360 722700 : c1 = DOT_PRODUCT(cv, MATMUL(cp_ddapc_env%AmI, bv)) - ch_dens
361 1764 : c1 = c1/cp_ddapc_env%c0
362 508884 : qv(:) = -MATMUL(cp_ddapc_env%AmI, (bv - c1*cv))
363 6776 : j = 0
364 6776 : qtot = 0.0_dp
365 6776 : DO i = 1, ndim, num_gauss
366 5012 : j = j + 1
367 19436 : DO ii = 1, num_gauss
368 17672 : qtot(j) = qtot(j) + qv((i - 1) + ii)
369 : END DO
370 : END DO
371 1764 : IF (PRESENT(qout1)) THEN
372 1662 : IF (ASSOCIATED(qout1)) THEN
373 0 : CPASSERT(SIZE(qout1) == SIZE(qv))
374 : ELSE
375 4986 : ALLOCATE (qout1(SIZE(qv)))
376 : END IF
377 13188 : qout1 = qv
378 : END IF
379 1764 : IF (PRESENT(qout2)) THEN
380 0 : IF (ASSOCIATED(qout2)) THEN
381 0 : CPASSERT(SIZE(qout2) == SIZE(qtot))
382 : ELSE
383 0 : ALLOCATE (qout2(SIZE(qtot)))
384 : END IF
385 0 : qout2 = qtot
386 : END IF
387 : CALL print_atomic_charges(particle_set, qs_kind_set, iw, title=" DDAP "// &
388 1764 : TRIM(type_of_density)//" charges:", atomic_charges=qtot)
389 1764 : CALL timestop(handle2)
390 : !
391 : ! If requested evaluate also the correction to derivatives due to Pulay Forces
392 : !
393 1764 : IF (need_f) THEN
394 140 : CALL timeset(routineN//"-forces", handle3)
395 140 : IF (iw > 0) THEN
396 18 : WRITE (iw, '(T3,A)') " Evaluating DDAPC atomic derivatives .."
397 : END IF
398 700 : ALLOCATE (dAm(ndim, ndim, 3))
399 420 : ALLOCATE (dbv(ndim, 3))
400 700 : ALLOCATE (dqv(ndim, SIZE(particle_set), 3))
401 : !NB refactor math in inner loop - no more dqv0, but new temporaries instead
402 280 : ALLOCATE (cvT_AmI(ndim))
403 280 : ALLOCATE (cvT_AmI_dAmj(ndim))
404 420 : ALLOCATE (tv(ndim, SIZE(particle_set), 3))
405 280 : ALLOCATE (AmI_cv(ndim))
406 25772 : cvT_AmI(:) = MATMUL(cv, cp_ddapc_env%AmI)
407 25772 : AmI_cv(:) = MATMUL(cp_ddapc_env%AmI, cv)
408 280 : ALLOCATE (dAmj_qv(ndim))
409 280 : ALLOCATE (AmI_bv(ndim))
410 37508 : AmI_bv(:) = MATMUL(cp_ddapc_env%AmI, bv)
411 :
412 : !NB call routine to precompute sin(g.r) and cos(g.r),
413 : ! so it doesn't have to be done for each r_i-r_j pair in build_der_A_matrix_rows()
414 140 : CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
415 : !NB do build_der_A_matrix_rows in blocks, for more efficient use of DGEMM
416 : #define NPSET 100
417 280 : DO iparticle0 = 1, SIZE(particle_set), NPSET
418 140 : nparticles = MIN(NPSET, SIZE(particle_set) - iparticle0 + 1)
419 : !NB each dAm is supposed to have one block of rows and one block of columns
420 : !NB for derivatives with respect to each atom. build_der_A_matrix_rows()
421 : !NB just returns rows, since dAm is symmetric, and missing columns can be
422 : !NB reconstructed with a simple transpose, as below
423 : CALL build_der_A_matrix_rows(dAm, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
424 : particle_set, radii, rho_tot_g, gcut, iparticle0, &
425 140 : nparticles, g_dot_rvec_sin, g_dot_rvec_cos)
426 : !NB no more reduction of dbv and dAm - instead we go through with each node's contribution
427 : !NB and reduce resulting charges/forces once, at the end. Intermediate speedup can be
428 : !NB had by reducing dqv after the inner loop, and then other routines don't need to know
429 : !NB that contributions to dqv are distributed over the nodes.
430 : !NB also get rid of zeroing of dAm and division by Vol**2 - it's slow, and can be done
431 : !NB more quickly later, to a scalar or vector rather than a matrix
432 676 : DO iparticle = iparticle0, iparticle0 + nparticles - 1
433 : IF (debug_this_module) THEN
434 : CALL debug_der_A_matrix(dAm, particle_set, radii, rho_tot_g, &
435 : gcut, iparticle, Vol, qs_env)
436 : cp_ddapc_env => qs_env%cp_ddapc_env
437 : END IF
438 13572 : dbv(:, :) = 0.0_dp
439 : CALL build_der_b_vector(dbv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
440 396 : particle_set, radii, rho_tot_g, gcut, iparticle)
441 13572 : dbv(:, :) = dbv(:, :)/Vol
442 : IF (debug_this_module) THEN
443 : CALL debug_der_b_vector(dbv, particle_set, radii, rho_tot_g, &
444 : gcut, iparticle, Vol, qs_env)
445 : cp_ddapc_env => qs_env%cp_ddapc_env
446 : END IF
447 1584 : DO j = 1, 3
448 : !NB dAmj is actually pretty sparse - one block of cols + one block of rows - use this here:
449 1188 : pmin = (iparticle - 1)*SIZE(radii) + 1
450 1188 : pmax = iparticle*SIZE(radii)
451 : !NB multiply by block of columns that aren't explicitly in dAm, but can be reconstructured
452 : !NB as transpose of relevant block of rows
453 1188 : IF (pmin > 1) THEN
454 768 : dAmj_qv(:pmin - 1) = MATMUL(TRANSPOSE(dAm(pmin:pmax, :pmin - 1, j)), qv(pmin:pmax))
455 768 : cvT_AmI_dAmj(:pmin - 1) = MATMUL(TRANSPOSE(dAm(pmin:pmax, :pmin - 1, j)), cvT_AmI(pmin:pmax))
456 : END IF
457 : !NB multiply by block of rows that are explicitly in dAm
458 51624 : dAmj_qv(pmin:pmax) = MATMUL(dAm(pmin:pmax, :, j), qv(:))
459 51624 : cvT_AmI_dAmj(pmin:pmax) = MATMUL(dAm(pmin:pmax, :, j), cvT_AmI(:))
460 : !NB multiply by block of columns that aren't explicitly in dAm, but can be reconstructured
461 : !NB as transpose of relevant block of rows
462 1188 : IF (pmax < SIZE(particle_set)*SIZE(radii)) THEN
463 768 : dAmj_qv(pmax + 1:) = MATMUL(TRANSPOSE(dAm(pmin:pmax, pmax + 1:, j)), qv(pmin:pmax))
464 768 : cvT_AmI_dAmj(pmax + 1:) = MATMUL(TRANSPOSE(dAm(pmin:pmax, pmax + 1:, j)), cvT_AmI(pmin:pmax))
465 : END IF
466 13176 : dAmj_qv(:) = dAmj_qv(:)/(Vol*Vol)
467 13176 : cvT_AmI_dAmj(:) = cvT_AmI_dAmj(:)/(Vol*Vol)
468 37152 : c3 = DOT_PRODUCT(cvT_AmI_dAmj, AmI_bv) - DOT_PRODUCT(cvT_AmI, dbv(:, j)) - c1*DOT_PRODUCT(cvT_AmI_dAmj, AmI_cv)
469 13572 : tv(:, iparticle, j) = -(dAmj_qv(:) + dbv(:, j) + c3/cp_ddapc_env%c0*cv)
470 : END DO ! j
471 : !NB zero relevant parts of dAm here
472 48920 : dAm((iparticle - 1)*SIZE(radii) + 1:iparticle*SIZE(radii), :, :) = 0.0_dp
473 : !! dAm(:,(iparticle-1)*SIZE(radii)+1:iparticle*SIZE(radii),:) = 0.0_dp
474 : END DO ! iparticle
475 : END DO ! iparticle0
476 : !NB final part of refactoring of math - one dgemm is faster than many dgemv
477 : CALL dgemm('N', 'N', SIZE(dqv, 1), SIZE(dqv, 2)*SIZE(dqv, 3), SIZE(cp_ddapc_env%AmI, 2), 1.0_dp, &
478 140 : cp_ddapc_env%AmI, SIZE(cp_ddapc_env%AmI, 1), tv, SIZE(tv, 1), 0.0_dp, dqv, SIZE(dqv, 1))
479 : !NB deallocate g_dot_rvec_sin and g_dot_rvec_cos
480 140 : CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
481 : !NB moved reduction out to where dqv is used to compute
482 : !NB a force contribution (smaller array to reduce, just size(particle_set) x 3)
483 : !NB namely ewald_ddapc_force(), cp_decl_ddapc_forces(), restraint_functional_force()
484 140 : CPASSERT(PRESENT(dq_out))
485 140 : IF (.NOT. ASSOCIATED(dq_out)) THEN
486 700 : ALLOCATE (dq_out(SIZE(dqv, 1), SIZE(dqv, 2), SIZE(dqv, 3)))
487 : ELSE
488 0 : CPASSERT(SIZE(dqv, 1) == SIZE(dq_out, 1))
489 0 : CPASSERT(SIZE(dqv, 2) == SIZE(dq_out, 2))
490 0 : CPASSERT(SIZE(dqv, 3) == SIZE(dq_out, 3))
491 : END IF
492 13736 : dq_out = dqv
493 : IF (debug_this_module) THEN
494 : CALL debug_charge(dqv, qs_env, density_fit_section, &
495 : particle_set, radii, rho_tot_g, type_of_density)
496 : cp_ddapc_env => qs_env%cp_ddapc_env
497 : END IF
498 140 : DEALLOCATE (dqv)
499 140 : DEALLOCATE (dAm)
500 140 : DEALLOCATE (dbv)
501 : !NB deallocate new temporaries
502 140 : DEALLOCATE (cvT_AmI)
503 140 : DEALLOCATE (cvT_AmI_dAmj)
504 140 : DEALLOCATE (AmI_cv)
505 140 : DEALLOCATE (tv)
506 140 : DEALLOCATE (dAmj_qv)
507 140 : DEALLOCATE (AmI_bv)
508 280 : CALL timestop(handle3)
509 : END IF
510 : !
511 : ! End of charge fit
512 : !
513 1764 : DEALLOCATE (radii)
514 1764 : DEALLOCATE (bv)
515 1764 : DEALLOCATE (cv)
516 1764 : DEALLOCATE (qv)
517 1764 : DEALLOCATE (qtot)
518 1764 : IF (.NOT. PRESENT(iwc)) THEN
519 : CALL cp_print_key_finished_output(iw, logger, density_fit_section, &
520 826 : "PROGRAM_RUN_INFO")
521 : END IF
522 1764 : CALL auxbas_pool%give_back_pw(rho_tot_g)
523 1764 : CALL timestop(handle)
524 8820 : END SUBROUTINE get_ddapc
525 :
526 : ! **************************************************************************************************
527 : !> \brief modify hartree potential to handle restraints in DDAPC scheme
528 : !> \param v_hartree_gspace ...
529 : !> \param density_fit_section ...
530 : !> \param particle_set ...
531 : !> \param AmI ...
532 : !> \param radii ...
533 : !> \param charges ...
534 : !> \param ddapc_restraint_control ...
535 : !> \param energy_res ...
536 : !> \par History
537 : !> 02.2006 modified [Teo]
538 : ! **************************************************************************************************
539 836 : SUBROUTINE restraint_functional_potential(v_hartree_gspace, &
540 : density_fit_section, particle_set, AmI, radii, charges, &
541 : ddapc_restraint_control, energy_res)
542 : TYPE(pw_c1d_gs_type), INTENT(IN) :: v_hartree_gspace
543 : TYPE(section_vals_type), POINTER :: density_fit_section
544 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
545 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: AmI
546 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii, charges
547 : TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control
548 : REAL(KIND=dp), INTENT(INOUT) :: energy_res
549 :
550 : CHARACTER(len=*), PARAMETER :: routineN = 'restraint_functional_potential'
551 :
552 : COMPLEX(KIND=dp) :: g_corr, phase
553 : INTEGER :: handle, idim, ig, igauss, iparticle, &
554 : n_gauss
555 : REAL(KIND=dp) :: arg, fac, fac2, g2, gcut, gcut2, gfunc, &
556 : gvec(3), rc, rc2, rvec(3), sfac, Vol, w
557 836 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
558 :
559 836 : CALL timeset(routineN, handle)
560 836 : n_gauss = SIZE(radii)
561 2508 : ALLOCATE (cv(n_gauss*SIZE(particle_set)))
562 1672 : ALLOCATE (uv(n_gauss*SIZE(particle_set)))
563 4460 : uv = 0.0_dp
564 : CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
565 836 : charges, energy_res)
566 : !
567 836 : CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
568 836 : gcut2 = gcut*gcut
569 : ASSOCIATE (pw_grid => v_hartree_gspace%pw_grid)
570 836 : Vol = pw_grid%vol
571 4460 : cv = 1.0_dp/Vol
572 836 : sfac = -1.0_dp/Vol
573 58740 : fac = DOT_PRODUCT(cv, MATMUL(AmI, cv))
574 81420 : fac2 = DOT_PRODUCT(cv, MATMUL(AmI, uv))
575 4460 : cv(:) = uv - cv*fac2/fac
576 57068 : cv(:) = MATMUL(AmI, cv)
577 836 : IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
578 144832 : DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
579 143996 : g2 = pw_grid%gsq(ig)
580 143996 : w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
581 143996 : IF (g2 > gcut2) EXIT
582 572640 : gvec = pw_grid%g(:, ig)
583 143160 : g_corr = 0.0_dp
584 143160 : idim = 0
585 507936 : DO iparticle = 1, SIZE(particle_set)
586 1401888 : DO igauss = 1, SIZE(radii)
587 893952 : idim = idim + 1
588 893952 : rc = radii(igauss)
589 893952 : rc2 = rc*rc
590 3575808 : rvec = particle_set(iparticle)%r
591 3575808 : arg = DOT_PRODUCT(gvec, rvec)
592 893952 : phase = CMPLX(COS(arg), -SIN(arg), KIND=dp)
593 893952 : gfunc = EXP(-g2*rc2/4.0_dp)
594 1258728 : g_corr = g_corr + gfunc*cv(idim)*phase
595 : END DO
596 : END DO
597 143160 : g_corr = g_corr*w
598 143996 : v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/Vol
599 : END DO
600 : END ASSOCIATE
601 836 : CALL timestop(handle)
602 2508 : END SUBROUTINE restraint_functional_potential
603 :
604 : ! **************************************************************************************************
605 : !> \brief Modify the Hartree potential
606 : !> \param v_hartree_gspace ...
607 : !> \param density_fit_section ...
608 : !> \param particle_set ...
609 : !> \param M ...
610 : !> \param AmI ...
611 : !> \param radii ...
612 : !> \param charges ...
613 : !> \par History
614 : !> 08.2005 created [tlaino]
615 : !> \author Teodoro Laino
616 : ! **************************************************************************************************
617 826 : SUBROUTINE modify_hartree_pot(v_hartree_gspace, density_fit_section, &
618 : particle_set, M, AmI, radii, charges)
619 : TYPE(pw_c1d_gs_type), INTENT(IN) :: v_hartree_gspace
620 : TYPE(section_vals_type), POINTER :: density_fit_section
621 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
622 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M, AmI
623 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii, charges
624 :
625 : CHARACTER(len=*), PARAMETER :: routineN = 'modify_hartree_pot'
626 :
627 : COMPLEX(KIND=dp) :: g_corr, phase
628 : INTEGER :: handle, idim, ig, igauss, iparticle
629 : REAL(kind=dp) :: arg, fac, fac2, g2, gcut, gcut2, gfunc, &
630 : gvec(3), rc, rc2, rvec(3), sfac, Vol, w
631 826 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
632 :
633 826 : CALL timeset(routineN, handle)
634 826 : CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
635 826 : gcut2 = gcut*gcut
636 : ASSOCIATE (pw_grid => v_hartree_gspace%pw_grid)
637 826 : Vol = pw_grid%vol
638 2478 : ALLOCATE (cv(SIZE(M, 1)))
639 1652 : ALLOCATE (uv(SIZE(M, 1)))
640 8728 : cv = 1.0_dp/Vol
641 492964 : uv(:) = MATMUL(M, charges)
642 826 : sfac = -1.0_dp/Vol
643 343740 : fac = DOT_PRODUCT(cv, MATMUL(AmI, cv))
644 343740 : fac2 = DOT_PRODUCT(cv, MATMUL(AmI, uv))
645 8728 : cv(:) = uv - cv*fac2/fac
646 342088 : cv(:) = MATMUL(AmI, cv)
647 826 : IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
648 708568 : DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
649 707742 : g2 = pw_grid%gsq(ig)
650 707742 : w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
651 707742 : IF (g2 > gcut2) EXIT
652 2827664 : gvec = pw_grid%g(:, ig)
653 706916 : g_corr = 0.0_dp
654 706916 : idim = 0
655 3035658 : DO iparticle = 1, SIZE(particle_set)
656 10021884 : DO igauss = 1, SIZE(radii)
657 6986226 : idim = idim + 1
658 6986226 : rc = radii(igauss)
659 6986226 : rc2 = rc*rc
660 27944904 : rvec = particle_set(iparticle)%r
661 27944904 : arg = DOT_PRODUCT(gvec, rvec)
662 6986226 : phase = CMPLX(COS(arg), -SIN(arg), KIND=dp)
663 6986226 : gfunc = EXP(-g2*rc2/4.0_dp)
664 9314968 : g_corr = g_corr + gfunc*cv(idim)*phase
665 : END DO
666 : END DO
667 706916 : g_corr = g_corr*w
668 707742 : v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/Vol
669 : END DO
670 : END ASSOCIATE
671 826 : CALL timestop(handle)
672 1652 : END SUBROUTINE modify_hartree_pot
673 :
674 : ! **************************************************************************************************
675 : !> \brief To Debug the derivative of the B vector for the solution of the
676 : !> linear system
677 : !> \param dbv ...
678 : !> \param particle_set ...
679 : !> \param radii ...
680 : !> \param rho_tot_g ...
681 : !> \param gcut ...
682 : !> \param iparticle ...
683 : !> \param Vol ...
684 : !> \param qs_env ...
685 : !> \par History
686 : !> 08.2005 created [tlaino]
687 : !> \author Teodoro Laino
688 : ! **************************************************************************************************
689 0 : SUBROUTINE debug_der_b_vector(dbv, particle_set, radii, &
690 : rho_tot_g, gcut, iparticle, Vol, qs_env)
691 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: dbv
692 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
693 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
694 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
695 : REAL(KIND=dp), INTENT(IN) :: gcut
696 : INTEGER, INTENT(in) :: iparticle
697 : REAL(KIND=dp), INTENT(IN) :: Vol
698 : TYPE(qs_environment_type), POINTER :: qs_env
699 :
700 : CHARACTER(len=*), PARAMETER :: routineN = 'debug_der_b_vector'
701 :
702 : INTEGER :: handle, i, kk, ndim
703 : REAL(KIND=dp) :: dx, rvec(3), v0
704 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bv1, bv2, ddbv
705 : TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env
706 :
707 0 : NULLIFY (cp_ddapc_env)
708 0 : CALL timeset(routineN, handle)
709 0 : dx = 0.01_dp
710 0 : ndim = SIZE(particle_set)*SIZE(radii)
711 0 : ALLOCATE (bv1(ndim))
712 0 : ALLOCATE (bv2(ndim))
713 0 : ALLOCATE (ddbv(ndim))
714 0 : rvec = particle_set(iparticle)%r
715 0 : cp_ddapc_env => qs_env%cp_ddapc_env
716 0 : DO i = 1, 3
717 0 : bv1(:) = 0.0_dp
718 0 : bv2(:) = 0.0_dp
719 0 : particle_set(iparticle)%r(i) = rvec(i) + dx
720 : CALL build_b_vector(bv1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
721 0 : particle_set, radii, rho_tot_g, gcut); bv1(:) = bv1(:)/Vol
722 0 : CALL rho_tot_g%pw_grid%para%group%sum(bv1)
723 0 : particle_set(iparticle)%r(i) = rvec(i) - dx
724 : CALL build_b_vector(bv2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
725 0 : particle_set, radii, rho_tot_g, gcut); bv2(:) = bv2(:)/Vol
726 0 : CALL rho_tot_g%pw_grid%para%group%sum(bv2)
727 0 : ddbv(:) = (bv1(:) - bv2(:))/(2.0_dp*dx)
728 0 : DO kk = 1, SIZE(ddbv)
729 0 : IF (ddbv(kk) > 1.0E-8_dp) THEN
730 0 : v0 = ABS(dbv(kk, i) - ddbv(kk))/ddbv(kk)*100.0_dp
731 0 : WRITE (*, *) "Error % on B ::", v0
732 0 : IF (v0 > 0.1_dp) THEN
733 0 : WRITE (*, '(A,2I5,2F15.9)') "ERROR IN DERIVATIVE OF B VECTOR, IPARTICLE, ICOORD:", iparticle, i, &
734 0 : dbv(kk, i), ddbv(kk)
735 0 : CPABORT("")
736 : END IF
737 : END IF
738 : END DO
739 0 : particle_set(iparticle)%r = rvec
740 : END DO
741 0 : DEALLOCATE (bv1)
742 0 : DEALLOCATE (bv2)
743 0 : DEALLOCATE (ddbv)
744 0 : CALL timestop(handle)
745 0 : END SUBROUTINE debug_der_b_vector
746 :
747 : ! **************************************************************************************************
748 : !> \brief To Debug the derivative of the A matrix for the solution of the
749 : !> linear system
750 : !> \param dAm ...
751 : !> \param particle_set ...
752 : !> \param radii ...
753 : !> \param rho_tot_g ...
754 : !> \param gcut ...
755 : !> \param iparticle ...
756 : !> \param Vol ...
757 : !> \param qs_env ...
758 : !> \par History
759 : !> 08.2005 created [tlaino]
760 : !> \author Teodoro Laino
761 : ! **************************************************************************************************
762 0 : SUBROUTINE debug_der_A_matrix(dAm, particle_set, radii, &
763 : rho_tot_g, gcut, iparticle, Vol, qs_env)
764 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dAm
765 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
766 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
767 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
768 : REAL(KIND=dp), INTENT(IN) :: gcut
769 : INTEGER, INTENT(in) :: iparticle
770 : REAL(KIND=dp), INTENT(IN) :: Vol
771 : TYPE(qs_environment_type), POINTER :: qs_env
772 :
773 : CHARACTER(len=*), PARAMETER :: routineN = 'debug_der_A_matrix'
774 :
775 : INTEGER :: handle, i, kk, ll, ndim
776 : REAL(KIND=dp) :: dx, rvec(3), v0
777 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Am1, Am2, ddAm, g_dot_rvec_cos, &
778 0 : g_dot_rvec_sin
779 : TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env
780 :
781 : !NB new temporaries sin(g.r) and cos(g.r), as used in get_ddapc, to speed up build_der_A_matrix()
782 :
783 0 : NULLIFY (cp_ddapc_env)
784 0 : CALL timeset(routineN, handle)
785 0 : dx = 0.01_dp
786 0 : ndim = SIZE(particle_set)*SIZE(radii)
787 0 : ALLOCATE (Am1(ndim, ndim))
788 0 : ALLOCATE (Am2(ndim, ndim))
789 0 : ALLOCATE (ddAm(ndim, ndim))
790 0 : rvec = particle_set(iparticle)%r
791 0 : cp_ddapc_env => qs_env%cp_ddapc_env
792 0 : CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
793 0 : DO i = 1, 3
794 0 : Am1 = 0.0_dp
795 0 : Am2 = 0.0_dp
796 0 : particle_set(iparticle)%r(i) = rvec(i) + dx
797 : CALL build_A_matrix(Am1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
798 0 : particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
799 0 : Am1(:, :) = Am1(:, :)/(Vol*Vol)
800 0 : CALL rho_tot_g%pw_grid%para%group%sum(Am1)
801 0 : particle_set(iparticle)%r(i) = rvec(i) - dx
802 : CALL build_A_matrix(Am2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
803 0 : particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
804 0 : Am2(:, :) = Am2(:, :)/(Vol*Vol)
805 0 : CALL rho_tot_g%pw_grid%para%group%sum(Am2)
806 0 : ddAm(:, :) = (Am1 - Am2)/(2.0_dp*dx)
807 0 : DO kk = 1, SIZE(ddAm, 1)
808 0 : DO ll = 1, SIZE(ddAm, 2)
809 0 : IF (ddAm(kk, ll) > 1.0E-8_dp) THEN
810 0 : v0 = ABS(dAm(kk, ll, i) - ddAm(kk, ll))/ddAm(kk, ll)*100.0_dp
811 0 : WRITE (*, *) "Error % on A ::", v0, Am1(kk, ll), Am2(kk, ll), iparticle, i, kk, ll
812 0 : IF (v0 > 0.1_dp) THEN
813 0 : WRITE (*, '(A,4I5,2F15.9)') "ERROR IN DERIVATIVE OF A MATRIX, IPARTICLE, ICOORD:", iparticle, i, kk, ll, &
814 0 : dAm(kk, ll, i), ddAm(kk, ll)
815 0 : CPABORT("")
816 : END IF
817 : END IF
818 : END DO
819 : END DO
820 0 : particle_set(iparticle)%r = rvec
821 : END DO
822 0 : CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
823 0 : DEALLOCATE (Am1)
824 0 : DEALLOCATE (Am2)
825 0 : DEALLOCATE (ddAm)
826 0 : CALL timestop(handle)
827 0 : END SUBROUTINE debug_der_A_matrix
828 :
829 : ! **************************************************************************************************
830 : !> \brief To Debug the fitted charges
831 : !> \param dqv ...
832 : !> \param qs_env ...
833 : !> \param density_fit_section ...
834 : !> \param particle_set ...
835 : !> \param radii ...
836 : !> \param rho_tot_g ...
837 : !> \param type_of_density ...
838 : !> \par History
839 : !> 08.2005 created [tlaino]
840 : !> \author Teodoro Laino
841 : ! **************************************************************************************************
842 0 : SUBROUTINE debug_charge(dqv, qs_env, density_fit_section, &
843 : particle_set, radii, rho_tot_g, type_of_density)
844 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dqv
845 : TYPE(qs_environment_type), POINTER :: qs_env
846 : TYPE(section_vals_type), POINTER :: density_fit_section
847 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
848 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
849 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
850 : CHARACTER(LEN=*) :: type_of_density
851 :
852 : CHARACTER(len=*), PARAMETER :: routineN = 'debug_charge'
853 :
854 : INTEGER :: handle, i, iparticle, kk, ndim
855 : REAL(KIND=dp) :: dx, rvec(3)
856 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ddqv
857 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: qtot1, qtot2
858 :
859 0 : CALL timeset(routineN, handle)
860 0 : WRITE (*, *) "DEBUG_CHARGE_ROUTINE"
861 0 : ndim = SIZE(particle_set)*SIZE(radii)
862 0 : NULLIFY (qtot1, qtot2)
863 0 : ALLOCATE (qtot1(ndim))
864 0 : ALLOCATE (qtot2(ndim))
865 0 : ALLOCATE (ddqv(ndim))
866 : !
867 : dx = 0.001_dp
868 0 : DO iparticle = 1, SIZE(particle_set)
869 0 : rvec = particle_set(iparticle)%r
870 0 : DO i = 1, 3
871 0 : particle_set(iparticle)%r(i) = rvec(i) + dx
872 : CALL get_ddapc(qs_env, .FALSE., density_fit_section, qout1=qtot1, &
873 0 : ext_rho_tot_g=rho_tot_g, Itype_of_density=type_of_density)
874 0 : particle_set(iparticle)%r(i) = rvec(i) - dx
875 : CALL get_ddapc(qs_env, .FALSE., density_fit_section, qout1=qtot2, &
876 0 : ext_rho_tot_g=rho_tot_g, Itype_of_density=type_of_density)
877 0 : ddqv(:) = (qtot1 - qtot2)/(2.0_dp*dx)
878 0 : DO kk = 1, SIZE(qtot1) - 1, SIZE(radii)
879 0 : IF (ANY(ddqv(kk:kk + 2) > 1.0E-8_dp)) THEN
880 0 : WRITE (*, '(A,2F12.6,F12.2)') "Error :", SUM(dqv(kk:kk + 2, iparticle, i)), SUM(ddqv(kk:kk + 2)), &
881 0 : ABS((SUM(ddqv(kk:kk + 2)) - SUM(dqv(kk:kk + 2, iparticle, i)))/SUM(ddqv(kk:kk + 2))*100.0_dp)
882 : END IF
883 : END DO
884 0 : particle_set(iparticle)%r = rvec
885 : END DO
886 : END DO
887 : !
888 0 : DEALLOCATE (qtot1)
889 0 : DEALLOCATE (qtot2)
890 0 : DEALLOCATE (ddqv)
891 0 : CALL timestop(handle)
892 0 : END SUBROUTINE debug_charge
893 :
894 8514 : END MODULE cp_ddapc_util
|
