Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 19548 : 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 19548 : 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 19548 : 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 19548 : CALL timeset(routineN, handle)
101 19548 : logger => cp_get_default_logger()
102 19548 : NULLIFY (dft_control, rho, pw_env, &
103 19548 : radii, inp_radii, particle_set, qs_charges, para_env)
104 :
105 19548 : 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 19548 : qs_env%cp_ddapc_ewald%do_restraint
110 : unimplemented = dft_control%qs_control%semi_empirical .OR. &
111 : dft_control%qs_control%dftb .OR. &
112 19548 : dft_control%qs_control%xtb
113 19548 : 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 19548 : qs_env%cp_ddapc_ewald%do_property
118 19548 : allocate_ddapc_env = allocate_ddapc_env .AND. (.NOT. unimplemented)
119 19548 : 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 246 : para_env=para_env)
129 246 : 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 246 : "PROGRAM_RUN_INFO", ".FitCharge")
132 246 : IF (iw > 0) THEN
133 41 : WRITE (iw, '(/,A)') " Initializing the DDAPC Environment"
134 : END IF
135 246 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pool)
136 246 : CALL auxbas_pool%create_pw(rho_tot_g)
137 246 : Vol = rho_tot_g%pw_grid%vol
138 : !
139 : ! Get Input Parameters
140 : !
141 246 : CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val)
142 246 : 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 246 : CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss)
149 246 : CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin)
150 246 : CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact)
151 738 : ALLOCATE (radii(num_gauss))
152 1182 : DO i = 1, num_gauss
153 936 : radii(i) = rcmin*pfact**(i - 1)
154 : END DO
155 : END IF
156 246 : 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 246 : "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 246 : 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 246 : qs_env%input)
175 : CALL cp_print_key_finished_output(iw2, logger, density_fit_section, &
176 246 : "PROGRAM_RUN_INFO/CONDITION_NUMBER")
177 246 : DEALLOCATE (radii)
178 246 : CALL auxbas_pool%give_back_pw(rho_tot_g)
179 : END IF
180 19548 : CALL timestop(handle)
181 19548 : 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 1342 : 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 1342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: AmI_bv, AmI_cv, bv, cv, cvT_AmI, &
224 1342 : cvT_AmI_dAmj, dAmj_qv, qtot, qv
225 1342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dbv, g_dot_rvec_cos, g_dot_rvec_sin
226 1342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dAm, dqv, tv
227 1342 : 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 1342 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
233 : TYPE(pw_c1d_gs_type) :: rho_tot_g
234 1342 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
235 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
236 : TYPE(pw_env_type), POINTER :: pw_env
237 : TYPE(pw_pool_type), POINTER :: auxbas_pool
238 1342 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
239 : TYPE(qs_charges_type), POINTER :: qs_charges
240 1342 : 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 1342 : CALL timeset(routineN, handle)
250 1342 : need_f = .FALSE.
251 1342 : IF (PRESENT(calc_force)) need_f = calc_force
252 1342 : logger => cp_get_default_logger()
253 1342 : NULLIFY (dft_control, rho, rho_core, rho0_s_gs, pw_env, rho_g, rho_r, &
254 1342 : 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 : pw_env=pw_env, &
261 : qs_charges=qs_charges, &
262 : particle_set=particle_set, &
263 : qs_kind_set=qs_kind_set, &
264 : cell=cell, &
265 1342 : super_cell=super_cell)
266 :
267 1342 : CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
268 :
269 1342 : IF (PRESENT(iwc)) THEN
270 602 : iw = iwc
271 : ELSE
272 : iw = cp_print_key_unit_nr(logger, density_fit_section, &
273 740 : "PROGRAM_RUN_INFO", ".FitCharge")
274 : END IF
275 : CALL pw_env_get(pw_env=pw_env, &
276 1342 : auxbas_pw_pool=auxbas_pool)
277 1342 : CALL auxbas_pool%create_pw(rho_tot_g)
278 1342 : IF (PRESENT(ext_rho_tot_g)) THEN
279 : ! If provided use the input density in g-space
280 740 : CALL pw_transfer(ext_rho_tot_g, rho_tot_g)
281 740 : type_of_density = Itype_of_density
282 : ELSE
283 602 : IF (PRESENT(density_type)) THEN
284 500 : myid = density_type
285 : ELSE
286 : CALL section_vals_val_get(qs_env%input, &
287 102 : "PROPERTIES%FIT_CHARGE%TYPE_OF_DENSITY", i_val=myid)
288 : END IF
289 502 : SELECT CASE (myid)
290 : CASE (do_full_density)
291 : ! Otherwise build the total QS density (electron+nuclei) in G-space
292 502 : IF (dft_control%qs_control%gapw) THEN
293 0 : CALL pw_transfer(rho0_s_gs, rho_tot_g)
294 : ELSE
295 502 : CALL pw_transfer(rho_core, rho_tot_g)
296 : END IF
297 1380 : DO ispin = 1, SIZE(rho_g)
298 1380 : CALL pw_axpy(rho_g(ispin), rho_tot_g)
299 : END DO
300 502 : type_of_density = "FULL DENSITY"
301 : CASE (do_spin_density)
302 100 : CALL pw_copy(rho_g(1), rho_tot_g)
303 100 : CALL pw_axpy(rho_g(2), rho_tot_g, alpha=-1._dp)
304 702 : type_of_density = "SPIN DENSITY"
305 : END SELECT
306 : END IF
307 1342 : Vol = rho_r(1)%pw_grid%vol
308 1342 : ch_dens = 0.0_dp
309 : ! should use pw_integrate
310 1342 : IF (rho_tot_g%pw_grid%have_g0) ch_dens = REAL(rho_tot_g%array(1), KIND=dp)
311 1342 : CALL logger%para_env%sum(ch_dens)
312 : !
313 : ! Get Input Parameters
314 : !
315 1342 : CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val)
316 1342 : IF (n_rep_val /= 0) THEN
317 0 : CALL section_vals_val_get(density_fit_section, "RADII", r_vals=inp_radii)
318 0 : num_gauss = SIZE(inp_radii)
319 0 : ALLOCATE (radii(num_gauss))
320 0 : radii = inp_radii
321 : ELSE
322 1342 : CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss)
323 1342 : CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin)
324 1342 : CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact)
325 4026 : ALLOCATE (radii(num_gauss))
326 6086 : DO i = 1, num_gauss
327 4744 : radii(i) = rcmin*pfact**(i - 1)
328 : END DO
329 : END IF
330 1342 : IF (PRESENT(out_radii)) THEN
331 1240 : IF (ASSOCIATED(out_radii)) THEN
332 0 : DEALLOCATE (out_radii)
333 : END IF
334 3720 : ALLOCATE (out_radii(SIZE(radii)))
335 7432 : out_radii = radii
336 : END IF
337 1342 : CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
338 1342 : cp_ddapc_env => qs_env%cp_ddapc_env
339 : !
340 : ! Start with the linear system
341 : !
342 1342 : ndim = SIZE(particle_set)*SIZE(radii)
343 4026 : ALLOCATE (bv(ndim))
344 2684 : ALLOCATE (qv(ndim))
345 4026 : ALLOCATE (qtot(SIZE(particle_set)))
346 2684 : ALLOCATE (cv(ndim))
347 1342 : CALL timeset(routineN//"-charges", handle2)
348 11500 : bv(:) = 0.0_dp
349 11500 : cv(:) = 1.0_dp/Vol
350 : CALL build_b_vector(bv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
351 11500 : particle_set, radii, rho_tot_g, gcut); bv(:) = bv(:)/Vol
352 1342 : CALL rho_tot_g%pw_grid%para%group%sum(bv)
353 632226 : c1 = DOT_PRODUCT(cv, MATMUL(cp_ddapc_env%AmI, bv)) - ch_dens
354 1342 : c1 = c1/cp_ddapc_env%c0
355 443142 : qv(:) = -MATMUL(cp_ddapc_env%AmI, (bv - c1*cv))
356 5196 : j = 0
357 5196 : qtot = 0.0_dp
358 5196 : DO i = 1, ndim, num_gauss
359 3854 : j = j + 1
360 15354 : DO ii = 1, num_gauss
361 14012 : qtot(j) = qtot(j) + qv((i - 1) + ii)
362 : END DO
363 : END DO
364 1342 : IF (PRESENT(qout1)) THEN
365 1240 : IF (ASSOCIATED(qout1)) THEN
366 0 : CPASSERT(SIZE(qout1) == SIZE(qv))
367 : ELSE
368 3720 : ALLOCATE (qout1(SIZE(qv)))
369 : END IF
370 10264 : qout1 = qv
371 : END IF
372 1342 : IF (PRESENT(qout2)) THEN
373 0 : IF (ASSOCIATED(qout2)) THEN
374 0 : CPASSERT(SIZE(qout2) == SIZE(qtot))
375 : ELSE
376 0 : ALLOCATE (qout2(SIZE(qtot)))
377 : END IF
378 0 : qout2 = qtot
379 : END IF
380 : CALL print_atomic_charges(particle_set, qs_kind_set, iw, title=" DDAP "// &
381 1342 : TRIM(type_of_density)//" charges:", atomic_charges=qtot)
382 1342 : CALL timestop(handle2)
383 : !
384 : ! If requested evaluate also the correction to derivatives due to Pulay Forces
385 : !
386 1342 : IF (need_f) THEN
387 140 : CALL timeset(routineN//"-forces", handle3)
388 140 : IF (iw > 0) THEN
389 18 : WRITE (iw, '(T3,A)') " Evaluating DDAPC atomic derivatives .."
390 : END IF
391 700 : ALLOCATE (dAm(ndim, ndim, 3))
392 420 : ALLOCATE (dbv(ndim, 3))
393 700 : ALLOCATE (dqv(ndim, SIZE(particle_set), 3))
394 : !NB refactor math in inner loop - no more dqv0, but new temporaries instead
395 280 : ALLOCATE (cvT_AmI(ndim))
396 280 : ALLOCATE (cvT_AmI_dAmj(ndim))
397 420 : ALLOCATE (tv(ndim, SIZE(particle_set), 3))
398 280 : ALLOCATE (AmI_cv(ndim))
399 25772 : cvT_AmI(:) = MATMUL(cv, cp_ddapc_env%AmI)
400 25772 : AmI_cv(:) = MATMUL(cp_ddapc_env%AmI, cv)
401 280 : ALLOCATE (dAmj_qv(ndim))
402 280 : ALLOCATE (AmI_bv(ndim))
403 37508 : AmI_bv(:) = MATMUL(cp_ddapc_env%AmI, bv)
404 :
405 : !NB call routine to precompute sin(g.r) and cos(g.r),
406 : ! so it doesn't have to be done for each r_i-r_j pair in build_der_A_matrix_rows()
407 140 : CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
408 : !NB do build_der_A_matrix_rows in blocks, for more efficient use of DGEMM
409 : #define NPSET 100
410 280 : DO iparticle0 = 1, SIZE(particle_set), NPSET
411 140 : nparticles = MIN(NPSET, SIZE(particle_set) - iparticle0 + 1)
412 : !NB each dAm is supposed to have one block of rows and one block of columns
413 : !NB for derivatives with respect to each atom. build_der_A_matrix_rows()
414 : !NB just returns rows, since dAm is symmetric, and missing columns can be
415 : !NB reconstructed with a simple transpose, as below
416 : CALL build_der_A_matrix_rows(dAm, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
417 : particle_set, radii, rho_tot_g, gcut, iparticle0, &
418 140 : nparticles, g_dot_rvec_sin, g_dot_rvec_cos)
419 : !NB no more reduction of dbv and dAm - instead we go through with each node's contribution
420 : !NB and reduce resulting charges/forces once, at the end. Intermediate speedup can be
421 : !NB had by reducing dqv after the inner loop, and then other routines don't need to know
422 : !NB that contributions to dqv are distributed over the nodes.
423 : !NB also get rid of zeroing of dAm and division by Vol**2 - it's slow, and can be done
424 : !NB more quickly later, to a scalar or vector rather than a matrix
425 676 : DO iparticle = iparticle0, iparticle0 + nparticles - 1
426 : IF (debug_this_module) THEN
427 : CALL debug_der_A_matrix(dAm, particle_set, radii, rho_tot_g, &
428 : gcut, iparticle, Vol, qs_env)
429 : cp_ddapc_env => qs_env%cp_ddapc_env
430 : END IF
431 13572 : dbv(:, :) = 0.0_dp
432 : CALL build_der_b_vector(dbv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
433 396 : particle_set, radii, rho_tot_g, gcut, iparticle)
434 13572 : dbv(:, :) = dbv(:, :)/Vol
435 : IF (debug_this_module) THEN
436 : CALL debug_der_b_vector(dbv, particle_set, radii, rho_tot_g, &
437 : gcut, iparticle, Vol, qs_env)
438 : cp_ddapc_env => qs_env%cp_ddapc_env
439 : END IF
440 1584 : DO j = 1, 3
441 : !NB dAmj is actually pretty sparse - one block of cols + one block of rows - use this here:
442 1188 : pmin = (iparticle - 1)*SIZE(radii) + 1
443 1188 : pmax = iparticle*SIZE(radii)
444 : !NB multiply by block of columns that aren't explicitly in dAm, but can be reconstructured
445 : !NB as transpose of relevant block of rows
446 1188 : IF (pmin > 1) THEN
447 768 : dAmj_qv(:pmin - 1) = MATMUL(TRANSPOSE(dAm(pmin:pmax, :pmin - 1, j)), qv(pmin:pmax))
448 768 : cvT_AmI_dAmj(:pmin - 1) = MATMUL(TRANSPOSE(dAm(pmin:pmax, :pmin - 1, j)), cvT_AmI(pmin:pmax))
449 : END IF
450 : !NB multiply by block of rows that are explicitly in dAm
451 51624 : dAmj_qv(pmin:pmax) = MATMUL(dAm(pmin:pmax, :, j), qv(:))
452 51624 : cvT_AmI_dAmj(pmin:pmax) = MATMUL(dAm(pmin:pmax, :, j), cvT_AmI(:))
453 : !NB multiply by block of columns that aren't explicitly in dAm, but can be reconstructured
454 : !NB as transpose of relevant block of rows
455 1188 : IF (pmax < SIZE(particle_set)*SIZE(radii)) THEN
456 768 : dAmj_qv(pmax + 1:) = MATMUL(TRANSPOSE(dAm(pmin:pmax, pmax + 1:, j)), qv(pmin:pmax))
457 768 : cvT_AmI_dAmj(pmax + 1:) = MATMUL(TRANSPOSE(dAm(pmin:pmax, pmax + 1:, j)), cvT_AmI(pmin:pmax))
458 : END IF
459 13176 : dAmj_qv(:) = dAmj_qv(:)/(Vol*Vol)
460 13176 : cvT_AmI_dAmj(:) = cvT_AmI_dAmj(:)/(Vol*Vol)
461 37152 : c3 = DOT_PRODUCT(cvT_AmI_dAmj, AmI_bv) - DOT_PRODUCT(cvT_AmI, dbv(:, j)) - c1*DOT_PRODUCT(cvT_AmI_dAmj, AmI_cv)
462 13572 : tv(:, iparticle, j) = -(dAmj_qv(:) + dbv(:, j) + c3/cp_ddapc_env%c0*cv)
463 : END DO ! j
464 : !NB zero relevant parts of dAm here
465 48920 : dAm((iparticle - 1)*SIZE(radii) + 1:iparticle*SIZE(radii), :, :) = 0.0_dp
466 : !! dAm(:,(iparticle-1)*SIZE(radii)+1:iparticle*SIZE(radii),:) = 0.0_dp
467 : END DO ! iparticle
468 : END DO ! iparticle0
469 : !NB final part of refactoring of math - one dgemm is faster than many dgemv
470 : CALL dgemm('N', 'N', SIZE(dqv, 1), SIZE(dqv, 2)*SIZE(dqv, 3), SIZE(cp_ddapc_env%AmI, 2), 1.0_dp, &
471 140 : cp_ddapc_env%AmI, SIZE(cp_ddapc_env%AmI, 1), tv, SIZE(tv, 1), 0.0_dp, dqv, SIZE(dqv, 1))
472 : !NB deallocate g_dot_rvec_sin and g_dot_rvec_cos
473 140 : CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
474 : !NB moved reduction out to where dqv is used to compute
475 : !NB a force contribution (smaller array to reduce, just size(particle_set) x 3)
476 : !NB namely ewald_ddapc_force(), cp_decl_ddapc_forces(), restraint_functional_force()
477 140 : CPASSERT(PRESENT(dq_out))
478 140 : IF (.NOT. ASSOCIATED(dq_out)) THEN
479 700 : ALLOCATE (dq_out(SIZE(dqv, 1), SIZE(dqv, 2), SIZE(dqv, 3)))
480 : ELSE
481 0 : CPASSERT(SIZE(dqv, 1) == SIZE(dq_out, 1))
482 0 : CPASSERT(SIZE(dqv, 2) == SIZE(dq_out, 2))
483 0 : CPASSERT(SIZE(dqv, 3) == SIZE(dq_out, 3))
484 : END IF
485 13736 : dq_out = dqv
486 : IF (debug_this_module) THEN
487 : CALL debug_charge(dqv, qs_env, density_fit_section, &
488 : particle_set, radii, rho_tot_g, type_of_density)
489 : cp_ddapc_env => qs_env%cp_ddapc_env
490 : END IF
491 140 : DEALLOCATE (dqv)
492 140 : DEALLOCATE (dAm)
493 140 : DEALLOCATE (dbv)
494 : !NB deallocate new temporaries
495 140 : DEALLOCATE (cvT_AmI)
496 140 : DEALLOCATE (cvT_AmI_dAmj)
497 140 : DEALLOCATE (AmI_cv)
498 140 : DEALLOCATE (tv)
499 140 : DEALLOCATE (dAmj_qv)
500 140 : DEALLOCATE (AmI_bv)
501 280 : CALL timestop(handle3)
502 : END IF
503 : !
504 : ! End of charge fit
505 : !
506 1342 : DEALLOCATE (radii)
507 1342 : DEALLOCATE (bv)
508 1342 : DEALLOCATE (cv)
509 1342 : DEALLOCATE (qv)
510 1342 : DEALLOCATE (qtot)
511 1342 : IF (.NOT. PRESENT(iwc)) THEN
512 : CALL cp_print_key_finished_output(iw, logger, density_fit_section, &
513 740 : "PROGRAM_RUN_INFO")
514 : END IF
515 1342 : CALL auxbas_pool%give_back_pw(rho_tot_g)
516 1342 : CALL timestop(handle)
517 6710 : END SUBROUTINE get_ddapc
518 :
519 : ! **************************************************************************************************
520 : !> \brief modify hartree potential to handle restraints in DDAPC scheme
521 : !> \param v_hartree_gspace ...
522 : !> \param density_fit_section ...
523 : !> \param particle_set ...
524 : !> \param AmI ...
525 : !> \param radii ...
526 : !> \param charges ...
527 : !> \param ddapc_restraint_control ...
528 : !> \param energy_res ...
529 : !> \par History
530 : !> 02.2006 modified [Teo]
531 : ! **************************************************************************************************
532 500 : SUBROUTINE restraint_functional_potential(v_hartree_gspace, &
533 : density_fit_section, particle_set, AmI, radii, charges, &
534 : ddapc_restraint_control, energy_res)
535 : TYPE(pw_c1d_gs_type), INTENT(IN) :: v_hartree_gspace
536 : TYPE(section_vals_type), POINTER :: density_fit_section
537 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
538 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: AmI
539 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii, charges
540 : TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control
541 : REAL(KIND=dp), INTENT(INOUT) :: energy_res
542 :
543 : CHARACTER(len=*), PARAMETER :: routineN = 'restraint_functional_potential'
544 :
545 : COMPLEX(KIND=dp) :: g_corr, phase
546 : INTEGER :: handle, idim, ig, igauss, iparticle, &
547 : n_gauss
548 : REAL(KIND=dp) :: arg, fac, fac2, g2, gcut, gcut2, gfunc, &
549 : gvec(3), rc, rc2, rvec(3), sfac, Vol, w
550 500 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
551 :
552 500 : CALL timeset(routineN, handle)
553 500 : n_gauss = SIZE(radii)
554 1500 : ALLOCATE (cv(n_gauss*SIZE(particle_set)))
555 1000 : ALLOCATE (uv(n_gauss*SIZE(particle_set)))
556 2702 : uv = 0.0_dp
557 : CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
558 500 : charges, energy_res)
559 : !
560 500 : CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
561 500 : gcut2 = gcut*gcut
562 : ASSOCIATE (pw_grid => v_hartree_gspace%pw_grid)
563 500 : Vol = pw_grid%vol
564 2702 : cv = 1.0_dp/Vol
565 500 : sfac = -1.0_dp/Vol
566 36078 : fac = DOT_PRODUCT(cv, MATMUL(AmI, cv))
567 50064 : fac2 = DOT_PRODUCT(cv, MATMUL(AmI, uv))
568 2702 : cv(:) = uv - cv*fac2/fac
569 35078 : cv(:) = MATMUL(AmI, cv)
570 500 : IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
571 89314 : DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
572 88814 : g2 = pw_grid%gsq(ig)
573 88814 : w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
574 88814 : IF (g2 > gcut2) EXIT
575 353256 : gvec = pw_grid%g(:, ig)
576 88314 : g_corr = 0.0_dp
577 88314 : idim = 0
578 315696 : DO iparticle = 1, SIZE(particle_set)
579 879438 : DO igauss = 1, SIZE(radii)
580 563742 : idim = idim + 1
581 563742 : rc = radii(igauss)
582 563742 : rc2 = rc*rc
583 2254968 : rvec = particle_set(iparticle)%r
584 2254968 : arg = DOT_PRODUCT(gvec, rvec)
585 563742 : phase = CMPLX(COS(arg), -SIN(arg), KIND=dp)
586 563742 : gfunc = EXP(-g2*rc2/4.0_dp)
587 791124 : g_corr = g_corr + gfunc*cv(idim)*phase
588 : END DO
589 : END DO
590 88314 : g_corr = g_corr*w
591 88814 : v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/Vol
592 : END DO
593 : END ASSOCIATE
594 500 : CALL timestop(handle)
595 1500 : END SUBROUTINE restraint_functional_potential
596 :
597 : ! **************************************************************************************************
598 : !> \brief Modify the Hartree potential
599 : !> \param v_hartree_gspace ...
600 : !> \param density_fit_section ...
601 : !> \param particle_set ...
602 : !> \param M ...
603 : !> \param AmI ...
604 : !> \param radii ...
605 : !> \param charges ...
606 : !> \par History
607 : !> 08.2005 created [tlaino]
608 : !> \author Teodoro Laino
609 : ! **************************************************************************************************
610 740 : SUBROUTINE modify_hartree_pot(v_hartree_gspace, density_fit_section, &
611 : particle_set, M, AmI, radii, charges)
612 : TYPE(pw_c1d_gs_type), INTENT(IN) :: v_hartree_gspace
613 : TYPE(section_vals_type), POINTER :: density_fit_section
614 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
615 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M, AmI
616 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii, charges
617 :
618 : CHARACTER(len=*), PARAMETER :: routineN = 'modify_hartree_pot'
619 :
620 : COMPLEX(KIND=dp) :: g_corr, phase
621 : INTEGER :: handle, idim, ig, igauss, iparticle
622 : REAL(kind=dp) :: arg, fac, fac2, g2, gcut, gcut2, gfunc, &
623 : gvec(3), rc, rc2, rvec(3), sfac, Vol, w
624 740 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
625 :
626 740 : CALL timeset(routineN, handle)
627 740 : CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
628 740 : gcut2 = gcut*gcut
629 : ASSOCIATE (pw_grid => v_hartree_gspace%pw_grid)
630 740 : Vol = pw_grid%vol
631 2220 : ALLOCATE (cv(SIZE(M, 1)))
632 1480 : ALLOCATE (uv(SIZE(M, 1)))
633 7562 : cv = 1.0_dp/Vol
634 435098 : uv(:) = MATMUL(M, charges)
635 740 : sfac = -1.0_dp/Vol
636 303162 : fac = DOT_PRODUCT(cv, MATMUL(AmI, cv))
637 303162 : fac2 = DOT_PRODUCT(cv, MATMUL(AmI, uv))
638 7562 : cv(:) = uv - cv*fac2/fac
639 301682 : cv(:) = MATMUL(AmI, cv)
640 740 : IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
641 695612 : DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
642 694872 : g2 = pw_grid%gsq(ig)
643 694872 : w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
644 694872 : IF (g2 > gcut2) EXIT
645 2776528 : gvec = pw_grid%g(:, ig)
646 694132 : g_corr = 0.0_dp
647 694132 : idim = 0
648 2972942 : DO iparticle = 1, SIZE(particle_set)
649 9809372 : DO igauss = 1, SIZE(radii)
650 6836430 : idim = idim + 1
651 6836430 : rc = radii(igauss)
652 6836430 : rc2 = rc*rc
653 27345720 : rvec = particle_set(iparticle)%r
654 27345720 : arg = DOT_PRODUCT(gvec, rvec)
655 6836430 : phase = CMPLX(COS(arg), -SIN(arg), KIND=dp)
656 6836430 : gfunc = EXP(-g2*rc2/4.0_dp)
657 9115240 : g_corr = g_corr + gfunc*cv(idim)*phase
658 : END DO
659 : END DO
660 694132 : g_corr = g_corr*w
661 694872 : v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/Vol
662 : END DO
663 : END ASSOCIATE
664 740 : CALL timestop(handle)
665 1480 : END SUBROUTINE modify_hartree_pot
666 :
667 : ! **************************************************************************************************
668 : !> \brief To Debug the derivative of the B vector for the solution of the
669 : !> linear system
670 : !> \param dbv ...
671 : !> \param particle_set ...
672 : !> \param radii ...
673 : !> \param rho_tot_g ...
674 : !> \param gcut ...
675 : !> \param iparticle ...
676 : !> \param Vol ...
677 : !> \param qs_env ...
678 : !> \par History
679 : !> 08.2005 created [tlaino]
680 : !> \author Teodoro Laino
681 : ! **************************************************************************************************
682 0 : SUBROUTINE debug_der_b_vector(dbv, particle_set, radii, &
683 : rho_tot_g, gcut, iparticle, Vol, qs_env)
684 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: dbv
685 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
686 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
687 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
688 : REAL(KIND=dp), INTENT(IN) :: gcut
689 : INTEGER, INTENT(in) :: iparticle
690 : REAL(KIND=dp), INTENT(IN) :: Vol
691 : TYPE(qs_environment_type), POINTER :: qs_env
692 :
693 : CHARACTER(len=*), PARAMETER :: routineN = 'debug_der_b_vector'
694 :
695 : INTEGER :: handle, i, kk, ndim
696 : REAL(KIND=dp) :: dx, rvec(3), v0
697 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bv1, bv2, ddbv
698 : TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env
699 :
700 0 : NULLIFY (cp_ddapc_env)
701 0 : CALL timeset(routineN, handle)
702 0 : dx = 0.01_dp
703 0 : ndim = SIZE(particle_set)*SIZE(radii)
704 0 : ALLOCATE (bv1(ndim))
705 0 : ALLOCATE (bv2(ndim))
706 0 : ALLOCATE (ddbv(ndim))
707 0 : rvec = particle_set(iparticle)%r
708 0 : cp_ddapc_env => qs_env%cp_ddapc_env
709 0 : DO i = 1, 3
710 0 : bv1(:) = 0.0_dp
711 0 : bv2(:) = 0.0_dp
712 0 : particle_set(iparticle)%r(i) = rvec(i) + dx
713 : CALL build_b_vector(bv1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
714 0 : particle_set, radii, rho_tot_g, gcut); bv1(:) = bv1(:)/Vol
715 0 : CALL rho_tot_g%pw_grid%para%group%sum(bv1)
716 0 : particle_set(iparticle)%r(i) = rvec(i) - dx
717 : CALL build_b_vector(bv2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
718 0 : particle_set, radii, rho_tot_g, gcut); bv2(:) = bv2(:)/Vol
719 0 : CALL rho_tot_g%pw_grid%para%group%sum(bv2)
720 0 : ddbv(:) = (bv1(:) - bv2(:))/(2.0_dp*dx)
721 0 : DO kk = 1, SIZE(ddbv)
722 0 : IF (ddbv(kk) .GT. 1.0E-8_dp) THEN
723 0 : v0 = ABS(dbv(kk, i) - ddbv(kk))/ddbv(kk)*100.0_dp
724 0 : WRITE (*, *) "Error % on B ::", v0
725 0 : IF (v0 .GT. 0.1_dp) THEN
726 0 : WRITE (*, '(A,2I5,2F15.9)') "ERROR IN DERIVATIVE OF B VECTOR, IPARTICLE, ICOORD:", iparticle, i, &
727 0 : dbv(kk, i), ddbv(kk)
728 0 : CPABORT("")
729 : END IF
730 : END IF
731 : END DO
732 0 : particle_set(iparticle)%r = rvec
733 : END DO
734 0 : DEALLOCATE (bv1)
735 0 : DEALLOCATE (bv2)
736 0 : DEALLOCATE (ddbv)
737 0 : CALL timestop(handle)
738 0 : END SUBROUTINE debug_der_b_vector
739 :
740 : ! **************************************************************************************************
741 : !> \brief To Debug the derivative of the A matrix for the solution of the
742 : !> linear system
743 : !> \param dAm ...
744 : !> \param particle_set ...
745 : !> \param radii ...
746 : !> \param rho_tot_g ...
747 : !> \param gcut ...
748 : !> \param iparticle ...
749 : !> \param Vol ...
750 : !> \param qs_env ...
751 : !> \par History
752 : !> 08.2005 created [tlaino]
753 : !> \author Teodoro Laino
754 : ! **************************************************************************************************
755 0 : SUBROUTINE debug_der_A_matrix(dAm, particle_set, radii, &
756 : rho_tot_g, gcut, iparticle, Vol, qs_env)
757 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dAm
758 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
759 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
760 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
761 : REAL(KIND=dp), INTENT(IN) :: gcut
762 : INTEGER, INTENT(in) :: iparticle
763 : REAL(KIND=dp), INTENT(IN) :: Vol
764 : TYPE(qs_environment_type), POINTER :: qs_env
765 :
766 : CHARACTER(len=*), PARAMETER :: routineN = 'debug_der_A_matrix'
767 :
768 : INTEGER :: handle, i, kk, ll, ndim
769 : REAL(KIND=dp) :: dx, rvec(3), v0
770 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Am1, Am2, ddAm, g_dot_rvec_cos, &
771 0 : g_dot_rvec_sin
772 : TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env
773 :
774 : !NB new temporaries sin(g.r) and cos(g.r), as used in get_ddapc, to speed up build_der_A_matrix()
775 :
776 0 : NULLIFY (cp_ddapc_env)
777 0 : CALL timeset(routineN, handle)
778 0 : dx = 0.01_dp
779 0 : ndim = SIZE(particle_set)*SIZE(radii)
780 0 : ALLOCATE (Am1(ndim, ndim))
781 0 : ALLOCATE (Am2(ndim, ndim))
782 0 : ALLOCATE (ddAm(ndim, ndim))
783 0 : rvec = particle_set(iparticle)%r
784 0 : cp_ddapc_env => qs_env%cp_ddapc_env
785 0 : CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
786 0 : DO i = 1, 3
787 0 : Am1 = 0.0_dp
788 0 : Am2 = 0.0_dp
789 0 : particle_set(iparticle)%r(i) = rvec(i) + dx
790 : CALL build_A_matrix(Am1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
791 0 : particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
792 0 : Am1(:, :) = Am1(:, :)/(Vol*Vol)
793 0 : CALL rho_tot_g%pw_grid%para%group%sum(Am1)
794 0 : particle_set(iparticle)%r(i) = rvec(i) - dx
795 : CALL build_A_matrix(Am2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
796 0 : particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
797 0 : Am2(:, :) = Am2(:, :)/(Vol*Vol)
798 0 : CALL rho_tot_g%pw_grid%para%group%sum(Am2)
799 0 : ddAm(:, :) = (Am1 - Am2)/(2.0_dp*dx)
800 0 : DO kk = 1, SIZE(ddAm, 1)
801 0 : DO ll = 1, SIZE(ddAm, 2)
802 0 : IF (ddAm(kk, ll) .GT. 1.0E-8_dp) THEN
803 0 : v0 = ABS(dAm(kk, ll, i) - ddAm(kk, ll))/ddAm(kk, ll)*100.0_dp
804 0 : WRITE (*, *) "Error % on A ::", v0, Am1(kk, ll), Am2(kk, ll), iparticle, i, kk, ll
805 0 : IF (v0 .GT. 0.1_dp) THEN
806 0 : WRITE (*, '(A,4I5,2F15.9)') "ERROR IN DERIVATIVE OF A MATRIX, IPARTICLE, ICOORD:", iparticle, i, kk, ll, &
807 0 : dAm(kk, ll, i), ddAm(kk, ll)
808 0 : CPABORT("")
809 : END IF
810 : END IF
811 : END DO
812 : END DO
813 0 : particle_set(iparticle)%r = rvec
814 : END DO
815 0 : CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
816 0 : DEALLOCATE (Am1)
817 0 : DEALLOCATE (Am2)
818 0 : DEALLOCATE (ddAm)
819 0 : CALL timestop(handle)
820 0 : END SUBROUTINE debug_der_A_matrix
821 :
822 : ! **************************************************************************************************
823 : !> \brief To Debug the fitted charges
824 : !> \param dqv ...
825 : !> \param qs_env ...
826 : !> \param density_fit_section ...
827 : !> \param particle_set ...
828 : !> \param radii ...
829 : !> \param rho_tot_g ...
830 : !> \param type_of_density ...
831 : !> \par History
832 : !> 08.2005 created [tlaino]
833 : !> \author Teodoro Laino
834 : ! **************************************************************************************************
835 0 : SUBROUTINE debug_charge(dqv, qs_env, density_fit_section, &
836 : particle_set, radii, rho_tot_g, type_of_density)
837 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dqv
838 : TYPE(qs_environment_type), POINTER :: qs_env
839 : TYPE(section_vals_type), POINTER :: density_fit_section
840 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
841 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
842 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
843 : CHARACTER(LEN=*) :: type_of_density
844 :
845 : CHARACTER(len=*), PARAMETER :: routineN = 'debug_charge'
846 :
847 : INTEGER :: handle, i, iparticle, kk, ndim
848 : REAL(KIND=dp) :: dx, rvec(3)
849 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ddqv
850 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: qtot1, qtot2
851 :
852 0 : CALL timeset(routineN, handle)
853 0 : WRITE (*, *) "DEBUG_CHARGE_ROUTINE"
854 0 : ndim = SIZE(particle_set)*SIZE(radii)
855 0 : NULLIFY (qtot1, qtot2)
856 0 : ALLOCATE (qtot1(ndim))
857 0 : ALLOCATE (qtot2(ndim))
858 0 : ALLOCATE (ddqv(ndim))
859 : !
860 : dx = 0.001_dp
861 0 : DO iparticle = 1, SIZE(particle_set)
862 0 : rvec = particle_set(iparticle)%r
863 0 : DO i = 1, 3
864 0 : particle_set(iparticle)%r(i) = rvec(i) + dx
865 : CALL get_ddapc(qs_env, .FALSE., density_fit_section, qout1=qtot1, &
866 0 : ext_rho_tot_g=rho_tot_g, Itype_of_density=type_of_density)
867 0 : particle_set(iparticle)%r(i) = rvec(i) - dx
868 : CALL get_ddapc(qs_env, .FALSE., density_fit_section, qout1=qtot2, &
869 0 : ext_rho_tot_g=rho_tot_g, Itype_of_density=type_of_density)
870 0 : ddqv(:) = (qtot1 - qtot2)/(2.0_dp*dx)
871 0 : DO kk = 1, SIZE(qtot1) - 1, SIZE(radii)
872 0 : IF (ANY(ddqv(kk:kk + 2) .GT. 1.0E-8_dp)) THEN
873 0 : WRITE (*, '(A,2F12.6,F12.2)') "Error :", SUM(dqv(kk:kk + 2, iparticle, i)), SUM(ddqv(kk:kk + 2)), &
874 0 : ABS((SUM(ddqv(kk:kk + 2)) - SUM(dqv(kk:kk + 2, iparticle, i)))/SUM(ddqv(kk:kk + 2))*100.0_dp)
875 : END IF
876 : END DO
877 0 : particle_set(iparticle)%r = rvec
878 : END DO
879 : END DO
880 : !
881 0 : DEALLOCATE (qtot1)
882 0 : DEALLOCATE (qtot2)
883 0 : DEALLOCATE (ddqv)
884 0 : CALL timestop(handle)
885 0 : END SUBROUTINE debug_charge
886 :
887 6404 : END MODULE cp_ddapc_util
|
