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