(git:50ddb19)
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-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! **************************************************************************************************
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)
359 bv(:) = bv(:)/vol
360 CALL rho_tot_g%pw_grid%para%group%sum(bv)
361 c1 = dot_product(cv, matmul(cp_ddapc_env%AmI, bv)) - ch_dens
362 c1 = c1/cp_ddapc_env%c0
363 qv(:) = -matmul(cp_ddapc_env%AmI, (bv - c1*cv))
364 j = 0
365 qtot = 0.0_dp
366 DO i = 1, ndim, num_gauss
367 j = j + 1
368 DO ii = 1, num_gauss
369 qtot(j) = qtot(j) + qv((i - 1) + ii)
370 END DO
371 END DO
372 IF (PRESENT(qout1)) THEN
373 IF (ASSOCIATED(qout1)) THEN
374 cpassert(SIZE(qout1) == SIZE(qv))
375 ELSE
376 ALLOCATE (qout1(SIZE(qv)))
377 END IF
378 qout1 = qv
379 END IF
380 IF (PRESENT(qout2)) THEN
381 IF (ASSOCIATED(qout2)) THEN
382 cpassert(SIZE(qout2) == SIZE(qtot))
383 ELSE
384 ALLOCATE (qout2(SIZE(qtot)))
385 END IF
386 qout2 = qtot
387 END IF
388 CALL print_atomic_charges(particle_set, qs_kind_set, iw, title=" DDAP "// &
389 trim(type_of_density)//" charges:", atomic_charges=qtot)
390 CALL timestop(handle2)
391 !
392 ! If requested evaluate also the correction to derivatives due to Pulay Forces
393 !
394 IF (need_f) THEN
395 CALL timeset(routinen//"-forces", handle3)
396 IF (iw > 0) THEN
397 WRITE (iw, '(T3,A)') " Evaluating DDAPC atomic derivatives .."
398 END IF
399 ALLOCATE (dam(ndim, ndim, 3))
400 ALLOCATE (dbv(ndim, 3))
401 ALLOCATE (dqv(ndim, SIZE(particle_set), 3))
402 !NB refactor math in inner loop - no more dqv0, but new temporaries instead
403 ALLOCATE (cvt_ami(ndim))
404 ALLOCATE (cvt_ami_damj(ndim))
405 ALLOCATE (tv(ndim, SIZE(particle_set), 3))
406 ALLOCATE (ami_cv(ndim))
407 cvt_ami(:) = matmul(cv, cp_ddapc_env%AmI)
408 ami_cv(:) = matmul(cp_ddapc_env%AmI, cv)
409 ALLOCATE (damj_qv(ndim))
410 ALLOCATE (ami_bv(ndim))
411 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 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 DO iparticle0 = 1, SIZE(particle_set), npset
419 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 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 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 dbv(:, :) = 0.0_dp
440 CALL build_der_b_vector(dbv, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
441 particle_set, radii, rho_tot_g, gcut, iparticle)
442 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 DO j = 1, 3
449 !NB dAmj is actually pretty sparse - one block of cols + one block of rows - use this here:
450 pmin = (iparticle - 1)*SIZE(radii) + 1
451 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 IF (pmin > 1) THEN
455 damj_qv(:pmin - 1) = matmul(transpose(dam(pmin:pmax, :pmin - 1, j)), qv(pmin:pmax))
456 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 damj_qv(pmin:pmax) = matmul(dam(pmin:pmax, :, j), qv(:))
460 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 IF (pmax < SIZE(particle_set)*SIZE(radii)) THEN
464 damj_qv(pmax + 1:) = matmul(transpose(dam(pmin:pmax, pmax + 1:, j)), qv(pmin:pmax))
465 cvt_ami_damj(pmax + 1:) = matmul(transpose(dam(pmin:pmax, pmax + 1:, j)), cvt_ami(pmin:pmax))
466 END IF
467 damj_qv(:) = damj_qv(:)/(vol*vol)
468 cvt_ami_damj(:) = cvt_ami_damj(:)/(vol*vol)
469 c3 = dot_product(cvt_ami_damj, ami_bv) - dot_product(cvt_ami, dbv(:, j)) - c1*dot_product(cvt_ami_damj, ami_cv)
470 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 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 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 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 cpassert(PRESENT(dq_out))
486 IF (.NOT. ASSOCIATED(dq_out)) THEN
487 ALLOCATE (dq_out(SIZE(dqv, 1), SIZE(dqv, 2), SIZE(dqv, 3)))
488 ELSE
489 cpassert(SIZE(dqv, 1) == SIZE(dq_out, 1))
490 cpassert(SIZE(dqv, 2) == SIZE(dq_out, 2))
491 cpassert(SIZE(dqv, 3) == SIZE(dq_out, 3))
492 END IF
493 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 DEALLOCATE (dqv)
500 DEALLOCATE (dam)
501 DEALLOCATE (dbv)
502 !NB deallocate new temporaries
503 DEALLOCATE (cvt_ami)
504 DEALLOCATE (cvt_ami_damj)
505 DEALLOCATE (ami_cv)
506 DEALLOCATE (tv)
507 DEALLOCATE (damj_qv)
508 DEALLOCATE (ami_bv)
509 CALL timestop(handle3)
510 END IF
511 !
512 ! End of charge fit
513 !
514 DEALLOCATE (radii)
515 DEALLOCATE (bv)
516 DEALLOCATE (cv)
517 DEALLOCATE (qv)
518 DEALLOCATE (qtot)
519 IF (.NOT. PRESENT(iwc)) THEN
520 CALL cp_print_key_finished_output(iw, logger, density_fit_section, &
521 "PROGRAM_RUN_INFO")
522 END IF
523 CALL auxbas_pool%give_back_pw(rho_tot_g)
524 CALL timestop(handle)
525 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
559
560 CALL timeset(routinen, handle)
561 n_gauss = SIZE(radii)
562 ALLOCATE (cv(n_gauss*SIZE(particle_set)))
563 ALLOCATE (uv(n_gauss*SIZE(particle_set)))
564 uv = 0.0_dp
565 CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
566 charges, energy_res)
567 !
568 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
569 gcut2 = gcut*gcut
570 associate(pw_grid => v_hartree_gspace%pw_grid)
571 vol = pw_grid%vol
572 cv = 1.0_dp/vol
573 sfac = -1.0_dp/vol
574 fac = dot_product(cv, matmul(ami, cv))
575 fac2 = dot_product(cv, matmul(ami, uv))
576 cv(:) = uv - cv*fac2/fac
577 cv(:) = matmul(ami, cv)
578 IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
579 DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
580 g2 = pw_grid%gsq(ig)
581 w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
582 IF (g2 > gcut2) EXIT
583 gvec = pw_grid%g(:, ig)
584 g_corr = 0.0_dp
585 idim = 0
586 DO iparticle = 1, SIZE(particle_set)
587 DO igauss = 1, SIZE(radii)
588 idim = idim + 1
589 rc = radii(igauss)
590 rc2 = rc*rc
591 rvec = particle_set(iparticle)%r
592 arg = dot_product(gvec, rvec)
593 phase = cmplx(cos(arg), -sin(arg), kind=dp)
594 gfunc = exp(-g2*rc2/4.0_dp)
595 g_corr = g_corr + gfunc*cv(idim)*phase
596 END DO
597 END DO
598 g_corr = g_corr*w
599 v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/vol
600 END DO
601 END associate
602 CALL timestop(handle)
603 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cv, uv
633
634 CALL timeset(routinen, handle)
635 CALL section_vals_val_get(density_fit_section, "GCUT", r_val=gcut)
636 gcut2 = gcut*gcut
637 associate(pw_grid => v_hartree_gspace%pw_grid)
638 vol = pw_grid%vol
639 ALLOCATE (cv(SIZE(m, 1)))
640 ALLOCATE (uv(SIZE(m, 1)))
641 cv = 1.0_dp/vol
642 uv(:) = matmul(m, charges)
643 sfac = -1.0_dp/vol
644 fac = dot_product(cv, matmul(ami, cv))
645 fac2 = dot_product(cv, matmul(ami, uv))
646 cv(:) = uv - cv*fac2/fac
647 cv(:) = matmul(ami, cv)
648 IF (pw_grid%have_g0) v_hartree_gspace%array(1) = v_hartree_gspace%array(1) + sfac*fac2/fac
649 DO ig = pw_grid%first_gne0, pw_grid%ngpts_cut_local
650 g2 = pw_grid%gsq(ig)
651 w = 4.0_dp*pi*(g2 - gcut2)**2.0_dp/(g2*gcut2)
652 IF (g2 > gcut2) EXIT
653 gvec = pw_grid%g(:, ig)
654 g_corr = 0.0_dp
655 idim = 0
656 DO iparticle = 1, SIZE(particle_set)
657 DO igauss = 1, SIZE(radii)
658 idim = idim + 1
659 rc = radii(igauss)
660 rc2 = rc*rc
661 rvec = particle_set(iparticle)%r
662 arg = dot_product(gvec, rvec)
663 phase = cmplx(cos(arg), -sin(arg), kind=dp)
664 gfunc = exp(-g2*rc2/4.0_dp)
665 g_corr = g_corr + gfunc*cv(idim)*phase
666 END DO
667 END DO
668 g_corr = g_corr*w
669 v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + sfac*g_corr/vol
670 END DO
671 END associate
672 CALL timestop(handle)
673 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: bv1, bv2, ddbv
706 TYPE(cp_ddapc_type), POINTER :: cp_ddapc_env
707
708 NULLIFY (cp_ddapc_env)
709 CALL timeset(routinen, handle)
710 dx = 0.01_dp
711 ndim = SIZE(particle_set)*SIZE(radii)
712 ALLOCATE (bv1(ndim))
713 ALLOCATE (bv2(ndim))
714 ALLOCATE (ddbv(ndim))
715 rvec = particle_set(iparticle)%r
716 cp_ddapc_env => qs_env%cp_ddapc_env
717 DO i = 1, 3
718 bv1(:) = 0.0_dp
719 bv2(:) = 0.0_dp
720 particle_set(iparticle)%r(i) = rvec(i) + dx
721 CALL build_b_vector(bv1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
722 particle_set, radii, rho_tot_g, gcut)
723 bv1(:) = bv1(:)/vol
724 CALL rho_tot_g%pw_grid%para%group%sum(bv1)
725 particle_set(iparticle)%r(i) = rvec(i) - dx
726 CALL build_b_vector(bv2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
727 particle_set, radii, rho_tot_g, gcut)
728 bv2(:) = bv2(:)/vol
729 CALL rho_tot_g%pw_grid%para%group%sum(bv2)
730 ddbv(:) = (bv1(:) - bv2(:))/(2.0_dp*dx)
731 DO kk = 1, SIZE(ddbv)
732 IF (ddbv(kk) > 1.0e-8_dp) THEN
733 v0 = abs(dbv(kk, i) - ddbv(kk))/ddbv(kk)*100.0_dp
734 WRITE (*, *) "Error % on B ::", v0
735 IF (v0 > 0.1_dp) THEN
736 WRITE (*, '(A,2I5,2F15.9)') "ERROR IN DERIVATIVE OF B VECTOR, IPARTICLE, ICOORD:", iparticle, i, &
737 dbv(kk, i), ddbv(kk)
738 cpabort("Error on B large than 0.1")
739 END IF
740 END IF
741 END DO
742 particle_set(iparticle)%r = rvec
743 END DO
744 DEALLOCATE (bv1)
745 DEALLOCATE (bv2)
746 DEALLOCATE (ddbv)
747 CALL timestop(handle)
748 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: am1, am2, ddam, g_dot_rvec_cos, &
781 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 NULLIFY (cp_ddapc_env)
787 CALL timeset(routinen, handle)
788 dx = 0.01_dp
789 ndim = SIZE(particle_set)*SIZE(radii)
790 ALLOCATE (am1(ndim, ndim))
791 ALLOCATE (am2(ndim, ndim))
792 ALLOCATE (ddam(ndim, ndim))
793 rvec = particle_set(iparticle)%r
794 cp_ddapc_env => qs_env%cp_ddapc_env
795 CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
796 DO i = 1, 3
797 am1 = 0.0_dp
798 am2 = 0.0_dp
799 particle_set(iparticle)%r(i) = rvec(i) + dx
800 CALL build_a_matrix(am1, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
801 particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
802 am1(:, :) = am1(:, :)/(vol*vol)
803 CALL rho_tot_g%pw_grid%para%group%sum(am1)
804 particle_set(iparticle)%r(i) = rvec(i) - dx
805 CALL build_a_matrix(am2, cp_ddapc_env%gfunc, cp_ddapc_env%w, &
806 particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
807 am2(:, :) = am2(:, :)/(vol*vol)
808 CALL rho_tot_g%pw_grid%para%group%sum(am2)
809 ddam(:, :) = (am1 - am2)/(2.0_dp*dx)
810 DO kk = 1, SIZE(ddam, 1)
811 DO ll = 1, SIZE(ddam, 2)
812 IF (ddam(kk, ll) > 1.0e-8_dp) THEN
813 v0 = abs(dam(kk, ll, i) - ddam(kk, ll))/ddam(kk, ll)*100.0_dp
814 WRITE (*, *) "Error % on A ::", v0, am1(kk, ll), am2(kk, ll), iparticle, i, kk, ll
815 IF (v0 > 0.1_dp) THEN
816 WRITE (*, '(A,4I5,2F15.9)') "ERROR IN DERIVATIVE OF A MATRIX, IPARTICLE, ICOORD:", iparticle, i, kk, ll, &
817 dam(kk, ll, i), ddam(kk, ll)
818 cpabort("Error on A larger than 0.1")
819 END IF
820 END IF
821 END DO
822 END DO
823 particle_set(iparticle)%r = rvec
824 END DO
825 CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
826 DEALLOCATE (am1)
827 DEALLOCATE (am2)
828 DEALLOCATE (ddam)
829 CALL timestop(handle)
830 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ddqv
860 REAL(kind=dp), DIMENSION(:), POINTER :: qtot1, qtot2
861
862 CALL timeset(routinen, handle)
863 WRITE (*, *) "DEBUG_CHARGE_ROUTINE"
864 ndim = SIZE(particle_set)*SIZE(radii)
865 NULLIFY (qtot1, qtot2)
866 ALLOCATE (qtot1(ndim))
867 ALLOCATE (qtot2(ndim))
868 ALLOCATE (ddqv(ndim))
869 !
870 dx = 0.001_dp
871 DO iparticle = 1, SIZE(particle_set)
872 rvec = particle_set(iparticle)%r
873 DO i = 1, 3
874 particle_set(iparticle)%r(i) = rvec(i) + dx
875 CALL get_ddapc(qs_env, .false., density_fit_section, qout1=qtot1, &
876 ext_rho_tot_g=rho_tot_g, itype_of_density=type_of_density)
877 particle_set(iparticle)%r(i) = rvec(i) - dx
878 CALL get_ddapc(qs_env, .false., density_fit_section, qout1=qtot2, &
879 ext_rho_tot_g=rho_tot_g, itype_of_density=type_of_density)
880 ddqv(:) = (qtot1 - qtot2)/(2.0_dp*dx)
881 DO kk = 1, SIZE(qtot1) - 1, SIZE(radii)
882 IF (any(ddqv(kk:kk + 2) > 1.0e-8_dp)) THEN
883 WRITE (*, '(A,2F12.6,F12.2)') "Error :", sum(dqv(kk:kk + 2, iparticle, i)), sum(ddqv(kk:kk + 2)), &
884 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 particle_set(iparticle)%r = rvec
888 END DO
889 END DO
890 !
891 DEALLOCATE (qtot1)
892 DEALLOCATE (qtot2)
893 DEALLOCATE (ddqv)
894 CALL timestop(handle)
895 END SUBROUTINE debug_charge
896
897END 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:56
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:225
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_transfer
Definition pw_methods.F:304
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:338
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:620
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:543
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:1073
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_spin_density
integer, parameter, public do_spin_density
Definition input_constants.F:43
input_constants::do_full_density
integer, parameter, public do_full_density
Definition input_constants.F:43
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:735
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:1052
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, mimic, 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, matrix_vhxc, 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, xcint_weights, 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:534
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:60
cp_control_types::ddapc_restraint_type
Definition cp_control_types.F:432
cp_control_types::dft_control_type
Definition cp_control_types.F:757
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:746
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:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:185
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