(git:1f285aa)
qs_cdft_methods.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 Subroutines for building CDFT constraints
10 !> \par History
11 !> separated from et_coupling [03.2017]
12 !> \author Nico Holmberg [03.2017]
13 ! **************************************************************************************************
14 MODULE qs_cdft_methods
15  USE ao_util, ONLY: exp_radius_very_extended
16  USE atomic_kind_types, ONLY: atomic_kind_type,&
17  get_atomic_kind,&
18  get_atomic_kind_set
19  USE cell_types, ONLY: cell_type,&
20  pbc
21  USE cp_control_types, ONLY: dft_control_type
22  USE cp_log_handling, ONLY: cp_get_default_logger,&
23  cp_logger_type
24  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
25  cp_print_key_unit_nr
26  USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw
27  USE grid_api, ONLY: grid_func_ab,&
28  collocate_pgf_product
29  USE hirshfeld_types, ONLY: hirshfeld_type
30  USE input_constants, ONLY: cdft_alpha_constraint,&
31  cdft_beta_constraint,&
32  cdft_charge_constraint,&
33  cdft_magnetization_constraint,&
34  outer_scf_becke_constraint,&
35  outer_scf_hirshfeld_constraint
36  USE input_section_types, ONLY: section_vals_get_subs_vals,&
37  section_vals_type
38  USE kahan_sum, ONLY: accurate_dot_product
39  USE kinds, ONLY: dp
40  USE message_passing, ONLY: mp_para_env_type
41  USE particle_types, ONLY: particle_type
42  USE pw_env_types, ONLY: pw_env_get,&
43  pw_env_type
44  USE pw_methods, ONLY: pw_axpy,&
45  pw_copy,&
46  pw_integral_ab,&
47  pw_integrate_function,&
48  pw_set,&
49  pw_zero
50  USE pw_pool_types, ONLY: pw_pool_type
51  USE pw_types, ONLY: pw_r3d_rs_type
52  USE qs_cdft_types, ONLY: becke_constraint_type,&
53  cdft_control_type,&
54  cdft_group_type,&
55  hirshfeld_constraint_type
56  USE qs_cdft_utils, ONLY: becke_constraint_init,&
57  cdft_constraint_print,&
58  cdft_print_hirshfeld_density,&
59  hfun_scale,&
60  hirshfeld_constraint_init
61  USE qs_energy_types, ONLY: qs_energy_type
62  USE qs_environment_types, ONLY: get_qs_env,&
63  qs_environment_type
64  USE qs_force_types, ONLY: qs_force_type
65  USE qs_kind_types, ONLY: get_qs_kind,&
66  qs_kind_type
67  USE qs_rho0_types, ONLY: get_rho0_mpole,&
68  mpole_rho_atom,&
69  rho0_mpole_type
70  USE qs_rho_types, ONLY: qs_rho_get,&
71  qs_rho_type
72  USE qs_subsys_types, ONLY: qs_subsys_get,&
73  qs_subsys_type
74  USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
75  realspace_grid_type,&
76  rs_grid_create,&
77  rs_grid_release,&
78  rs_grid_zero,&
79  transfer_rs2pw
80 #include "./base/base_uses.f90"
81 
82  IMPLICIT NONE
83 
84  PRIVATE
85 
86  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_methods'
87  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
88 
89 ! *** Public subroutines ***
90 
91  PUBLIC :: becke_constraint, hirshfeld_constraint
92 
93 CONTAINS
94 
95 ! **************************************************************************************************
96 !> \brief Driver routine for calculating a Becke constraint
97 !> \param qs_env the qs_env where to build the constraint
98 !> \param calc_pot if the potential needs to be recalculated or just integrated
99 !> \param calculate_forces logical if potential has to be calculated or only_energy
100 !> \par History
101 !> Created 01.2007 [fschiff]
102 !> Extended functionality 12/15-12/16 [Nico Holmberg]
103 ! **************************************************************************************************
104  SUBROUTINE becke_constraint(qs_env, calc_pot, calculate_forces)
105  TYPE(qs_environment_type), POINTER :: qs_env
106  LOGICAL :: calc_pot, calculate_forces
107 
108  CHARACTER(len=*), PARAMETER :: routinen = 'becke_constraint'
109 
110  INTEGER :: handle
111  TYPE(cdft_control_type), POINTER :: cdft_control
112  TYPE(dft_control_type), POINTER :: dft_control
113 
114  CALL timeset(routinen, handle)
115  CALL get_qs_env(qs_env, dft_control=dft_control)
116  cdft_control => dft_control%qs_control%cdft_control
117  IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_becke_constraint) THEN
118  IF (calc_pot) THEN
119  ! Initialize the Becke constraint environment
120  CALL becke_constraint_init(qs_env)
121  ! Calculate the Becke weight function and possibly the gradients
122  CALL becke_constraint_low(qs_env)
123  END IF
124  ! Integrate the Becke constraint
125  CALL cdft_constraint_integrate(qs_env)
126  ! Calculate forces
127  IF (calculate_forces) CALL cdft_constraint_force(qs_env)
128  END IF
129  CALL timestop(handle)
130 
131  END SUBROUTINE becke_constraint
132 
133 ! **************************************************************************************************
134 !> \brief Low level routine to build a Becke weight function and its gradients
135 !> \param qs_env the qs_env where to build the constraint
136 !> \param just_gradients optional logical which determines if only the gradients should be calculated
137 !> \par History
138 !> Created 03.2017 [Nico Holmberg]
139 ! **************************************************************************************************
140  SUBROUTINE becke_constraint_low(qs_env, just_gradients)
141  TYPE(qs_environment_type), POINTER :: qs_env
142  LOGICAL, OPTIONAL :: just_gradients
143 
144  CHARACTER(len=*), PARAMETER :: routinen = 'becke_constraint_low'
145 
146  INTEGER :: handle, i, iatom, igroup, ind(3), ip, j, &
147  jatom, jp, k, natom, np(3), nskipped
148  INTEGER, ALLOCATABLE, DIMENSION(:) :: catom
149  INTEGER, DIMENSION(2, 3) :: bo, bo_conf
150  LOGICAL :: in_memory, my_just_gradients
151  LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint, skip_me
152  LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: atom_in_group
153  REAL(kind=dp) :: dist1, dist2, dmyexp, dvol, eps_cavity, &
154  my1, my1_homo, myexp, sum_cell_f_all, &
155  th, tmp_const
156  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cell_functions, ds_dr_i, ds_dr_j, &
157  sum_cell_f_group
158  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_sum_pm_dr, dp_i_dri
159  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dp_i_drj
160  REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dmy_dr_i, dmy_dr_j, &
161  dr, dr1_r2, dr_i_dr, dr_ij_dr, &
162  dr_j_dr, grid_p, r, r1, shift
163  REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
164  TYPE(becke_constraint_type), POINTER :: becke_control
165  TYPE(cdft_control_type), POINTER :: cdft_control
166  TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
167  TYPE(cell_type), POINTER :: cell
168  TYPE(dft_control_type), POINTER :: dft_control
169  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
170  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge
171 
172  NULLIFY (cutoffs, cell, dft_control, particle_set, group, charge, cdft_control)
173  CALL timeset(routinen, handle)
174  ! Get simulation environment
175  CALL get_qs_env(qs_env, &
176  cell=cell, &
177  particle_set=particle_set, &
178  natom=natom, &
179  dft_control=dft_control)
180  cdft_control => dft_control%qs_control%cdft_control
181  becke_control => cdft_control%becke_control
182  group => cdft_control%group
183  cutoffs => becke_control%cutoffs
184  IF (cdft_control%atomic_charges) &
185  charge => cdft_control%charge
186  in_memory = .false.
187  IF (cdft_control%save_pot) THEN
188  in_memory = becke_control%in_memory
189  END IF
190  eps_cavity = becke_control%eps_cavity
191  ! Decide if only gradients need to be calculated
192  my_just_gradients = .false.
193  IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
194  IF (my_just_gradients) THEN
195  in_memory = .true.
196  ! Pairwise distances need to be recalculated
197  IF (becke_control%vector_buffer%store_vectors) THEN
198  ALLOCATE (becke_control%vector_buffer%distances(natom))
199  ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom))
200  IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom))
201  ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom))
202  END IF
203  ALLOCATE (becke_control%vector_buffer%R12(natom, natom))
204  DO i = 1, 3
205  cell_v(i) = cell%hmat(i, i)
206  END DO
207  DO iatom = 1, natom - 1
208  DO jatom = iatom + 1, natom
209  r = particle_set(iatom)%r
210  r1 = particle_set(jatom)%r
211  DO i = 1, 3
212  r(i) = modulo(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
213  r1(i) = modulo(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
214  END DO
215  dist_vec = (r - r1) - anint((r - r1)/cell_v)*cell_v
216  IF (becke_control%vector_buffer%store_vectors) THEN
217  becke_control%vector_buffer%position_vecs(:, iatom) = r(:)
218  IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:)
219  IF (in_memory) THEN
220  becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
221  becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
222  END IF
223  END IF
224  becke_control%vector_buffer%R12(iatom, jatom) = sqrt(dot_product(dist_vec, dist_vec))
225  becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom)
226  END DO
227  END DO
228  END IF
229  ALLOCATE (catom(cdft_control%natoms))
230  IF (cdft_control%save_pot .OR. &
231  becke_control%cavity_confine .OR. &
232  becke_control%should_skip) THEN
233  ALLOCATE (is_constraint(natom))
234  is_constraint = .false.
235  END IF
236  ! This boolean is needed to prevent calculation of atom pairs ji when the pair ij has
237  ! already been calculated (data for pair ji is set using symmetry)
238  ! With gradient precomputation, symmetry exploited for both weight function and gradients
239  ALLOCATE (skip_me(natom))
240  DO i = 1, cdft_control%natoms
241  catom(i) = cdft_control%atoms(i)
242  ! Notice that here is_constraint=.TRUE. also for dummy atoms to properly compute their Becke charges
243  ! A subsequent check (atom_in_group) ensures that the gradients of these dummy atoms are correct
244  IF (cdft_control%save_pot .OR. &
245  becke_control%cavity_confine .OR. &
246  becke_control%should_skip) &
247  is_constraint(catom(i)) = .true.
248  END DO
249  bo = group(1)%weight%pw_grid%bounds_local
250  dvol = group(1)%weight%pw_grid%dvol
251  dr = group(1)%weight%pw_grid%dr
252  np = group(1)%weight%pw_grid%npts
253  shift = -real(modulo(np, 2), dp)*dr/2.0_dp
254  DO i = 1, 3
255  cell_v(i) = cell%hmat(i, i)
256  END DO
257  ! If requested, allocate storage for gradients
258  IF (in_memory) THEN
259  bo_conf = bo
260  ! With confinement active, we dont need to store gradients outside
261  ! the confinement bounds since they vanish for all particles
262  IF (becke_control%cavity_confine) THEN
263  bo_conf(1, 3) = becke_control%confine_bounds(1)
264  bo_conf(2, 3) = becke_control%confine_bounds(2)
265  END IF
266  ALLOCATE (atom_in_group(SIZE(group), natom))
267  atom_in_group = .false.
268  DO igroup = 1, SIZE(group)
269  ALLOCATE (group(igroup)%gradients(3*natom, bo_conf(1, 1):bo_conf(2, 1), &
270  bo_conf(1, 2):bo_conf(2, 2), &
271  bo_conf(1, 3):bo_conf(2, 3)))
272  group(igroup)%gradients = 0.0_dp
273  ALLOCATE (group(igroup)%d_sum_const_dR(3, natom))
274  group(igroup)%d_sum_const_dR = 0.0_dp
275  DO ip = 1, SIZE(group(igroup)%atoms)
276  atom_in_group(igroup, group(igroup)%atoms(ip)) = .true.
277  END DO
278  END DO
279  END IF
280  ! Allocate remaining work
281  ALLOCATE (sum_cell_f_group(SIZE(group)))
282  ALLOCATE (cell_functions(natom))
283  IF (in_memory) THEN
284  ALLOCATE (ds_dr_j(3))
285  ALLOCATE (ds_dr_i(3))
286  ALLOCATE (d_sum_pm_dr(3, natom))
287  ALLOCATE (dp_i_drj(3, natom, natom))
288  ALLOCATE (dp_i_dri(3, natom))
289  th = 1.0e-8_dp
290  END IF
291  ! Build constraint
292  DO k = bo(1, 1), bo(2, 1)
293  DO j = bo(1, 2), bo(2, 2)
294  DO i = bo(1, 3), bo(2, 3)
295  ! If the grid point is too far from all constraint atoms and cavity confinement is active,
296  ! we can skip this grid point as it does not contribute to the weight or gradients
297  IF (becke_control%cavity_confine) THEN
298  IF (becke_control%cavity%array(k, j, i) < eps_cavity) cycle
299  END IF
300  ind = (/k, j, i/)
301  grid_p(1) = k*dr(1) + shift(1)
302  grid_p(2) = j*dr(2) + shift(2)
303  grid_p(3) = i*dr(3) + shift(3)
304  nskipped = 0
305  cell_functions = 1.0_dp
306  skip_me = .false.
307  IF (becke_control%vector_buffer%store_vectors) becke_control%vector_buffer%distances = 0.0_dp
308  IF (in_memory) THEN
309  d_sum_pm_dr = 0.0_dp
310  DO igroup = 1, SIZE(group)
311  group(igroup)%d_sum_const_dR = 0.0_dp
312  END DO
313  dp_i_dri = 0.0_dp
314  END IF
315  ! Iterate over all atoms in the system
316  DO iatom = 1, natom
317  IF (skip_me(iatom)) THEN
318  cell_functions(iatom) = 0.0_dp
319  IF (becke_control%should_skip) THEN
320  IF (is_constraint(iatom)) nskipped = nskipped + 1
321  IF (nskipped == cdft_control%natoms) THEN
322  IF (in_memory) THEN
323  IF (becke_control%cavity_confine) THEN
324  becke_control%cavity%array(k, j, i) = 0.0_dp
325  END IF
326  END IF
327  EXIT
328  END IF
329  END IF
330  cycle
331  END IF
332  IF (becke_control%vector_buffer%store_vectors) THEN
333  IF (becke_control%vector_buffer%distances(iatom) .EQ. 0.0_dp) THEN
334  r = becke_control%vector_buffer%position_vecs(:, iatom)
335  dist_vec = (r - grid_p) - anint((r - grid_p)/cell_v)*cell_v
336  dist1 = sqrt(dot_product(dist_vec, dist_vec))
337  becke_control%vector_buffer%distance_vecs(:, iatom) = dist_vec
338  becke_control%vector_buffer%distances(iatom) = dist1
339  ELSE
340  dist_vec = becke_control%vector_buffer%distance_vecs(:, iatom)
341  dist1 = becke_control%vector_buffer%distances(iatom)
342  END IF
343  ELSE
344  r = particle_set(iatom)%r
345  DO ip = 1, 3
346  r(ip) = modulo(r(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
347  END DO
348  dist_vec = (r - grid_p) - anint((r - grid_p)/cell_v)*cell_v
349  dist1 = sqrt(dot_product(dist_vec, dist_vec))
350  END IF
351  IF (dist1 .LE. cutoffs(iatom)) THEN
352  IF (in_memory) THEN
353  IF (dist1 .LE. th) dist1 = th
354  dr_i_dr(:) = dist_vec(:)/dist1
355  END IF
356  DO jatom = 1, natom
357  IF (jatom .NE. iatom) THEN
358  ! Using pairwise symmetry, execute block only for such j<i
359  ! that have previously not been looped over
360  ! Note that if skip_me(jatom) = .TRUE., this means that the outer
361  ! loop over iatom skipped this index when iatom=jatom, but we still
362  ! need to compute the pair for iatom>jatom
363  IF (jatom < iatom) THEN
364  IF (.NOT. skip_me(jatom)) cycle
365  END IF
366  IF (becke_control%vector_buffer%store_vectors) THEN
367  IF (becke_control%vector_buffer%distances(jatom) .EQ. 0.0_dp) THEN
368  r1 = becke_control%vector_buffer%position_vecs(:, jatom)
369  dist_vec = (r1 - grid_p) - anint((r1 - grid_p)/cell_v)*cell_v
370  dist2 = sqrt(dot_product(dist_vec, dist_vec))
371  becke_control%vector_buffer%distance_vecs(:, jatom) = dist_vec
372  becke_control%vector_buffer%distances(jatom) = dist2
373  ELSE
374  dist_vec = becke_control%vector_buffer%distance_vecs(:, jatom)
375  dist2 = becke_control%vector_buffer%distances(jatom)
376  END IF
377  ELSE
378  r1 = particle_set(jatom)%r
379  DO ip = 1, 3
380  r1(ip) = modulo(r1(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
381  END DO
382  dist_vec = (r1 - grid_p) - anint((r1 - grid_p)/cell_v)*cell_v
383  dist2 = sqrt(dot_product(dist_vec, dist_vec))
384  END IF
385  IF (in_memory) THEN
386  IF (becke_control%vector_buffer%store_vectors) THEN
387  dr1_r2 = becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom)
388  ELSE
389  dr1_r2 = (r - r1) - anint((r - r1)/cell_v)*cell_v
390  END IF
391  IF (dist2 .LE. th) dist2 = th
392  tmp_const = (becke_control%vector_buffer%R12(iatom, jatom)**3)
393  dr_ij_dr(:) = dr1_r2(:)/tmp_const
394  !derivative w.r.t. Rj
395  dr_j_dr = dist_vec(:)/dist2
396  dmy_dr_j(:) = -(dr_j_dr(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dr(:))
397  !derivative w.r.t. Ri
398  dmy_dr_i(:) = dr_i_dr(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dr(:)
399  END IF
400  ! myij
401  my1 = (dist1 - dist2)/becke_control%vector_buffer%R12(iatom, jatom)
402  IF (becke_control%adjust) THEN
403  my1_homo = my1 ! Homonuclear quantity needed for gradient
404  my1 = my1 + becke_control%aij(iatom, jatom)*(1.0_dp - my1**2)
405  END IF
406  ! f(myij)
407  myexp = 1.5_dp*my1 - 0.5_dp*my1**3
408  IF (in_memory) THEN
409  dmyexp = 1.5_dp - 1.5_dp*my1**2
410  tmp_const = (1.5_dp**2)*dmyexp*(1 - myexp**2)* &
411  (1.0_dp - ((1.5_dp*myexp - 0.5_dp*(myexp**3))**2))
412  ! d s(myij)/d R_i
413  ds_dr_i(:) = -0.5_dp*tmp_const*dmy_dr_i(:)
414  ! d s(myij)/d R_j
415  ds_dr_j(:) = -0.5_dp*tmp_const*dmy_dr_j(:)
416  IF (becke_control%adjust) THEN
417  tmp_const = 1.0_dp - 2.0_dp*my1_homo* &
418  becke_control%aij(iatom, jatom)
419  ds_dr_i(:) = ds_dr_i(:)*tmp_const
420  ! tmp_const is same for both since aij=-aji and myij=-myji
421  ds_dr_j(:) = ds_dr_j(:)*tmp_const
422  END IF
423  END IF
424  ! s(myij) = f[f(f{myij})]
425  myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
426  myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
427  tmp_const = 0.5_dp*(1.0_dp - myexp)
428  cell_functions(iatom) = cell_functions(iatom)*tmp_const
429  IF (in_memory) THEN
430  IF (abs(tmp_const) .LE. th) tmp_const = tmp_const + th
431  ! P_i independent part of dP_i/dR_i
432  dp_i_dri(:, iatom) = dp_i_dri(:, iatom) + ds_dr_i(:)/tmp_const
433  ! P_i independent part of dP_i/dR_j
434  dp_i_drj(:, iatom, jatom) = ds_dr_j(:)/tmp_const
435  END IF
436 
437  IF (dist2 .LE. cutoffs(jatom)) THEN
438  tmp_const = 0.5_dp*(1.0_dp + myexp) ! s(myji)
439  cell_functions(jatom) = cell_functions(jatom)*tmp_const
440  IF (in_memory) THEN
441  IF (abs(tmp_const) .LE. th) tmp_const = tmp_const + th
442  ! P_j independent part of dP_j/dR_i
443  ! d s(myji)/d R_i = -d s(myij)/d R_i
444  dp_i_drj(:, jatom, iatom) = -ds_dr_i(:)/tmp_const
445  ! P_j independent part of dP_j/dR_j
446  ! d s(myji)/d R_j = -d s(myij)/d R_j
447  dp_i_dri(:, jatom) = dp_i_dri(:, jatom) - ds_dr_j(:)/tmp_const
448  END IF
449  ELSE
450  skip_me(jatom) = .true.
451  END IF
452  END IF
453  END DO ! jatom
454  IF (in_memory) THEN
455  ! Final value of dP_i_dRi
456  dp_i_dri(:, iatom) = cell_functions(iatom)*dp_i_dri(:, iatom)
457  ! Update relevant sums with value
458  d_sum_pm_dr(:, iatom) = d_sum_pm_dr(:, iatom) + dp_i_dri(:, iatom)
459  IF (is_constraint(iatom)) THEN
460  DO igroup = 1, SIZE(group)
461  IF (.NOT. atom_in_group(igroup, iatom)) cycle
462  DO jp = 1, SIZE(group(igroup)%atoms)
463  IF (iatom == group(igroup)%atoms(jp)) THEN
464  ip = jp
465  EXIT
466  END IF
467  END DO
468  group(igroup)%d_sum_const_dR(1:3, iatom) = group(igroup)%d_sum_const_dR(1:3, iatom) + &
469  group(igroup)%coeff(ip)*dp_i_dri(:, iatom)
470  END DO
471  END IF
472  DO jatom = 1, natom
473  IF (jatom .NE. iatom) THEN
474  ! Final value of dP_i_dRj
475  dp_i_drj(:, iatom, jatom) = cell_functions(iatom)*dp_i_drj(:, iatom, jatom)
476  ! Update where needed
477  d_sum_pm_dr(:, jatom) = d_sum_pm_dr(:, jatom) + dp_i_drj(:, iatom, jatom)
478  IF (is_constraint(iatom)) THEN
479  DO igroup = 1, SIZE(group)
480  IF (.NOT. atom_in_group(igroup, iatom)) cycle
481  ip = -1
482  DO jp = 1, SIZE(group(igroup)%atoms)
483  IF (iatom == group(igroup)%atoms(jp)) THEN
484  ip = jp
485  EXIT
486  END IF
487  END DO
488  group(igroup)%d_sum_const_dR(1:3, jatom) = group(igroup)%d_sum_const_dR(1:3, jatom) + &
489  group(igroup)%coeff(ip)* &
490  dp_i_drj(:, iatom, jatom)
491  END DO
492  END IF
493  END IF
494  END DO
495  END IF
496  ELSE
497  cell_functions(iatom) = 0.0_dp
498  skip_me(iatom) = .true.
499  IF (becke_control%should_skip) THEN
500  IF (is_constraint(iatom)) nskipped = nskipped + 1
501  IF (nskipped == cdft_control%natoms) THEN
502  IF (in_memory) THEN
503  IF (becke_control%cavity_confine) THEN
504  becke_control%cavity%array(k, j, i) = 0.0_dp
505  END IF
506  END IF
507  EXIT
508  END IF
509  END IF
510  END IF
511  END DO !iatom
512  IF (nskipped == cdft_control%natoms) cycle
513  ! Sum up cell functions
514  sum_cell_f_group = 0.0_dp
515  DO igroup = 1, SIZE(group)
516  DO ip = 1, SIZE(group(igroup)%atoms)
517  sum_cell_f_group(igroup) = sum_cell_f_group(igroup) + group(igroup)%coeff(ip)* &
518  cell_functions(group(igroup)%atoms(ip))
519  END DO
520  END DO
521  sum_cell_f_all = 0.0_dp
522  DO ip = 1, natom
523  sum_cell_f_all = sum_cell_f_all + cell_functions(ip)
524  END DO
525  ! Gradients at (k,j,i)
526  IF (in_memory .AND. abs(sum_cell_f_all) .GT. 0.0_dp) THEN
527  DO igroup = 1, SIZE(group)
528  DO iatom = 1, natom
529  group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) = &
530  group(igroup)%d_sum_const_dR(1:3, iatom)/sum_cell_f_all - sum_cell_f_group(igroup)* &
531  d_sum_pm_dr(1:3, iatom)/(sum_cell_f_all**2)
532  END DO
533  END DO
534  END IF
535  ! Weight function(s) at (k,j,i)
536  IF (.NOT. my_just_gradients .AND. abs(sum_cell_f_all) .GT. 0.000001) THEN
537  DO igroup = 1, SIZE(group)
538  group(igroup)%weight%array(k, j, i) = sum_cell_f_group(igroup)/sum_cell_f_all
539  END DO
540  IF (cdft_control%atomic_charges) THEN
541  DO iatom = 1, cdft_control%natoms
542  charge(iatom)%array(k, j, i) = cell_functions(catom(iatom))/sum_cell_f_all
543  END DO
544  END IF
545  END IF
546  END DO
547  END DO
548  END DO
549  ! Release storage
550  IF (in_memory) THEN
551  DEALLOCATE (ds_dr_j)
552  DEALLOCATE (ds_dr_i)
553  DEALLOCATE (d_sum_pm_dr)
554  DEALLOCATE (dp_i_drj)
555  DEALLOCATE (dp_i_dri)
556  DO igroup = 1, SIZE(group)
557  DEALLOCATE (group(igroup)%d_sum_const_dR)
558  END DO
559  DEALLOCATE (atom_in_group)
560  IF (becke_control%vector_buffer%store_vectors) THEN
561  DEALLOCATE (becke_control%vector_buffer%pair_dist_vecs)
562  END IF
563  END IF
564  NULLIFY (cutoffs)
565  IF (ALLOCATED(is_constraint)) &
566  DEALLOCATE (is_constraint)
567  DEALLOCATE (catom)
568  DEALLOCATE (cell_functions)
569  DEALLOCATE (skip_me)
570  DEALLOCATE (sum_cell_f_group)
571  DEALLOCATE (becke_control%vector_buffer%R12)
572  IF (becke_control%vector_buffer%store_vectors) THEN
573  DEALLOCATE (becke_control%vector_buffer%distances)
574  DEALLOCATE (becke_control%vector_buffer%distance_vecs)
575  DEALLOCATE (becke_control%vector_buffer%position_vecs)
576  END IF
577  CALL timestop(handle)
578 
579  END SUBROUTINE becke_constraint_low
580 
581 ! **************************************************************************************************
582 !> \brief Driver routine for calculating a Hirshfeld constraint
583 !> \param qs_env ...
584 !> \param calc_pot ...
585 !> \param calculate_forces ...
586 ! **************************************************************************************************
587  SUBROUTINE hirshfeld_constraint(qs_env, calc_pot, calculate_forces)
588  TYPE(qs_environment_type), POINTER :: qs_env
589  LOGICAL :: calc_pot, calculate_forces
590 
591  CHARACTER(len=*), PARAMETER :: routinen = 'hirshfeld_constraint'
592 
593  INTEGER :: handle
594  TYPE(cdft_control_type), POINTER :: cdft_control
595  TYPE(dft_control_type), POINTER :: dft_control
596 
597  CALL timeset(routinen, handle)
598  CALL get_qs_env(qs_env, dft_control=dft_control)
599  cdft_control => dft_control%qs_control%cdft_control
600  IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_hirshfeld_constraint) THEN
601  IF (calc_pot) THEN
602  ! Initialize the Hirshfeld constraint environment
603  CALL hirshfeld_constraint_init(qs_env)
604  ! Calculate the Hirshfeld weight function and possibly the gradients
605  CALL hirshfeld_constraint_low(qs_env)
606  END IF
607  ! Integrate the Hirshfeld constraint
608  CALL cdft_constraint_integrate(qs_env)
609  ! Calculate forces
610  IF (calculate_forces) CALL cdft_constraint_force(qs_env)
611  END IF
612  CALL timestop(handle)
613 
614  END SUBROUTINE hirshfeld_constraint
615 
616 ! **************************************************************************************************
617 !> \brief Calculates Hirshfeld constraints
618 !> \param qs_env ...
619 !> \param just_gradients ...
620 ! **************************************************************************************************
621  SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
622  TYPE(qs_environment_type), POINTER :: qs_env
623  LOGICAL, OPTIONAL :: just_gradients
624 
625  CHARACTER(len=*), PARAMETER :: routinen = 'hirshfeld_constraint_low'
626 
627  INTEGER :: atom_a, atoms_memory, atoms_memory_num, handle, i, iatom, iex, igroup, ikind, &
628  ithread, j, k, natom, npme, nthread, num_atoms, num_species, numexp, subpatch_pattern
629  INTEGER, ALLOCATABLE, DIMENSION(:) :: num_species_small
630  INTEGER, DIMENSION(2, 3) :: bo
631  INTEGER, DIMENSION(3) :: lb_pw, lb_rs, npts, ub_pw, ub_rs
632  INTEGER, DIMENSION(:), POINTER :: atom_list, cores
633  LOGICAL :: my_just_gradients
634  LOGICAL, ALLOCATABLE, DIMENSION(:) :: compute_charge, is_constraint
635  REAL(kind=dp) :: alpha, coef, eps_rho_rspace, exp_eval, &
636  prefactor, radius
637  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients
638  REAL(kind=dp), DIMENSION(3) :: dr_pw, dr_rs, origin, r2, r_pbc, ra
639  REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
640  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
641  TYPE(cdft_control_type), POINTER :: cdft_control
642  TYPE(cell_type), POINTER :: cell
643  TYPE(dft_control_type), POINTER :: dft_control
644  TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control
645  TYPE(hirshfeld_type), POINTER :: hirshfeld_env
646  TYPE(mp_para_env_type), POINTER :: para_env
647  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
648  TYPE(pw_env_type), POINTER :: pw_env
649  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
650  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: pw_single_dr
651  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
652  TYPE(qs_rho_type), POINTER :: rho
653  TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
654  TYPE(realspace_grid_type) :: rs_rho_all, rs_rho_constr
655  TYPE(realspace_grid_type), ALLOCATABLE, &
656  DIMENSION(:) :: rs_single, rs_single_charge, rs_single_dr
657 
658  NULLIFY (atom_list, atomic_kind_set, dft_control, &
659  hirshfeld_env, particle_set, pw_env, auxbas_pw_pool, para_env, &
660  auxbas_rs_desc, cdft_control, pab, &
661  hirshfeld_control, cell, rho_r, rho)
662 
663  CALL timeset(routinen, handle)
664  CALL get_qs_env(qs_env, &
665  atomic_kind_set=atomic_kind_set, &
666  particle_set=particle_set, &
667  natom=natom, &
668  cell=cell, &
669  rho=rho, &
670  dft_control=dft_control, &
671  para_env=para_env, &
672  pw_env=pw_env)
673  CALL qs_rho_get(rho, rho_r=rho_r)
674 
675  num_atoms = natom
676 
677  cdft_control => dft_control%qs_control%cdft_control
678  hirshfeld_control => cdft_control%hirshfeld_control
679  hirshfeld_env => hirshfeld_control%hirshfeld_env
680 
681  ! Check if only gradient should be calculated, if gradients should be precomputed
682  my_just_gradients = .false.
683  IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
684  IF (my_just_gradients) THEN
685  cdft_control%in_memory = .true.
686  cdft_control%atomic_charges = .false.
687  hirshfeld_control%print_density = .false.
688  END IF
689 
690  ALLOCATE (coefficients(natom))
691  ALLOCATE (is_constraint(natom))
692 
693  subpatch_pattern = 0
694  eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
695  radius = 100.0_dp
696 
697  dr_pw(1) = rho_r(1)%pw_grid%dr(1)
698  dr_pw(2) = rho_r(1)%pw_grid%dr(2)
699  dr_pw(3) = rho_r(1)%pw_grid%dr(3)
700  lb_pw(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
701  ub_pw(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
702  npts = rho_r(1)%pw_grid%npts
703  origin(1) = (dr_pw(1)*npts(1))*0.5_dp
704  origin(2) = (dr_pw(2)*npts(2))*0.5_dp
705  origin(3) = (dr_pw(3)*npts(3))*0.5_dp
706 
707  CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
708  auxbas_pw_pool=auxbas_pw_pool)
709  CALL rs_grid_create(rs_rho_all, auxbas_rs_desc)
710  CALL rs_grid_zero(rs_rho_all)
711 
712  dr_rs(1) = rs_rho_all%desc%dh(1, 1)
713  dr_rs(2) = rs_rho_all%desc%dh(2, 2)
714  dr_rs(3) = rs_rho_all%desc%dh(3, 3)
715  lb_rs(1) = lbound(rs_rho_all%r(:, :, :), 1)
716  lb_rs(2) = lbound(rs_rho_all%r(:, :, :), 2)
717  lb_rs(3) = lbound(rs_rho_all%r(:, :, :), 3)
718  ub_rs(1) = ubound(rs_rho_all%r(:, :, :), 1)
719  ub_rs(2) = ubound(rs_rho_all%r(:, :, :), 2)
720  ub_rs(3) = ubound(rs_rho_all%r(:, :, :), 3)
721 
722  ! For each CDFT group
723  DO igroup = 1, SIZE(cdft_control%group)
724 
725  IF (igroup == 2 .AND. .NOT. cdft_control%in_memory) THEN
726  CALL rs_grid_zero(rs_rho_all)
727  END IF
728  bo = cdft_control%group(igroup)%weight%pw_grid%bounds_local
729 
730  ! Coefficients
731  coefficients(:) = 0.0_dp
732  is_constraint = .false.
733  DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
734  coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
735  is_constraint(cdft_control%group(igroup)%atoms(i)) = .true.
736  END DO
737 
738  ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group
739  CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc)
740  CALL rs_grid_zero(rs_rho_constr)
741 
742  ! rs_single: Gaussian density over single atoms when required
743  IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
744  ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic(cdft_control%natoms))
745  ALLOCATE (rs_single(cdft_control%natoms))
746  DO i = 1, cdft_control%natoms
747  CALL rs_grid_create(rs_single(i), auxbas_rs_desc)
748  CALL rs_grid_zero(rs_single(i))
749  END DO
750  END IF
751 
752  ! Setup pw
753  CALL pw_zero(cdft_control%group(igroup)%weight)
754 
755  CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
756  CALL pw_set(cdft_control%group(igroup)%hw_rho_total_constraint, 1.0_dp)
757 
758  IF (igroup == 1) THEN
759  CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total)
760  CALL pw_set(cdft_control%group(igroup)%hw_rho_total, 1.0_dp)
761 
762  IF (hirshfeld_control%print_density) THEN
763  DO iatom = 1, cdft_control%natoms
764  CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic(iatom))
765  CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic(iatom), 1.0_dp)
766  END DO
767  END IF
768  END IF
769 
770  IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
771  ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge(cdft_control%natoms))
772  ALLOCATE (rs_single_charge(cdft_control%natoms))
773  ALLOCATE (compute_charge(natom))
774  compute_charge = .false.
775 
776  DO i = 1, cdft_control%natoms
777  CALL rs_grid_create(rs_single_charge(i), auxbas_rs_desc)
778  CALL rs_grid_zero(rs_single_charge(i))
779  compute_charge(cdft_control%atoms(i)) = .true.
780  END DO
781 
782  DO iatom = 1, cdft_control%natoms
783  CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom))
784  CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom), 1.0_dp)
785  END DO
786  END IF
787 
788  ALLOCATE (pab(1, 1))
789  nthread = 1
790  ithread = 0
791 
792  DO ikind = 1, SIZE(atomic_kind_set)
793  numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
794  IF (numexp <= 0) cycle
795  CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
796  ALLOCATE (cores(num_species))
797 
798  DO iex = 1, numexp
799  alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
800  coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
801  npme = 0
802  cores = 0
803  DO iatom = 1, num_species
804  atom_a = atom_list(iatom)
805  ra(:) = pbc(particle_set(atom_a)%r, cell)
806  IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
807  IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
808  npme = npme + 1
809  cores(npme) = iatom
810  END IF
811  ELSE
812  npme = npme + 1
813  cores(npme) = iatom
814  END IF
815  END DO
816  DO j = 1, npme
817  iatom = cores(j)
818  atom_a = atom_list(iatom)
819  pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
820  ra(:) = pbc(particle_set(atom_a)%r, cell)
821 
822  IF (hirshfeld_control%use_atomic_cutoff) THEN
823  radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
824  ra=ra, rb=ra, rp=ra, &
825  zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
826  pab=pab, o1=0, o2=0, & ! without map_consistent
827  prefactor=1.0_dp, cutoff=0.0_dp)
828  ! IF (para_env%is_source()) PRINT *, "radius", radius
829  END IF
830 
831  IF (igroup == 1) THEN
832  CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
833  (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
834  rs_rho_all, radius=radius, &
835  ga_gb_function=grid_func_ab, use_subpatch=.true., &
836  subpatch_pattern=subpatch_pattern)
837  END IF
838 
839  IF (is_constraint(atom_a)) THEN
840  CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
841  (/0.0_dp, 0.0_dp, 0.0_dp/), coefficients(atom_a), &
842  pab, 0, 0, rs_rho_constr, &
843  radius=radius, &
844  ga_gb_function=grid_func_ab, use_subpatch=.true., &
845  subpatch_pattern=subpatch_pattern)
846  END IF
847 
848  IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
849  IF (is_constraint(atom_a)) THEN
850  DO iatom = 1, cdft_control%natoms
851  IF (atom_a == cdft_control%atoms(iatom)) EXIT
852  END DO
853  cpassert(iatom <= cdft_control%natoms)
854  CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
855  (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
856  rs_single(iatom), radius=radius, &
857  ga_gb_function=grid_func_ab, use_subpatch=.true., &
858  subpatch_pattern=subpatch_pattern)
859  END IF
860  END IF
861 
862  IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
863  IF (compute_charge(atom_a)) THEN
864  DO iatom = 1, cdft_control%natoms
865  IF (atom_a == cdft_control%atoms(iatom)) EXIT
866  END DO
867  cpassert(iatom <= cdft_control%natoms)
868  CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
869  (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
870  rs_single_charge(iatom), radius=radius, &
871  ga_gb_function=grid_func_ab, use_subpatch=.true., &
872  subpatch_pattern=subpatch_pattern)
873  END IF
874  END IF
875 
876  END DO
877  END DO
878  DEALLOCATE (cores)
879  END DO
880  DEALLOCATE (pab)
881 
882  IF (igroup == 1) THEN
883  CALL transfer_rs2pw(rs_rho_all, cdft_control%group(igroup)%hw_rho_total)
884  END IF
885 
886  CALL transfer_rs2pw(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint)
887  CALL rs_grid_release(rs_rho_constr)
888 
889  ! Calculate weight function
890  CALL hfun_scale(cdft_control%group(igroup)%weight%array, &
891  cdft_control%group(igroup)%hw_rho_total_constraint%array, &
892  cdft_control%group(1)%hw_rho_total%array, divide=.true., &
893  small=hirshfeld_control%eps_cutoff)
894 
895  ! Calculate charges
896  IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
897  DO i = 1, cdft_control%natoms
898  CALL transfer_rs2pw(rs_single_charge(i), cdft_control%group(igroup)%hw_rho_atomic_charge(i))
899  CALL hfun_scale(cdft_control%charge(i)%array, &
900  cdft_control%group(igroup)%hw_rho_atomic_charge(i)%array, &
901  cdft_control%group(igroup)%hw_rho_total%array, divide=.true., &
902  small=hirshfeld_control%eps_cutoff)
903  END DO
904  END IF
905 
906  ! Print atomic densities if requested
907  IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
908  DO i = 1, cdft_control%natoms
909  CALL transfer_rs2pw(rs_single(i), cdft_control%group(igroup)%hw_rho_atomic(i))
910  END DO
911  CALL cdft_print_hirshfeld_density(qs_env)
912  END IF
913 
914  END DO
915 
916  DO igroup = 1, SIZE(cdft_control%group)
917 
918  CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
919 
920  IF (.NOT. cdft_control%in_memory .AND. igroup == 1) THEN
921  CALL auxbas_pw_pool%give_back_pw(cdft_control%group(1)%hw_rho_total)
922  END IF
923 
924  IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
925  DO i = 1, cdft_control%natoms
926  CALL rs_grid_release(rs_single(i))
927  CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic(i))
928  END DO
929  DEALLOCATE (rs_single)
930  DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic)
931  END IF
932 
933  IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
934  DO i = 1, cdft_control%natoms
935  CALL rs_grid_release(rs_single_charge(i))
936  CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(i))
937  END DO
938  DEALLOCATE (rs_single_charge)
939  DEALLOCATE (compute_charge)
940  DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge)
941  END IF
942 
943  END DO
944 
945  IF (cdft_control%in_memory) THEN
946  DO igroup = 1, SIZE(cdft_control%group)
947  ALLOCATE (cdft_control%group(igroup)%gradients_x(1*natom, lb_pw(1):ub_pw(1), &
948  lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
949  cdft_control%group(igroup)%gradients_x(:, :, :, :) = 0.0_dp
950  END DO
951  END IF
952 
953  IF (cdft_control%in_memory) THEN
954  DO igroup = 1, SIZE(cdft_control%group)
955 
956  ALLOCATE (pab(1, 1))
957  nthread = 1
958  ithread = 0
959  atoms_memory = hirshfeld_control%atoms_memory
960 
961  DO ikind = 1, SIZE(atomic_kind_set)
962  numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
963  IF (numexp <= 0) cycle
964  CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
965 
966  ALLOCATE (pw_single_dr(num_species))
967  ALLOCATE (rs_single_dr(num_species))
968 
969  DO i = 1, num_species
970  CALL auxbas_pw_pool%create_pw(pw_single_dr(i))
971  CALL pw_zero(pw_single_dr(i))
972  END DO
973 
974  atoms_memory_num = SIZE((/(j, j=1, num_species, atoms_memory)/))
975 
976  ! Can't store all pw grids, therefore split into groups of size atom_memory
977  ! Ideally this code should be re-written to be more memory efficient
978  IF (num_species > atoms_memory) THEN
979  ALLOCATE (num_species_small(atoms_memory_num + 1))
980  num_species_small(1:atoms_memory_num) = (/(j, j=1, num_species, atoms_memory)/)
981  num_species_small(atoms_memory_num + 1) = num_species
982  ELSE
983  ALLOCATE (num_species_small(2))
984  num_species_small(:) = (/1, num_species/)
985  END IF
986 
987  DO k = 1, SIZE(num_species_small) - 1
988  IF (num_species > atoms_memory) THEN
989  ALLOCATE (cores(num_species_small(k + 1) - (num_species_small(k) - 1)))
990  ELSE
991  ALLOCATE (cores(num_species))
992  END IF
993 
994  DO i = num_species_small(k), num_species_small(k + 1)
995  CALL rs_grid_create(rs_single_dr(i), auxbas_rs_desc)
996  CALL rs_grid_zero(rs_single_dr(i))
997  END DO
998  DO iex = 1, numexp
999 
1000  alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
1001  coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
1002  prefactor = 2.0_dp*alpha
1003  npme = 0
1004  cores = 0
1005 
1006  DO iatom = 1, SIZE(cores)
1007  atom_a = atom_list(iatom + (num_species_small(k) - 1))
1008  ra(:) = pbc(particle_set(atom_a)%r, cell)
1009 
1010  IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
1011  IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
1012  npme = npme + 1
1013  cores(npme) = iatom
1014  END IF
1015  ELSE
1016  npme = npme + 1
1017  cores(npme) = iatom
1018  END IF
1019  END DO
1020  DO j = 1, npme
1021  iatom = cores(j)
1022  atom_a = atom_list(iatom + (num_species_small(k) - 1))
1023  pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
1024  ra(:) = pbc(particle_set(atom_a)%r, cell)
1025  subpatch_pattern = 0
1026 
1027  ! Calculate cutoff
1028  IF (hirshfeld_control%use_atomic_cutoff) THEN
1029  radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
1030  ra=ra, rb=ra, rp=ra, &
1031  zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
1032  pab=pab, o1=0, o2=0, & ! without map_consistent
1033  prefactor=1.0_dp, cutoff=0.0_dp)
1034  END IF
1035 
1036  CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1037  (/0.0_dp, 0.0_dp, 0.0_dp/), prefactor, &
1038  pab, 0, 0, rs_single_dr(iatom + (num_species_small(k) - 1)), &
1039  radius=radius, &
1040  ga_gb_function=grid_func_ab, use_subpatch=.true., &
1041  subpatch_pattern=subpatch_pattern)
1042 
1043  END DO
1044  END DO
1045 
1046  DO iatom = num_species_small(k), num_species_small(k + 1)
1047  CALL transfer_rs2pw(rs_single_dr(iatom), pw_single_dr(iatom))
1048  CALL rs_grid_release(rs_single_dr(iatom))
1049  END DO
1050 
1051  DEALLOCATE (cores)
1052  END DO
1053 
1054  DO iatom = 1, num_species
1055  atom_a = atom_list(iatom)
1056  cdft_control%group(igroup)%gradients_x(atom_a, :, :, :) = pw_single_dr(iatom)%array(:, :, :)
1057  CALL auxbas_pw_pool%give_back_pw(pw_single_dr(iatom))
1058  END DO
1059 
1060  DEALLOCATE (rs_single_dr)
1061  DEALLOCATE (num_species_small)
1062  DEALLOCATE (pw_single_dr)
1063  END DO
1064  DEALLOCATE (pab)
1065  END DO
1066  END IF
1067 
1068  IF (cdft_control%in_memory) THEN
1069  DO igroup = 1, SIZE(cdft_control%group)
1070  ALLOCATE (cdft_control%group(igroup)%gradients_y(1*num_atoms, lb_pw(1):ub_pw(1), &
1071  lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1072  ALLOCATE (cdft_control%group(igroup)%gradients_z(1*num_atoms, lb_pw(1):ub_pw(1), &
1073  lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1074  cdft_control%group(igroup)%gradients_y(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1075  cdft_control%group(igroup)%gradients_z(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1076  END DO
1077  END IF
1078 
1079  ! Calculate gradient if requested
1080  IF (cdft_control%in_memory) THEN
1081 
1082  DO igroup = 1, SIZE(cdft_control%group)
1083  DO k = lb_pw(3), ub_pw(3)
1084  DO j = lb_pw(2), ub_pw(2)
1085  DO i = lb_pw(1), ub_pw(1)
1086  DO iatom = 1, natom
1087 
1088  ra(:) = particle_set(iatom)%r
1089 
1090  IF (cdft_control%group(igroup)%hw_rho_total%array(i, j, k) > hirshfeld_control%eps_cutoff) THEN
1091 
1092  exp_eval = (coefficients(iatom) - &
1093  cdft_control%group(igroup)%weight%array(i, j, k))/ &
1094  cdft_control%group(igroup)%hw_rho_total%array(i, j, k)
1095 
1096  r2 = (/i*dr_pw(1), j*dr_pw(2), k*dr_pw(3)/) + origin
1097  r_pbc = pbc(ra, r2, cell)
1098 
1099  ! Store gradient d/dR_x w, including term: (r_x - R_x)
1100  cdft_control%group(igroup)%gradients_x(iatom, i, j, k) = &
1101  cdft_control%group(igroup)%gradients_x(iatom, i, j, k)* &
1102  r_pbc(1)*exp_eval
1103 
1104  ! Store gradient d/dR_y w, including term: (r_y - R_y)
1105  cdft_control%group(igroup)%gradients_y(iatom, i, j, k) = &
1106  cdft_control%group(igroup)%gradients_y(iatom, i, j, k)* &
1107  r_pbc(2)*exp_eval
1108 
1109  ! Store gradient d/dR_z w, including term:(r_z - R_z)
1110  cdft_control%group(igroup)%gradients_z(iatom, i, j, k) = &
1111  cdft_control%group(igroup)%gradients_z(iatom, i, j, k)* &
1112  r_pbc(3)*exp_eval
1113 
1114  END IF
1115  END DO
1116  END DO
1117  END DO
1118  END DO
1119 
1120  IF (igroup == 1) THEN
1121  CALL auxbas_pw_pool%give_back_pw(cdft_control%group(1)%hw_rho_total)
1122  END IF
1123  END DO
1124  END IF
1125 
1126  CALL rs_grid_release(rs_rho_all)
1127 
1128  IF (ALLOCATED(coefficients)) DEALLOCATE (coefficients)
1129  IF (ALLOCATED(is_constraint)) DEALLOCATE (is_constraint)
1130 
1131  CALL timestop(handle)
1132 
1133  END SUBROUTINE hirshfeld_constraint_low
1134 
1135 ! **************************************************************************************************
1136 !> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT
1137 !> weight function and the realspace electron density
1138 !> \param qs_env ...
1139 ! **************************************************************************************************
1140  SUBROUTINE cdft_constraint_integrate(qs_env)
1141  TYPE(qs_environment_type), POINTER :: qs_env
1142 
1143  CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_integrate'
1144 
1145  INTEGER :: handle, i, iatom, igroup, ikind, ivar, &
1146  iw, jatom, natom, nvar
1147  LOGICAL :: is_becke, paw_atom
1148  REAL(kind=dp) :: dvol, eps_cavity, sign
1149  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: de, strength, target_val
1150  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge, gapw_offset
1151  TYPE(becke_constraint_type), POINTER :: becke_control
1152  TYPE(cdft_control_type), POINTER :: cdft_control
1153  TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1154  TYPE(cp_logger_type), POINTER :: logger
1155  TYPE(dft_control_type), POINTER :: dft_control
1156  TYPE(mp_para_env_type), POINTER :: para_env
1157  TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
1158  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1159  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge, rho_r
1160  TYPE(qs_energy_type), POINTER :: energy
1161  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1162  TYPE(qs_rho_type), POINTER :: rho
1163  TYPE(rho0_mpole_type), POINTER :: rho0_mpole
1164  TYPE(section_vals_type), POINTER :: cdft_constraint_section
1165 
1166  NULLIFY (para_env, dft_control, particle_set, rho_r, energy, rho, &
1167  logger, cdft_constraint_section, qs_kind_set, mp_rho, &
1168  rho0_mpole, group, charge)
1169  CALL timeset(routinen, handle)
1170  logger => cp_get_default_logger()
1171  CALL get_qs_env(qs_env, &
1172  particle_set=particle_set, &
1173  rho=rho, &
1174  natom=natom, &
1175  dft_control=dft_control, &
1176  para_env=para_env, &
1177  qs_kind_set=qs_kind_set)
1178  CALL qs_rho_get(rho, rho_r=rho_r)
1179  cpassert(ASSOCIATED(qs_kind_set))
1180  cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1181  iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
1182  cdft_control => dft_control%qs_control%cdft_control
1183  is_becke = (cdft_control%type == outer_scf_becke_constraint)
1184  becke_control => cdft_control%becke_control
1185  IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
1186  cpabort("Becke control has not been allocated.")
1187  group => cdft_control%group
1188  ! Initialize
1189  nvar = SIZE(cdft_control%target)
1190  ALLOCATE (strength(nvar))
1191  ALLOCATE (target_val(nvar))
1192  ALLOCATE (de(nvar))
1193  strength(:) = cdft_control%strength(:)
1194  target_val(:) = cdft_control%target(:)
1195  sign = 1.0_dp
1196  de = 0.0_dp
1197  dvol = group(1)%weight%pw_grid%dvol
1198  IF (cdft_control%atomic_charges) THEN
1199  charge => cdft_control%charge
1200  ALLOCATE (electronic_charge(cdft_control%natoms, dft_control%nspins))
1201  electronic_charge = 0.0_dp
1202  END IF
1203  ! Calculate value of constraint i.e. int ( rho(r) w(r) dr)
1204  DO i = 1, dft_control%nspins
1205  DO igroup = 1, SIZE(group)
1206  SELECT CASE (group(igroup)%constraint_type)
1207  CASE (cdft_charge_constraint)
1208  sign = 1.0_dp
1209  CASE (cdft_magnetization_constraint)
1210  IF (i == 1) THEN
1211  sign = 1.0_dp
1212  ELSE
1213  sign = -1.0_dp
1214  END IF
1215  CASE (cdft_alpha_constraint)
1216  sign = 1.0_dp
1217  IF (i == 2) cycle
1218  CASE (cdft_beta_constraint)
1219  sign = 1.0_dp
1220  IF (i == 1) cycle
1221  CASE DEFAULT
1222  cpabort("Unknown constraint type.")
1223  END SELECT
1224  IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1225  ! With external control, we can use cavity_mat as a mask to kahan sum
1226  eps_cavity = becke_control%eps_cavity
1227  IF (igroup /= 1) &
1228  CALL cp_abort(__location__, &
1229  "Multiple constraints not yet supported by parallel mixed calculations.")
1230  de(igroup) = de(igroup) + sign*accurate_dot_product(group(igroup)%weight%array, rho_r(i)%array, &
1231  becke_control%cavity_mat, eps_cavity)*dvol
1232  ELSE
1233  de(igroup) = de(igroup) + sign*pw_integral_ab(group(igroup)%weight, rho_r(i), local_only=.true.)
1234  END IF
1235  END DO
1236  IF (cdft_control%atomic_charges) THEN
1237  DO iatom = 1, cdft_control%natoms
1238  electronic_charge(iatom, i) = pw_integral_ab(charge(iatom), rho_r(i), local_only=.true.)
1239  END DO
1240  END IF
1241  END DO
1242  CALL get_qs_env(qs_env, energy=energy)
1243  CALL para_env%sum(de)
1244  IF (cdft_control%atomic_charges) THEN
1245  CALL para_env%sum(electronic_charge)
1246  END IF
1247  ! Use fragment densities as reference value (= Becke deformation density)
1248  IF (cdft_control%fragment_density .AND. .NOT. cdft_control%fragments_integrated) THEN
1249  CALL prepare_fragment_constraint(qs_env)
1250  END IF
1251  IF (dft_control%qs_control%gapw) THEN
1252  ! GAPW: add core charges (rho_hard - rho_soft)
1253  IF (cdft_control%fragment_density) &
1254  CALL cp_abort(__location__, &
1255  "Fragment constraints not yet compatible with GAPW.")
1256  ALLOCATE (gapw_offset(nvar, dft_control%nspins))
1257  gapw_offset = 0.0_dp
1258  CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole)
1259  CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho)
1260  DO i = 1, dft_control%nspins
1261  DO igroup = 1, SIZE(group)
1262  DO iatom = 1, SIZE(group(igroup)%atoms)
1263  SELECT CASE (group(igroup)%constraint_type)
1264  CASE (cdft_charge_constraint)
1265  sign = 1.0_dp
1266  CASE (cdft_magnetization_constraint)
1267  IF (i == 1) THEN
1268  sign = 1.0_dp
1269  ELSE
1270  sign = -1.0_dp
1271  END IF
1272  CASE (cdft_alpha_constraint)
1273  sign = 1.0_dp
1274  IF (i == 2) cycle
1275  CASE (cdft_beta_constraint)
1276  sign = 1.0_dp
1277  IF (i == 1) cycle
1278  CASE DEFAULT
1279  cpabort("Unknown constraint type.")
1280  END SELECT
1281  jatom = group(igroup)%atoms(iatom)
1282  CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
1283  CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
1284  IF (paw_atom) THEN
1285  gapw_offset(igroup, i) = gapw_offset(igroup, i) + sign*group(igroup)%coeff(iatom)*mp_rho(jatom)%q0(i)
1286  END IF
1287  END DO
1288  END DO
1289  END DO
1290  IF (cdft_control%atomic_charges) THEN
1291  DO iatom = 1, cdft_control%natoms
1292  jatom = cdft_control%atoms(iatom)
1293  CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
1294  CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
1295  IF (paw_atom) THEN
1296  DO i = 1, dft_control%nspins
1297  electronic_charge(iatom, i) = electronic_charge(iatom, i) + mp_rho(jatom)%q0(i)
1298  END DO
1299  END IF
1300  END DO
1301  END IF
1302  DO i = 1, dft_control%nspins
1303  DO ivar = 1, nvar
1304  de(ivar) = de(ivar) + gapw_offset(ivar, i)
1305  END DO
1306  END DO
1307  DEALLOCATE (gapw_offset)
1308  END IF
1309  ! Update constraint value and energy
1310  cdft_control%value(:) = de(:)
1311  energy%cdft = 0.0_dp
1312  DO ivar = 1, nvar
1313  energy%cdft = energy%cdft + (de(ivar) - target_val(ivar))*strength(ivar)
1314  END DO
1315  ! Print constraint info and atomic CDFT charges
1316  CALL cdft_constraint_print(qs_env, electronic_charge)
1317  ! Deallocate tmp storage
1318  DEALLOCATE (de, strength, target_val)
1319  IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge)
1320  CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO")
1321  CALL timestop(handle)
1322 
1323  END SUBROUTINE cdft_constraint_integrate
1324 
1325 ! **************************************************************************************************
1326 !> \brief Calculates atomic forces due to a CDFT constraint (Becke or Hirshfeld)
1327 !> \param qs_env ...
1328 ! **************************************************************************************************
1329  SUBROUTINE cdft_constraint_force(qs_env)
1330  TYPE(qs_environment_type), POINTER :: qs_env
1331 
1332  CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_force'
1333 
1334  INTEGER :: handle, i, iatom, igroup, ikind, ispin, &
1335  j, k, natom, nvar
1336  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
1337  INTEGER, DIMENSION(2, 3) :: bo
1338  INTEGER, DIMENSION(3) :: lb, ub
1339  REAL(kind=dp) :: dvol, eps_cavity, sign
1340  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength
1341  REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
1342  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1343  TYPE(becke_constraint_type), POINTER :: becke_control
1344  TYPE(cdft_control_type), POINTER :: cdft_control
1345  TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1346  TYPE(cell_type), POINTER :: cell
1347  TYPE(dft_control_type), POINTER :: dft_control
1348  TYPE(mp_para_env_type), POINTER :: para_env
1349  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1350  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1351  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1352  TYPE(qs_rho_type), POINTER :: rho
1353 
1354  CALL timeset(routinen, handle)
1355  NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, &
1356  rho, rho_r, force, cutoffs, becke_control, group)
1357 
1358  CALL get_qs_env(qs_env, &
1359  atomic_kind_set=atomic_kind_set, &
1360  natom=natom, &
1361  particle_set=particle_set, &
1362  cell=cell, &
1363  rho=rho, &
1364  force=force, &
1365  dft_control=dft_control, &
1366  para_env=para_env)
1367  CALL qs_rho_get(rho, rho_r=rho_r)
1368 
1369  cdft_control => dft_control%qs_control%cdft_control
1370  becke_control => cdft_control%becke_control
1371  group => cdft_control%group
1372  nvar = SIZE(cdft_control%target)
1373  ALLOCATE (strength(nvar))
1374  strength(:) = cdft_control%strength(:)
1375  cutoffs => cdft_control%becke_control%cutoffs
1376  eps_cavity = cdft_control%becke_control%eps_cavity
1377 
1378  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1379  atom_of_kind=atom_of_kind, &
1380  kind_of=kind_of)
1381  DO igroup = 1, SIZE(cdft_control%group)
1382  ALLOCATE (cdft_control%group(igroup)%integrated(3, natom))
1383  cdft_control%group(igroup)%integrated = 0.0_dp
1384  END DO
1385 
1386  lb(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
1387  ub(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
1388  bo = cdft_control%group(1)%weight%pw_grid%bounds_local
1389  dvol = cdft_control%group(1)%weight%pw_grid%dvol
1390  sign = 1.0_dp
1391 
1392  IF (cdft_control%type == outer_scf_becke_constraint) THEN
1393  IF (.NOT. cdft_control%becke_control%in_memory) THEN
1394  CALL becke_constraint_low(qs_env, just_gradients=.true.)
1395  END IF
1396 
1397  ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1398  IF (.NOT. cdft_control%in_memory) THEN
1399  CALL hirshfeld_constraint_low(qs_env, just_gradients=.true.)
1400  END IF
1401  END IF
1402 
1403  ! If no Becke Gaussian confinement
1404  IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN
1405  ! No external control
1406  DO k = bo(1, 1), bo(2, 1)
1407  DO j = bo(1, 2), bo(2, 2)
1408  DO i = bo(1, 3), bo(2, 3)
1409  ! First check if this grid point should be skipped
1410  IF (cdft_control%becke_control%cavity_confine) THEN
1411  IF (cdft_control%becke_control%cavity%array(k, j, i) < eps_cavity) cycle
1412  END IF
1413 
1414  DO igroup = 1, SIZE(cdft_control%group)
1415  DO iatom = 1, natom
1416  DO ispin = 1, dft_control%nspins
1417 
1418  SELECT CASE (cdft_control%group(igroup)%constraint_type)
1419  CASE (cdft_charge_constraint)
1420  sign = 1.0_dp
1421  CASE (cdft_magnetization_constraint)
1422  IF (ispin == 1) THEN
1423  sign = 1.0_dp
1424  ELSE
1425  sign = -1.0_dp
1426  END IF
1427  CASE (cdft_alpha_constraint)
1428  sign = 1.0_dp
1429  IF (ispin == 2) cycle
1430  CASE (cdft_beta_constraint)
1431  sign = 1.0_dp
1432  IF (ispin == 1) cycle
1433  CASE DEFAULT
1434  cpabort("Unknown constraint type.")
1435  END SELECT
1436 
1437  IF (cdft_control%type == outer_scf_becke_constraint) THEN
1438 
1439  cdft_control%group(igroup)%integrated(:, iatom) = &
1440  cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1441  cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1442  *rho_r(ispin)%array(k, j, i) &
1443  *dvol
1444 
1445  ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1446 
1447  cdft_control%group(igroup)%integrated(1, iatom) = &
1448  cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1449  cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1450  *rho_r(ispin)%array(k, j, i) &
1451  *dvol
1452 
1453  cdft_control%group(igroup)%integrated(2, iatom) = &
1454  cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1455  cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1456  *rho_r(ispin)%array(k, j, i) &
1457  *dvol
1458 
1459  cdft_control%group(igroup)%integrated(3, iatom) = &
1460  cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1461  cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1462  *rho_r(ispin)%array(k, j, i) &
1463  *dvol
1464 
1465  END IF
1466 
1467  END DO
1468  END DO
1469  END DO
1470  END DO
1471  END DO
1472  END DO
1473 
1474  ! If Becke Gaussian confinement
1475  ELSE
1476  DO k = lbound(cdft_control%becke_control%cavity_mat, 1), ubound(cdft_control%becke_control%cavity_mat, 1)
1477  DO j = lbound(cdft_control%becke_control%cavity_mat, 2), ubound(cdft_control%becke_control%cavity_mat, 2)
1478  DO i = lbound(cdft_control%becke_control%cavity_mat, 3), ubound(cdft_control%becke_control%cavity_mat, 3)
1479 
1480  ! First check if this grid point should be skipped
1481  IF (cdft_control%becke_control%cavity_mat(k, j, i) < eps_cavity) cycle
1482 
1483  DO igroup = 1, SIZE(group)
1484  DO iatom = 1, natom
1485  DO ispin = 1, dft_control%nspins
1486  SELECT CASE (group(igroup)%constraint_type)
1487  CASE (cdft_charge_constraint)
1488  sign = 1.0_dp
1489  CASE (cdft_magnetization_constraint)
1490  IF (ispin == 1) THEN
1491  sign = 1.0_dp
1492  ELSE
1493  sign = -1.0_dp
1494  END IF
1495  CASE (cdft_alpha_constraint)
1496  sign = 1.0_dp
1497  IF (ispin == 2) cycle
1498  CASE (cdft_beta_constraint)
1499  sign = 1.0_dp
1500  IF (ispin == 1) cycle
1501  CASE DEFAULT
1502  cpabort("Unknown constraint type.")
1503  END SELECT
1504 
1505  ! Integrate gradient of weight function
1506  IF (cdft_control%type == outer_scf_becke_constraint) THEN
1507 
1508  cdft_control%group(igroup)%integrated(:, iatom) = &
1509  cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1510  cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1511  *rho_r(ispin)%array(k, j, i) &
1512  *dvol
1513 
1514  ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1515 
1516  cdft_control%group(igroup)%integrated(1, iatom) = &
1517  cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1518  cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1519  *rho_r(ispin)%array(k, j, i) &
1520  *dvol
1521 
1522  cdft_control%group(igroup)%integrated(2, iatom) = &
1523  cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1524  cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1525  *rho_r(ispin)%array(k, j, i) &
1526  *dvol
1527 
1528  cdft_control%group(igroup)%integrated(3, iatom) = &
1529  cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1530  cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1531  *rho_r(ispin)%array(k, j, i) &
1532  *dvol
1533 
1534  END IF
1535 
1536  END DO
1537  END DO
1538  END DO
1539  END DO
1540  END DO
1541  END DO
1542  END IF
1543 
1544  IF (.NOT. cdft_control%transfer_pot) THEN
1545  IF (cdft_control%type == outer_scf_becke_constraint) THEN
1546  DO igroup = 1, SIZE(group)
1547  DEALLOCATE (cdft_control%group(igroup)%gradients)
1548  END DO
1549  ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1550  DO igroup = 1, SIZE(group)
1551  DEALLOCATE (cdft_control%group(igroup)%gradients_x)
1552  DEALLOCATE (cdft_control%group(igroup)%gradients_y)
1553  DEALLOCATE (cdft_control%group(igroup)%gradients_z)
1554  END DO
1555  END IF
1556  END IF
1557 
1558  DO igroup = 1, SIZE(group)
1559  CALL para_env%sum(group(igroup)%integrated)
1560  END DO
1561 
1562  ! Update force only on master process. Otherwise force due to constraint becomes multiplied
1563  ! by the number of processes when the final force%rho_elec is constructed in qs_force
1564  ! by mp_summing [the final integrated(:,:) is distributed on all processors]
1565  IF (para_env%is_source()) THEN
1566  DO igroup = 1, SIZE(group)
1567  DO iatom = 1, natom
1568  ikind = kind_of(iatom)
1569  i = atom_of_kind(iatom)
1570  force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup)
1571  END DO
1572  END DO
1573  END IF
1574 
1575  DEALLOCATE (strength)
1576  DO igroup = 1, SIZE(group)
1577  DEALLOCATE (group(igroup)%integrated)
1578  END DO
1579  NULLIFY (group)
1580 
1581  CALL timestop(handle)
1582 
1583  END SUBROUTINE cdft_constraint_force
1584 
1585 ! **************************************************************************************************
1586 !> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied
1587 !> by the CDFT weight functions and integrated over the realspace grid.
1588 !> \param qs_env ...
1589 ! **************************************************************************************************
1590  SUBROUTINE prepare_fragment_constraint(qs_env)
1591  TYPE(qs_environment_type), POINTER :: qs_env
1592 
1593  CHARACTER(len=*), PARAMETER :: routinen = 'prepare_fragment_constraint'
1594 
1595  INTEGER :: handle, i, iatom, igroup, natom, &
1596  nelectron_total, nfrag_spins
1597  LOGICAL :: is_becke, needs_spin_density
1598  REAL(kind=dp) :: dvol, multiplier(2), nelectron_frag
1599  TYPE(becke_constraint_type), POINTER :: becke_control
1600  TYPE(cdft_control_type), POINTER :: cdft_control
1601  TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1602  TYPE(cp_logger_type), POINTER :: logger
1603  TYPE(dft_control_type), POINTER :: dft_control
1604  TYPE(mp_para_env_type), POINTER :: para_env
1605  TYPE(pw_env_type), POINTER :: pw_env
1606  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1607  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: rho_frag
1608  TYPE(qs_subsys_type), POINTER :: subsys
1609 
1610  NULLIFY (para_env, dft_control, logger, subsys, pw_env, auxbas_pw_pool, group)
1611  CALL timeset(routinen, handle)
1612  logger => cp_get_default_logger()
1613  CALL get_qs_env(qs_env, &
1614  natom=natom, &
1615  dft_control=dft_control, &
1616  para_env=para_env)
1617 
1618  cdft_control => dft_control%qs_control%cdft_control
1619  is_becke = (cdft_control%type == outer_scf_becke_constraint)
1620  becke_control => cdft_control%becke_control
1621  IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
1622  cpabort("Becke control has not been allocated.")
1623  group => cdft_control%group
1624  dvol = group(1)%weight%pw_grid%dvol
1625  ! Fragment densities are meaningful only for some calculation types
1626  IF (.NOT. qs_env%single_point_run) &
1627  CALL cp_abort(__location__, &
1628  "CDFT fragment constraints are only compatible with single "// &
1629  "point calculations (run_type ENERGY or ENERGY_FORCE).")
1630  IF (dft_control%qs_control%gapw) &
1631  CALL cp_abort(__location__, &
1632  "CDFT fragment constraint not compatible with GAPW.")
1633  needs_spin_density = .false.
1634  multiplier = 1.0_dp
1635  nfrag_spins = 1
1636  DO igroup = 1, SIZE(group)
1637  SELECT CASE (group(igroup)%constraint_type)
1638  CASE (cdft_charge_constraint)
1639  ! Do nothing
1640  CASE (cdft_magnetization_constraint)
1641  needs_spin_density = .true.
1642  CASE (cdft_alpha_constraint, cdft_beta_constraint)
1643  CALL cp_abort(__location__, &
1644  "CDFT fragment constraint not yet compatible with "// &
1645  "spin specific constraints.")
1646  CASE DEFAULT
1647  cpabort("Unknown constraint type.")
1648  END SELECT
1649  END DO
1650  IF (needs_spin_density) THEN
1651  nfrag_spins = 2
1652  DO i = 1, 2
1653  IF (cdft_control%flip_fragment(i)) multiplier(i) = -1.0_dp
1654  END DO
1655  END IF
1656  ! Read fragment reference densities
1657  ALLOCATE (cdft_control%fragments(nfrag_spins, 2))
1658  ALLOCATE (rho_frag(nfrag_spins))
1659  CALL get_qs_env(qs_env, pw_env=pw_env)
1660  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1661  ! Total density (rho_alpha + rho_beta)
1662  CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 1))
1663  CALL cp_cube_to_pw(cdft_control%fragments(1, 1), &
1664  cdft_control%fragment_a_fname, 1.0_dp)
1665  CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 2))
1666  CALL cp_cube_to_pw(cdft_control%fragments(1, 2), &
1667  cdft_control%fragment_b_fname, 1.0_dp)
1668  ! Spin difference density (rho_alpha - rho_beta) if needed
1669  IF (needs_spin_density) THEN
1670  CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 1))
1671  CALL cp_cube_to_pw(cdft_control%fragments(2, 1), &
1672  cdft_control%fragment_a_spin_fname, multiplier(1))
1673  CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 2))
1674  CALL cp_cube_to_pw(cdft_control%fragments(2, 2), &
1675  cdft_control%fragment_b_spin_fname, multiplier(2))
1676  END IF
1677  ! Sum up fragments
1678  DO i = 1, nfrag_spins
1679  CALL auxbas_pw_pool%create_pw(rho_frag(i))
1680  CALL pw_copy(cdft_control%fragments(i, 1), rho_frag(i))
1681  CALL pw_axpy(cdft_control%fragments(i, 2), rho_frag(i), 1.0_dp)
1682  CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 1))
1683  CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 2))
1684  END DO
1685  DEALLOCATE (cdft_control%fragments)
1686  ! Check that the number of electrons is consistent
1687  CALL get_qs_env(qs_env, subsys=subsys)
1688  CALL qs_subsys_get(subsys, nelectron_total=nelectron_total)
1689  nelectron_frag = pw_integrate_function(rho_frag(1))
1690  IF (nint(nelectron_frag) /= nelectron_total) &
1691  CALL cp_abort(__location__, &
1692  "The number of electrons in the reference and interacting "// &
1693  "configurations does not match. Check your fragment cube files.")
1694  ! Update constraint target value i.e. perform integration w_i*rho_frag_{tot/spin}*dr
1695  cdft_control%target = 0.0_dp
1696  DO igroup = 1, SIZE(group)
1697  IF (group(igroup)%constraint_type == cdft_charge_constraint) THEN
1698  i = 1
1699  ELSE
1700  i = 2
1701  END IF
1702  IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1703  cdft_control%target(igroup) = cdft_control%target(igroup) + &
1704  accurate_dot_product(group(igroup)%weight%array, rho_frag(i)%array, &
1705  becke_control%cavity_mat, becke_control%eps_cavity)*dvol
1706  ELSE
1707  cdft_control%target(igroup) = cdft_control%target(igroup) + &
1708  pw_integral_ab(group(igroup)%weight, rho_frag(i), local_only=.true.)
1709  END IF
1710  END DO
1711  CALL para_env%sum(cdft_control%target)
1712  ! Calculate reference atomic charges int( w_i * rho_frag * dr )
1713  IF (cdft_control%atomic_charges) THEN
1714  ALLOCATE (cdft_control%charges_fragment(cdft_control%natoms, nfrag_spins))
1715  DO i = 1, nfrag_spins
1716  DO iatom = 1, cdft_control%natoms
1717  cdft_control%charges_fragment(iatom, i) = &
1718  pw_integral_ab(cdft_control%charge(iatom), rho_frag(i), local_only=.true.)
1719  END DO
1720  END DO
1721  CALL para_env%sum(cdft_control%charges_fragment)
1722  END IF
1723  DO i = 1, nfrag_spins
1724  CALL auxbas_pw_pool%give_back_pw(rho_frag(i))
1725  END DO
1726  DEALLOCATE (rho_frag)
1727  cdft_control%fragments_integrated = .true.
1728 
1729  CALL timestop(handle)
1730 
1731  END SUBROUTINE prepare_fragment_constraint
1732 
1733 END MODULE qs_cdft_methods
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
ao_util
All kind of helpful little routines.
Definition: ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition: ao_util.F:209
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
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_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
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition: cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_cube_to_pw
subroutine, public cp_cube_to_pw(grid, filename, scaling, silent)
Thin wrapper around routine cube_to_pw.
Definition: cp_realspace_grid_cube.F:113
grid_api
Fortran API for the grid package, which is written in C.
Definition: grid_api.F:12
grid_api::grid_func_ab
integer, parameter, public grid_func_ab
Definition: grid_api.F:29
grid_api::collocate_pgf_product
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
Definition: grid_api.F:119
hirshfeld_types
The types needed for the calculation of Hirshfeld charges and related functions.
Definition: hirshfeld_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition: input_constants.F:766
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition: input_constants.F:766
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition: input_constants.F:766
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition: input_constants.F:718
input_constants::outer_scf_hirshfeld_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Definition: input_constants.F:718
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition: input_constants.F:766
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
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition: kahan_sum.F:29
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
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_cdft_methods
Subroutines for building CDFT constraints.
Definition: qs_cdft_methods.F:14
qs_cdft_methods::becke_constraint
subroutine, public becke_constraint(qs_env, calc_pot, calculate_forces)
Driver routine for calculating a Becke constraint.
Definition: qs_cdft_methods.F:105
qs_cdft_methods::hirshfeld_constraint
subroutine, public hirshfeld_constraint(qs_env, calc_pot, calculate_forces)
Driver routine for calculating a Hirshfeld constraint.
Definition: qs_cdft_methods.F:588
qs_cdft_types
Defines CDFT control structures.
Definition: qs_cdft_types.F:14
qs_cdft_utils
Utility subroutines for CDFT calculations.
Definition: qs_cdft_utils.F:14
qs_cdft_utils::cdft_constraint_print
subroutine, public cdft_constraint_print(qs_env, electronic_charge)
Prints information about CDFT constraints.
Definition: qs_cdft_utils.F:1103
qs_cdft_utils::hirshfeld_constraint_init
subroutine, public hirshfeld_constraint_init(qs_env)
Initializes Gaussian Hirshfeld constraints.
Definition: qs_cdft_utils.F:974
qs_cdft_utils::becke_constraint_init
subroutine, public becke_constraint_init(qs_env)
Initializes the Becke constraint environment.
Definition: qs_cdft_utils.F:97
qs_cdft_utils::hfun_scale
subroutine, public hfun_scale(fout, fun1, fun2, divide, small)
Calculate fout = fun1/fun2 or fout = fun1*fun2.
Definition: qs_cdft_utils.F:856
qs_cdft_utils::cdft_print_hirshfeld_density
subroutine, public cdft_print_hirshfeld_density(qs_env)
Prints Hirshfeld weight function and promolecule density.
Definition: qs_cdft_utils.F:1313
qs_energy_types
Definition: qs_energy_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_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_rho0_types
Definition: qs_rho0_types.F:8
qs_rho0_types::get_rho0_mpole
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, Qlm_gg, Qlm_car, Qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs)
...
Definition: qs_rho0_types.F:425
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
qs_subsys_types
types that represent a quickstep subsys
Definition: qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition: qs_subsys_types.F:134
realspace_grid_types
Definition: realspace_grid_types.F:18
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition: realspace_grid_types.F:482
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition: realspace_grid_types.F:656
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition: realspace_grid_types.F:2016
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition: realspace_grid_types.F:1909