(git:34ef472)
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-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Density Derived atomic point charges from a QM calculation
10 !> (see Bloechl, J. Chem. Phys. Vol. 103 pp. 7422-7428)
11 !> \par History
12 !> 08.2005 created [tlaino]
13 !> \author Teodoro Laino
14 ! **************************************************************************************************
15 MODULE cp_ddapc_util
16 
17  USE atomic_charges, ONLY: print_atomic_charges
18  USE cell_types, ONLY: cell_type
19  USE cp_control_types, ONLY: ddapc_restraint_type,&
20  dft_control_type
21  USE cp_ddapc_forces, ONLY: evaluate_restraint_functional
22  USE cp_ddapc_methods, ONLY: build_a_matrix,&
23  build_b_vector,&
24  build_der_a_matrix_rows,&
25  build_der_b_vector,&
26  cleanup_g_dot_rvec_sin_cos,&
27  prep_g_dot_rvec_sin_cos
28  USE cp_ddapc_types, ONLY: cp_ddapc_create,&
29  cp_ddapc_type
30  USE cp_log_handling, ONLY: cp_get_default_logger,&
31  cp_logger_type
32  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
33  cp_print_key_unit_nr
34  USE input_constants, ONLY: do_full_density,&
35  do_spin_density
36  USE input_section_types, ONLY: section_vals_get_subs_vals,&
37  section_vals_type,&
38  section_vals_val_get
39  USE kinds, ONLY: default_string_length,&
40  dp
41  USE mathconstants, ONLY: pi
42  USE message_passing, ONLY: mp_para_env_type
43  USE particle_types, ONLY: particle_type
44  USE pw_env_types, ONLY: pw_env_get,&
45  pw_env_type
46  USE pw_methods, ONLY: pw_axpy,&
47  pw_copy,&
48  pw_transfer
49  USE pw_pool_types, ONLY: pw_pool_type
50  USE pw_types, ONLY: pw_c1d_gs_type,&
51  pw_r3d_rs_type
52  USE qs_charges_types, ONLY: qs_charges_type
53  USE qs_environment_types, ONLY: get_qs_env,&
54  qs_environment_type
55  USE qs_kind_types, ONLY: qs_kind_type
56  USE qs_rho_types, ONLY: qs_rho_get,&
57  qs_rho_type
58 #include "./base/base_uses.f90"
59 
60  IMPLICIT NONE
61  PRIVATE
62 
63  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
64  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_util'
65  PUBLIC :: get_ddapc, &
66  restraint_functional_potential, &
67  modify_hartree_pot, &
68  cp_ddapc_init
69 
70 CONTAINS
71 
72 ! **************************************************************************************************
73 !> \brief Initialize the cp_ddapc_environment
74 !> \param qs_env ...
75 !> \par History
76 !> 08.2005 created [tlaino]
77 !> \author Teodoro Laino
78 ! **************************************************************************************************
79  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 
887 END 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:48
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:195
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_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_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_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:803
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:1013
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:724
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:1040
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
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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