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