(git:b17b328)
Loading...
Searching...
No Matches
cp_ddapc_util.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
70CONTAINS
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 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 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 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 CALL timeset(routinen, handle)
101 logger => cp_get_default_logger()
102 NULLIFY (dft_control, rho, pw_env, &
103 radii, inp_radii, particle_set, qs_charges, para_env)
104
105 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 qs_env%cp_ddapc_ewald%do_restraint
110 unimplemented = dft_control%qs_control%semi_empirical .OR. &
111 dft_control%qs_control%dftb .OR. &
112 dft_control%qs_control%xtb
113 IF (allocate_ddapc_env .AND. unimplemented) THEN
114 cpabort("DDAP charges work only with GPW/GAPW code.")
115 END IF
116 allocate_ddapc_env = allocate_ddapc_env .OR. &
117 qs_env%cp_ddapc_ewald%do_property
118 allocate_ddapc_env = allocate_ddapc_env .AND. (.NOT. unimplemented)
119 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 para_env=para_env)
129 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 "PROGRAM_RUN_INFO", ".FitCharge")
132 IF (iw > 0) THEN
133 WRITE (iw, '(/,A)') " Initializing the DDAPC Environment"
134 END IF
135 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pool)
136 CALL auxbas_pool%create_pw(rho_tot_g)
137 vol = rho_tot_g%pw_grid%vol
138 !
139 ! Get Input Parameters
140 !
141 CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val)
142 IF (n_rep_val /= 0) THEN
143 CALL section_vals_val_get(density_fit_section, "RADII", r_vals=inp_radii)
144 num_gauss = SIZE(inp_radii)
145 ALLOCATE (radii(num_gauss))
146 radii = inp_radii
147 ELSE
148 CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss)
149 CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin)
150 CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact)
151 ALLOCATE (radii(num_gauss))
152 DO i = 1, num_gauss
153 radii(i) = rcmin*pfact**(i - 1)
154 END DO
155 END IF
156 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 "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 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 qs_env%input)
175 CALL cp_print_key_finished_output(iw2, logger, density_fit_section, &
176 "PROGRAM_RUN_INFO/CONDITION_NUMBER")
177 DEALLOCATE (radii)
178 CALL auxbas_pool%give_back_pw(rho_tot_g)
179 END IF
180 CALL timestop(handle)
181 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ami_bv, ami_cv, bv, cv, cvt_ami, &
224 cvt_ami_damj, damj_qv, qtot, qv
225 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dbv, g_dot_rvec_cos, g_dot_rvec_sin
226 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dam, dqv, tv
227 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
233 TYPE(pw_c1d_gs_type) :: rho_tot_g
234 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 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
239 TYPE(qs_charges_type), POINTER :: qs_charges
240 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 CALL timeset(routinen, handle)
250 need_f = .false.
251 IF (PRESENT(calc_force)) need_f = calc_force
252 logger => cp_get_default_logger()
253 NULLIFY (dft_control, rho, rho_core, rho0_s_gs, rhoz_cneo_s_gs, pw_env, rho_g, rho_r, &
254 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 super_cell=super_cell)
267
268 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
269
270 IF (PRESENT(iwc)) THEN
271 iw = iwc
272 ELSE
273 iw = cp_print_key_unit_nr(logger, density_fit_section, &
274 "PROGRAM_RUN_INFO", ".FitCharge")
275 END IF
276 CALL pw_env_get(pw_env=pw_env, &
277 auxbas_pw_pool=auxbas_pool)
278 CALL auxbas_pool%create_pw(rho_tot_g)
279 IF (PRESENT(ext_rho_tot_g)) THEN
280 ! If provided use the input density in g-space
281 CALL pw_transfer(ext_rho_tot_g, rho_tot_g)
282 type_of_density = itype_of_density
283 ELSE
284 IF (PRESENT(density_type)) THEN
285 myid = density_type
286 ELSE
287 CALL section_vals_val_get(qs_env%input, &
288 "PROPERTIES%FIT_CHARGE%TYPE_OF_DENSITY", i_val=myid)
289 END IF
290 SELECT CASE (myid)
291 CASE (do_full_density)
292 ! Otherwise build the total QS density (electron+nuclei) in G-space
293 IF (dft_control%qs_control%gapw) THEN
294 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
295 CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
296 END IF
297 CALL pw_transfer(rho0_s_gs, rho_tot_g)
298 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
299 CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
300 END IF
301 ELSE
302 CALL pw_transfer(rho_core, rho_tot_g)
303 END IF
304 DO ispin = 1, SIZE(rho_g)
305 CALL pw_axpy(rho_g(ispin), rho_tot_g)
306 END DO
307 type_of_density = "FULL DENSITY"
308 CASE (do_spin_density)
309 CALL pw_copy(rho_g(1), rho_tot_g)
310 CALL pw_axpy(rho_g(2), rho_tot_g, alpha=-1._dp)
311 type_of_density = "SPIN DENSITY"
312 END SELECT
313 END IF
314 vol = rho_r(1)%pw_grid%vol
315 ch_dens = 0.0_dp
316 ! should use pw_integrate
317 IF (rho_tot_g%pw_grid%have_g0) ch_dens = real(rho_tot_g%array(1), kind=dp)
318 CALL logger%para_env%sum(ch_dens)
319 !
320 ! Get Input Parameters
321 !
322 CALL section_vals_val_get(density_fit_section, "RADII", n_rep_val=n_rep_val)
323 IF (n_rep_val /= 0) THEN
324 CALL section_vals_val_get(density_fit_section, "RADII", r_vals=inp_radii)
325 num_gauss = SIZE(inp_radii)
326 ALLOCATE (radii(num_gauss))
327 radii = inp_radii
328 ELSE
329 CALL section_vals_val_get(density_fit_section, "NUM_GAUSS", i_val=num_gauss)
330 CALL section_vals_val_get(density_fit_section, "MIN_RADIUS", r_val=rcmin)
331 CALL section_vals_val_get(density_fit_section, "PFACTOR", r_val=pfact)
332 ALLOCATE (radii(num_gauss))
333 DO i = 1, num_gauss
334 radii(i) = rcmin*pfact**(i - 1)
335 END DO
336 END IF
337 IF (PRESENT(out_radii)) THEN
338 IF (ASSOCIATED(out_radii)) THEN
339 DEALLOCATE (out_radii)
340 END IF
341 ALLOCATE (out_radii(SIZE(radii)))
342 out_radii = radii
343 END IF
344 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
345 cp_ddapc_env => qs_env%cp_ddapc_env
346 !
347 ! Start with the linear system
348 !
349 ndim = SIZE(particle_set)*SIZE(radii)
350 ALLOCATE (bv(ndim))
351 ALLOCATE (qv(ndim))
352 ALLOCATE (qtot(SIZE(particle_set)))
353 ALLOCATE (cv(ndim))
354 CALL timeset(routinen//"-charges", handle2)
355 bv(:) = 0.0_dp
356 cv(:) = 1.0_dp/vol
357 CALL build_b_vector(bv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
358 particle_set, radii, rho_tot_g, gcut); bv(:) = bv(:)/vol
359 CALL rho_tot_g%pw_grid%para%group%sum(bv)
360 c1 = dot_product(cv, matmul(cp_ddapc_env%AmI, bv)) - ch_dens
361 c1 = c1/cp_ddapc_env%c0
362 qv(:) = -matmul(cp_ddapc_env%AmI, (bv - c1*cv))
363 j = 0
364 qtot = 0.0_dp
365 DO i = 1, ndim, num_gauss
366 j = j + 1
367 DO ii = 1, num_gauss
368 qtot(j) = qtot(j) + qv((i - 1) + ii)
369 END DO
370 END DO
371 IF (PRESENT(qout1)) THEN
372 IF (ASSOCIATED(qout1)) THEN
373 cpassert(SIZE(qout1) == SIZE(qv))
374 ELSE
375 ALLOCATE (qout1(SIZE(qv)))
376 END IF
377 qout1 = qv
378 END IF
379 IF (PRESENT(qout2)) THEN
380 IF (ASSOCIATED(qout2)) THEN
381 cpassert(SIZE(qout2) == SIZE(qtot))
382 ELSE
383 ALLOCATE (qout2(SIZE(qtot)))
384 END IF
385 qout2 = qtot
386 END IF
387 CALL print_atomic_charges(particle_set, qs_kind_set, iw, title=" DDAP "// &
388 trim(type_of_density)//" charges:", atomic_charges=qtot)
389 CALL timestop(handle2)
390 !
391 ! If requested evaluate also the correction to derivatives due to Pulay Forces
392 !
393 IF (need_f) THEN
394 CALL timeset(routinen//"-forces", handle3)
395 IF (iw > 0) THEN
396 WRITE (iw, '(T3,A)') " Evaluating DDAPC atomic derivatives .."
397 END IF
398 ALLOCATE (dam(ndim, ndim, 3))
399 ALLOCATE (dbv(ndim, 3))
400 ALLOCATE (dqv(ndim, SIZE(particle_set), 3))
401 !NB refactor math in inner loop - no more dqv0, but new temporaries instead
402 ALLOCATE (cvt_ami(ndim))
403 ALLOCATE (cvt_ami_damj(ndim))
404 ALLOCATE (tv(ndim, SIZE(particle_set), 3))
405 ALLOCATE (ami_cv(ndim))
406 cvt_ami(:) = matmul(cv, cp_ddapc_env%AmI)
407 ami_cv(:) = matmul(cp_ddapc_env%AmI, cv)
408 ALLOCATE (damj_qv(ndim))
409 ALLOCATE (ami_bv(ndim))
410 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 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 DO iparticle0 = 1, SIZE(particle_set), npset
418 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 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 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 dbv(:, :) = 0.0_dp
439 CALL build_der_b_vector(dbv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
440 particle_set, radii, rho_tot_g, gcut, iparticle)
441 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 DO j = 1, 3
448 !NB dAmj is actually pretty sparse - one block of cols + one block of rows - use this here:
449 pmin = (iparticle - 1)*SIZE(radii) + 1
450 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 IF (pmin > 1) THEN
454 damj_qv(:pmin - 1) = matmul(transpose(dam(pmin:pmax, :pmin - 1, j)), qv(pmin:pmax))
455 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 damj_qv(pmin:pmax) = matmul(dam(pmin:pmax, :, j), qv(:))
459 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 IF (pmax < SIZE(particle_set)*SIZE(radii)) THEN
463 damj_qv(pmax + 1:) = matmul(transpose(dam(pmin:pmax, pmax + 1:, j)), qv(pmin:pmax))
464 cvt_ami_damj(pmax + 1:) = matmul(transpose(dam(pmin:pmax, pmax + 1:, j)), cvt_ami(pmin:pmax))
465 END IF
466 damj_qv(:) = damj_qv(:)/(vol*vol)
467 cvt_ami_damj(:) = cvt_ami_damj(:)/(vol*vol)
468 c3 = dot_product(cvt_ami_damj, ami_bv) - dot_product(cvt_ami, dbv(:, j)) - c1*dot_product(cvt_ami_damj, ami_cv)
469 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 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 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 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 cpassert(PRESENT(dq_out))
485 IF (.NOT. ASSOCIATED(dq_out)) THEN
486 ALLOCATE (dq_out(SIZE(dqv, 1), SIZE(dqv, 2), SIZE(dqv, 3)))
487 ELSE
488 cpassert(SIZE(dqv, 1) == SIZE(dq_out, 1))
489 cpassert(SIZE(dqv, 2) == SIZE(dq_out, 2))
490 cpassert(SIZE(dqv, 3) == SIZE(dq_out, 3))
491 END IF
492 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 DEALLOCATE (dqv)
499 DEALLOCATE (dam)
500 DEALLOCATE (dbv)
501 !NB deallocate new temporaries
502 DEALLOCATE (cvt_ami)
503 DEALLOCATE (cvt_ami_damj)
504 DEALLOCATE (ami_cv)
505 DEALLOCATE (tv)
506 DEALLOCATE (damj_qv)
507 DEALLOCATE (ami_bv)
508 CALL timestop(handle3)
509 END IF
510 !
511 ! End of charge fit
512 !
513 DEALLOCATE (radii)
514 DEALLOCATE (bv)
515 DEALLOCATE (cv)
516 DEALLOCATE (qv)
517 DEALLOCATE (qtot)
518 IF (.NOT. PRESENT(iwc)) THEN
519 CALL cp_print_key_finished_output(iw, logger, density_fit_section, &
520 "PROGRAM_RUN_INFO")
521 END IF
522 CALL auxbas_pool%give_back_pw(rho_tot_g)
523 CALL timestop(handle)
524 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
558
559 CALL timeset(routinen, handle)
560 n_gauss = SIZE(radii)
561 ALLOCATE (cv(n_gauss*SIZE(particle_set)))
562 ALLOCATE (uv(n_gauss*SIZE(particle_set)))
563 uv = 0.0_dp
564 CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
565 charges, energy_res)
566 !
567 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
568 gcut2 = gcut*gcut
569 associate(pw_grid => v_hartree_gspace%pw_grid)
570 vol = pw_grid%vol
571 cv = 1.0_dp/vol
572 sfac = -1.0_dp/vol
573 fac = dot_product(cv, matmul(ami, cv))
574 fac2 = dot_product(cv, matmul(ami, uv))
575 cv(:) = uv - cv*fac2/fac
576 cv(:) = matmul(ami, cv)
577 IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
578 DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
579 g2 = pw_grid%gsq(ig)
580 w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
581 IF (g2 > gcut2) EXIT
582 gvec = pw_grid%g(:, ig)
583 g_corr = 0.0_dp
584 idim = 0
585 DO iparticle = 1, SIZE(particle_set)
586 DO igauss = 1, SIZE(radii)
587 idim = idim + 1
588 rc = radii(igauss)
589 rc2 = rc*rc
590 rvec = particle_set(iparticle)%r
591 arg = dot_product(gvec, rvec)
592 phase = cmplx(cos(arg), -sin(arg), kind=dp)
593 gfunc = exp(-g2*rc2/4.0_dp)
594 g_corr = g_corr + gfunc*cv(idim)*phase
595 END DO
596 END DO
597 g_corr = g_corr*w
598 v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/vol
599 END DO
600 END associate
601 CALL timestop(handle)
602 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
632
633 CALL timeset(routinen, handle)
634 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
635 gcut2 = gcut*gcut
636 associate(pw_grid => v_hartree_gspace%pw_grid)
637 vol = pw_grid%vol
638 ALLOCATE (cv(SIZE(m, 1)))
639 ALLOCATE (uv(SIZE(m, 1)))
640 cv = 1.0_dp/vol
641 uv(:) = matmul(m, charges)
642 sfac = -1.0_dp/vol
643 fac = dot_product(cv, matmul(ami, cv))
644 fac2 = dot_product(cv, matmul(ami, uv))
645 cv(:) = uv - cv*fac2/fac
646 cv(:) = matmul(ami, cv)
647 IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
648 DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
649 g2 = pw_grid%gsq(ig)
650 w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
651 IF (g2 > gcut2) EXIT
652 gvec = pw_grid%g(:, ig)
653 g_corr = 0.0_dp
654 idim = 0
655 DO iparticle = 1, SIZE(particle_set)
656 DO igauss = 1, SIZE(radii)
657 idim = idim + 1
658 rc = radii(igauss)
659 rc2 = rc*rc
660 rvec = particle_set(iparticle)%r
661 arg = dot_product(gvec, rvec)
662 phase = cmplx(cos(arg), -sin(arg), kind=dp)
663 gfunc = exp(-g2*rc2/4.0_dp)
664 g_corr = g_corr + gfunc*cv(idim)*phase
665 END DO
666 END DO
667 g_corr = g_corr*w
668 v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/vol
669 END DO
670 END associate
671 CALL timestop(handle)
672 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: bv1, bv2, ddbv
705 TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env
706
707 NULLIFY (cp_ddapc_env)
708 CALL timeset(routinen, handle)
709 dx = 0.01_dp
710 ndim = SIZE(particle_set)*SIZE(radii)
711 ALLOCATE (bv1(ndim))
712 ALLOCATE (bv2(ndim))
713 ALLOCATE (ddbv(ndim))
714 rvec = particle_set(iparticle)%r
715 cp_ddapc_env => qs_env%cp_ddapc_env
716 DO i = 1, 3
717 bv1(:) = 0.0_dp
718 bv2(:) = 0.0_dp
719 particle_set(iparticle)%r(i) = rvec(i) + dx
720 CALL build_b_vector(bv1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
721 particle_set, radii, rho_tot_g, gcut); bv1(:) = bv1(:)/vol
722 CALL rho_tot_g%pw_grid%para%group%sum(bv1)
723 particle_set(iparticle)%r(i) = rvec(i) - dx
724 CALL build_b_vector(bv2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
725 particle_set, radii, rho_tot_g, gcut); bv2(:) = bv2(:)/vol
726 CALL rho_tot_g%pw_grid%para%group%sum(bv2)
727 ddbv(:) = (bv1(:) - bv2(:))/(2.0_dp*dx)
728 DO kk = 1, SIZE(ddbv)
729 IF (ddbv(kk) .GT. 1.0e-8_dp) THEN
730 v0 = abs(dbv(kk, i) - ddbv(kk))/ddbv(kk)*100.0_dp
731 WRITE (*, *) "Error % on B ::", v0
732 IF (v0 .GT. 0.1_dp) THEN
733 WRITE (*, '(A,2I5,2F15.9)') "ERROR IN DERIVATIVE OF B VECTOR, IPARTICLE, ICOORD:", iparticle, i, &
734 dbv(kk, i), ddbv(kk)
735 cpabort("")
736 END IF
737 END IF
738 END DO
739 particle_set(iparticle)%r = rvec
740 END DO
741 DEALLOCATE (bv1)
742 DEALLOCATE (bv2)
743 DEALLOCATE (ddbv)
744 CALL timestop(handle)
745 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: am1, am2, ddam, g_dot_rvec_cos, &
778 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 NULLIFY (cp_ddapc_env)
784 CALL timeset(routinen, handle)
785 dx = 0.01_dp
786 ndim = SIZE(particle_set)*SIZE(radii)
787 ALLOCATE (am1(ndim, ndim))
788 ALLOCATE (am2(ndim, ndim))
789 ALLOCATE (ddam(ndim, ndim))
790 rvec = particle_set(iparticle)%r
791 cp_ddapc_env => qs_env%cp_ddapc_env
792 CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
793 DO i = 1, 3
794 am1 = 0.0_dp
795 am2 = 0.0_dp
796 particle_set(iparticle)%r(i) = rvec(i) + dx
797 CALL build_a_matrix(am1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
798 particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
799 am1(:, :) = am1(:, :)/(vol*vol)
800 CALL rho_tot_g%pw_grid%para%group%sum(am1)
801 particle_set(iparticle)%r(i) = rvec(i) - dx
802 CALL build_a_matrix(am2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
803 particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
804 am2(:, :) = am2(:, :)/(vol*vol)
805 CALL rho_tot_g%pw_grid%para%group%sum(am2)
806 ddam(:, :) = (am1 - am2)/(2.0_dp*dx)
807 DO kk = 1, SIZE(ddam, 1)
808 DO ll = 1, SIZE(ddam, 2)
809 IF (ddam(kk, ll) .GT. 1.0e-8_dp) THEN
810 v0 = abs(dam(kk, ll, i) - ddam(kk, ll))/ddam(kk, ll)*100.0_dp
811 WRITE (*, *) "Error % on A ::", v0, am1(kk, ll), am2(kk, ll), iparticle, i, kk, ll
812 IF (v0 .GT. 0.1_dp) THEN
813 WRITE (*, '(A,4I5,2F15.9)') "ERROR IN DERIVATIVE OF A MATRIX, IPARTICLE, ICOORD:", iparticle, i, kk, ll, &
814 dam(kk, ll, i), ddam(kk, ll)
815 cpabort("")
816 END IF
817 END IF
818 END DO
819 END DO
820 particle_set(iparticle)%r = rvec
821 END DO
822 CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
823 DEALLOCATE (am1)
824 DEALLOCATE (am2)
825 DEALLOCATE (ddam)
826 CALL timestop(handle)
827 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ddqv
857 REAL(kind=dp), DIMENSION(:), POINTER :: qtot1, qtot2
858
859 CALL timeset(routinen, handle)
860 WRITE (*, *) "DEBUG_CHARGE_ROUTINE"
861 ndim = SIZE(particle_set)*SIZE(radii)
862 NULLIFY (qtot1, qtot2)
863 ALLOCATE (qtot1(ndim))
864 ALLOCATE (qtot2(ndim))
865 ALLOCATE (ddqv(ndim))
866 !
867 dx = 0.001_dp
868 DO iparticle = 1, SIZE(particle_set)
869 rvec = particle_set(iparticle)%r
870 DO i = 1, 3
871 particle_set(iparticle)%r(i) = rvec(i) + dx
872 CALL get_ddapc(qs_env, .false., density_fit_section, qout1=qtot1, &
873 ext_rho_tot_g=rho_tot_g, itype_of_density=type_of_density)
874 particle_set(iparticle)%r(i) = rvec(i) - dx
875 CALL get_ddapc(qs_env, .false., density_fit_section, qout1=qtot2, &
876 ext_rho_tot_g=rho_tot_g, itype_of_density=type_of_density)
877 ddqv(:) = (qtot1 - qtot2)/(2.0_dp*dx)
878 DO kk = 1, SIZE(qtot1) - 1, SIZE(radii)
879 IF (any(ddqv(kk:kk + 2) .GT. 1.0e-8_dp)) THEN
880 WRITE (*, '(A,2F12.6,F12.2)') "Error :", sum(dqv(kk:kk + 2, iparticle, i)), sum(ddqv(kk:kk + 2)), &
881 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 particle_set(iparticle)%r = rvec
885 END DO
886 END DO
887 !
888 DEALLOCATE (qtot1)
889 DEALLOCATE (qtot2)
890 DEALLOCATE (ddqv)
891 CALL timestop(handle)
892 END SUBROUTINE debug_charge
893
894END MODULE cp_ddapc_util
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:51
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:214
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_copy
Definition pw_methods.F:145
pw_methods::pw_transfer
Definition pw_methods.F:303
atomic_charges
simple routine to print charges for all atomic charge methods (currently mulliken,...
Definition atomic_charges.F:14
atomic_charges::print_atomic_charges
subroutine, public print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, atomic_charges)
generates a unified output format for atomic charges
Definition atomic_charges.F:48
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_ddapc_forces
Density Derived atomic point charges from a QM calculation (see J. Chem. Phys. Vol....
Definition cp_ddapc_forces.F:15
cp_ddapc_forces::evaluate_restraint_functional
subroutine, public evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, charges, energy_res)
computes energy and derivatives given a set of charges
Definition cp_ddapc_forces.F:330
cp_ddapc_methods
contains information regarding the decoupling/recoupling method of Bloechl
Definition cp_ddapc_methods.F:12
cp_ddapc_methods::cleanup_g_dot_rvec_sin_cos
subroutine, public cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
deallocate g_dot_rvec_* arrays
Definition cp_ddapc_methods.F:477
cp_ddapc_methods::build_der_b_vector
subroutine, public build_der_b_vector(dbv, gfunc, w, particle_set, radii, rho_tot_g, gcut, iparticle0)
Computes the derivative of B vector for the evaluation of the Pulay forces.
Definition cp_ddapc_methods.F:277
cp_ddapc_methods::build_b_vector
subroutine, public build_b_vector(bv, gfunc, w, particle_set, radii, rho_tot_g, gcut)
Computes the B vector for the solution of the linear system.
Definition cp_ddapc_methods.F:113
cp_ddapc_methods::build_a_matrix
subroutine, public build_a_matrix(am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
Computes the A matrix for the solution of the linear system.
Definition cp_ddapc_methods.F:191
cp_ddapc_methods::prep_g_dot_rvec_sin_cos
subroutine, public prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
precompute sin(g.r) and cos(g.r) for quicker evaluations of sin(g.(r1-r2)) and cos(g....
Definition cp_ddapc_methods.F:492
cp_ddapc_methods::build_der_a_matrix_rows
subroutine, public build_der_a_matrix_rows(dam, gfunc, w, particle_set, radii, rho_tot_g, gcut, iparticle0, nparticles, g_dot_rvec_sin, g_dot_rvec_cos)
Computes the derivative of the A matrix for the evaluation of the Pulay forces.
Definition cp_ddapc_methods.F:369
cp_ddapc_types
contains information regarding the decoupling/recoupling method of Bloechl
Definition cp_ddapc_types.F:12
cp_ddapc_types::cp_ddapc_create
subroutine, public cp_ddapc_create(cp_para_env, cp_ddapc_env, cp_ddapc_ewald, particle_set, radii, cell, super_cell, rho_tot_g, gcut, iw2, vol, force_env_section)
...
Definition cp_ddapc_types.F:99
cp_ddapc_util
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition cp_ddapc_util.F:15
cp_ddapc_util::cp_ddapc_init
subroutine, public cp_ddapc_init(qs_env)
Initialize the cp_ddapc_environment.
Definition cp_ddapc_util.F:80
cp_ddapc_util::get_ddapc
recursive subroutine, public get_ddapc(qs_env, calc_force, density_fit_section, density_type, qout1, qout2, out_radii, dq_out, ext_rho_tot_g, itype_of_density, iwc)
Computes the Density Derived Atomic Point Charges.
Definition cp_ddapc_util.F:203
cp_ddapc_util::modify_hartree_pot
subroutine, public modify_hartree_pot(v_hartree_gspace, density_fit_section, particle_set, m, ami, radii, charges)
Modify the Hartree potential.
Definition cp_ddapc_util.F:619
cp_ddapc_util::restraint_functional_potential
subroutine, public restraint_functional_potential(v_hartree_gspace, density_fit_section, particle_set, ami, radii, charges, ddapc_restraint_control, energy_res)
modify hartree potential to handle restraints in DDAPC scheme
Definition cp_ddapc_util.F:542
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_spin_density
integer, parameter, public do_spin_density
Definition input_constants.F:42
input_constants::do_full_density
integer, parameter, public do_full_density
Definition input_constants.F:42
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_charges_types
container for information about total charges on the grids
Definition qs_charges_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:526
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::ddapc_restraint_type
Definition cp_control_types.F:309
cp_control_types::dft_control_type
Definition cp_control_types.F:598
cp_ddapc_types::cp_ddapc_type
Definition cp_ddapc_types.F:53
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_charges_types::qs_charges_type
Container for information about total charges on the grids.
Definition qs_charges_types.F:41
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62