(git:6344619)
Loading...
Searching...
No Matches
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-2026 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! **************************************************************************************************
14MODULE 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_files, ONLY: close_file,&
23 open_file
24 USE cp_log_handling, ONLY: cp_get_default_logger,&
25 cp_logger_type
26 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
27 cp_print_key_unit_nr
28 USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw
29 USE grid_api, ONLY: grid_func_ab,&
30 collocate_pgf_product
31 USE hirshfeld_types, ONLY: hirshfeld_type
32 USE input_constants, ONLY: cdft_alpha_constraint,&
33 cdft_beta_constraint,&
34 cdft_charge_constraint,&
35 cdft_magnetization_constraint,&
36 outer_scf_becke_constraint,&
37 outer_scf_hirshfeld_constraint
38 USE input_section_types, ONLY: section_vals_get_subs_vals,&
39 section_vals_type
40 USE kahan_sum, ONLY: accurate_dot_product
41 USE kinds, ONLY: dp
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_integral_ab,&
49 pw_integrate_function,&
50 pw_set,&
51 pw_zero
52 USE pw_pool_types, ONLY: pw_pool_type
53 USE pw_types, ONLY: pw_r3d_rs_type
54 USE qs_cdft_types, ONLY: becke_constraint_type,&
55 cdft_control_type,&
56 cdft_group_type,&
57 hirshfeld_constraint_type
58 USE qs_cdft_utils, ONLY: becke_constraint_init,&
59 cdft_constraint_print,&
60 cdft_print_hirshfeld_density,&
61 hfun_scale,&
62 hirshfeld_constraint_init
63 USE qs_energy_types, ONLY: qs_energy_type
64 USE qs_environment_types, ONLY: get_qs_env,&
65 qs_environment_type
66 USE qs_force_types, ONLY: qs_force_type
67 USE qs_rho_types, ONLY: qs_rho_get,&
68 qs_rho_type
69 USE qs_subsys_types, ONLY: qs_subsys_get,&
70 qs_subsys_type
71 USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
72 realspace_grid_type,&
73 rs_grid_create,&
74 rs_grid_release,&
75 rs_grid_zero,&
76 transfer_rs2pw
77#include "./base/base_uses.f90"
78
79 IMPLICIT NONE
80
81 PRIVATE
82
83 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_methods'
84 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
85
86! *** Public subroutines ***
87
88 PUBLIC :: becke_constraint, hirshfeld_constraint
89
90CONTAINS
91
92! **************************************************************************************************
93!> \brief Driver routine for calculating a Becke constraint
94!> \param qs_env the qs_env where to build the constraint
95!> \param calc_pot if the potential needs to be recalculated or just integrated
96!> \param calculate_forces logical if potential has to be calculated or only_energy
97!> \par History
98!> Created 01.2007 [fschiff]
99!> Extended functionality 12/15-12/16 [Nico Holmberg]
100! **************************************************************************************************
101 SUBROUTINE becke_constraint(qs_env, calc_pot, calculate_forces)
102 TYPE(qs_environment_type), POINTER :: qs_env
103 LOGICAL :: calc_pot, calculate_forces
104
105 CHARACTER(len=*), PARAMETER :: routinen = 'becke_constraint'
106
107 INTEGER :: handle
108 TYPE(cdft_control_type), POINTER :: cdft_control
109 TYPE(dft_control_type), POINTER :: dft_control
110
111 CALL timeset(routinen, handle)
112 CALL get_qs_env(qs_env, dft_control=dft_control)
113 cdft_control => dft_control%qs_control%cdft_control
114 IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_becke_constraint) THEN
115 IF (calc_pot) THEN
116 ! Initialize the Becke constraint environment
117 CALL becke_constraint_init(qs_env)
118 ! Calculate the Becke weight function and possibly the gradients
119 CALL becke_constraint_low(qs_env)
120 END IF
121 ! Integrate the smooth density contribution. GAPW one-center terms are added
122 ! while the Kohn-Sham matrix is assembled.
123 CALL cdft_constraint_integrate(qs_env)
124 IF (calculate_forces) CALL cdft_constraint_force(qs_env)
125 END IF
126 CALL timestop(handle)
127
128 END SUBROUTINE becke_constraint
129
130! **************************************************************************************************
131!> \brief Low level routine to build a Becke weight function and its gradients
132!> \param qs_env the qs_env where to build the constraint
133!> \param just_gradients optional logical which determines if only the gradients should be calculated
134!> \par History
135!> Created 03.2017 [Nico Holmberg]
136! **************************************************************************************************
137 SUBROUTINE becke_constraint_low(qs_env, just_gradients)
138 TYPE(qs_environment_type), POINTER :: qs_env
139 LOGICAL, OPTIONAL :: just_gradients
140
141 CHARACTER(len=*), PARAMETER :: routinen = 'becke_constraint_low'
142 REAL(kind=dp), PARAMETER :: eps_sum_cell_f_all = 1.0e-06_dp
143
144 INTEGER :: handle, i, iatom, igroup, ind(3), ip, j, &
145 jatom, jp, k, natom, np(3), nskipped
146 INTEGER, ALLOCATABLE, DIMENSION(:) :: catom
147 INTEGER, DIMENSION(2, 3) :: bo, bo_conf
148 LOGICAL :: in_memory, my_just_gradients
149 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint, skip_me
150 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: atom_in_group
151 REAL(kind=dp) :: dist1, dist2, dmyexp, dvol, eps_cavity, &
152 my1, my1_homo, myexp, sum_cell_f_all, &
153 th, tmp_const
154 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cell_functions, ds_dr_i, ds_dr_j, &
155 sum_cell_f_group
156 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_sum_pm_dr, dp_i_dri
157 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dp_i_drj
158 REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dmy_dr_i, dmy_dr_j, &
159 dr, dr1_r2, dr_i_dr, dr_ij_dr, &
160 dr_j_dr, grid_p, r, r1, shift
161 REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
162 TYPE(becke_constraint_type), POINTER :: becke_control
163 TYPE(cdft_control_type), POINTER :: cdft_control
164 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
165 TYPE(cell_type), POINTER :: cell
166 TYPE(dft_control_type), POINTER :: dft_control
167 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
168 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge
169
170 NULLIFY (cutoffs, cell, dft_control, particle_set, group, charge, cdft_control)
171 CALL timeset(routinen, handle)
172 ! Get simulation environment
173 CALL get_qs_env(qs_env, &
174 cell=cell, &
175 particle_set=particle_set, &
176 natom=natom, &
177 dft_control=dft_control)
178 cdft_control => dft_control%qs_control%cdft_control
179 becke_control => cdft_control%becke_control
180 group => cdft_control%group
181 cutoffs => becke_control%cutoffs
182 IF (cdft_control%atomic_charges) THEN
183 charge => cdft_control%charge
184 END IF
185 in_memory = .false.
186 IF (cdft_control%save_pot) THEN
187 in_memory = becke_control%in_memory
188 END IF
189 eps_cavity = becke_control%eps_cavity
190 ! Decide if only gradients need to be calculated
191 my_just_gradients = .false.
192 IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
193 IF (my_just_gradients) THEN
194 in_memory = .true.
195 ! Pairwise distances need to be recalculated
196 IF (becke_control%vector_buffer%store_vectors) THEN
197 ALLOCATE (becke_control%vector_buffer%distances(natom))
198 ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom))
199 IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom))
200 ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom))
201 END IF
202 ALLOCATE (becke_control%vector_buffer%R12(natom, natom))
203 DO i = 1, 3
204 cell_v(i) = cell%hmat(i, i)
205 END DO
206 DO iatom = 1, natom - 1
207 DO jatom = iatom + 1, natom
208 r = particle_set(iatom)%r
209 r1 = particle_set(jatom)%r
210 DO i = 1, 3
211 r(i) = modulo(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
212 r1(i) = modulo(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
213 END DO
214 dist_vec = (r - r1) - anint((r - r1)/cell_v)*cell_v
215 IF (becke_control%vector_buffer%store_vectors) THEN
216 becke_control%vector_buffer%position_vecs(:, iatom) = r(:)
217 IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:)
218 IF (in_memory) THEN
219 becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
220 becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
221 END IF
222 END IF
223 becke_control%vector_buffer%R12(iatom, jatom) = norm2(dist_vec)
224 becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom)
225 END DO
226 END DO
227 END IF
228 ALLOCATE (catom(cdft_control%natoms))
229 IF (cdft_control%save_pot .OR. &
230 becke_control%cavity_confine .OR. &
231 becke_control%should_skip) THEN
232 ALLOCATE (is_constraint(natom))
233 is_constraint = .false.
234 END IF
235 ! This boolean is needed to prevent calculation of atom pairs ji when the pair ij has
236 ! already been calculated (data for pair ji is set using symmetry)
237 ! With gradient precomputation, symmetry exploited for both weight function and gradients
238 ALLOCATE (skip_me(natom))
239 DO i = 1, cdft_control%natoms
240 catom(i) = cdft_control%atoms(i)
241 ! Notice that here is_constraint=.TRUE. also for dummy atoms to properly compute their Becke charges
242 ! A subsequent check (atom_in_group) ensures that the gradients of these dummy atoms are correct
243 IF (cdft_control%save_pot .OR. &
244 becke_control%cavity_confine .OR. &
245 becke_control%should_skip) THEN
246 is_constraint(catom(i)) = .true.
247 END IF
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) == 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 = norm2(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 = norm2(dist_vec)
350 END IF
351 IF (dist1 <= cutoffs(iatom)) THEN
352 IF (in_memory) THEN
353 IF (dist1 <= th) dist1 = th
354 dr_i_dr(:) = dist_vec(:)/dist1
355 END IF
356 DO jatom = 1, natom
357 IF (jatom /= 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) == 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 = norm2(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 = norm2(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 <= 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) <= 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 <= 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) <= 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 /= 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) > 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) > eps_sum_cell_f_all) 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)) THEN
566 DEALLOCATE (is_constraint)
567 END IF
568 DEALLOCATE (catom)
569 DEALLOCATE (cell_functions)
570 DEALLOCATE (skip_me)
571 DEALLOCATE (sum_cell_f_group)
572 DEALLOCATE (becke_control%vector_buffer%R12)
573 IF (becke_control%vector_buffer%store_vectors) THEN
574 DEALLOCATE (becke_control%vector_buffer%distances)
575 DEALLOCATE (becke_control%vector_buffer%distance_vecs)
576 DEALLOCATE (becke_control%vector_buffer%position_vecs)
577 END IF
578 CALL timestop(handle)
579
580 END SUBROUTINE becke_constraint_low
581
582! **************************************************************************************************
583!> \brief Driver routine for calculating a Hirshfeld constraint
584!> \param qs_env ...
585!> \param calc_pot ...
586!> \param calculate_forces ...
587! **************************************************************************************************
588 SUBROUTINE hirshfeld_constraint(qs_env, calc_pot, calculate_forces)
589 TYPE(qs_environment_type), POINTER :: qs_env
590 LOGICAL :: calc_pot, calculate_forces
591
592 CHARACTER(len=*), PARAMETER :: routinen = 'hirshfeld_constraint'
593
594 INTEGER :: handle
595 TYPE(cdft_control_type), POINTER :: cdft_control
596 TYPE(dft_control_type), POINTER :: dft_control
597
598 CALL timeset(routinen, handle)
599 CALL get_qs_env(qs_env, dft_control=dft_control)
600 cdft_control => dft_control%qs_control%cdft_control
601 IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_hirshfeld_constraint) THEN
602 IF (calc_pot) THEN
603 ! Initialize the Hirshfeld constraint environment
604 CALL hirshfeld_constraint_init(qs_env)
605 ! Calculate the Hirshfeld weight function and possibly the gradients
606 CALL hirshfeld_constraint_low(qs_env)
607 END IF
608 ! Integrate the smooth density contribution. GAPW one-center terms are added
609 ! while the Kohn-Sham matrix is assembled.
610 CALL cdft_constraint_integrate(qs_env)
611 IF (calculate_forces) CALL cdft_constraint_force(qs_env)
612 END IF
613 CALL timestop(handle)
614
615 END SUBROUTINE hirshfeld_constraint
616
617! **************************************************************************************************
618!> \brief Calculates Hirshfeld constraints
619!> \param qs_env ...
620!> \param just_gradients ...
621! **************************************************************************************************
622 SUBROUTINE hirshfeld_constraint_low(qs_env, just_gradients)
623 TYPE(qs_environment_type), POINTER :: qs_env
624 LOGICAL, OPTIONAL :: just_gradients
625
626 CHARACTER(len=*), PARAMETER :: routinen = 'hirshfeld_constraint_low'
627
628 INTEGER :: atom_a, atoms_memory, atoms_memory_num, handle, i, iatom, iex, igroup, ikind, &
629 ithread, j, k, natom, npme, nthread, num_atoms, num_species, numexp, subpatch_pattern
630 INTEGER, ALLOCATABLE, DIMENSION(:) :: num_species_small
631 INTEGER, DIMENSION(2, 3) :: bo
632 INTEGER, DIMENSION(3) :: lb_pw, lb_rs, ub_pw, ub_rs
633 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
634 LOGICAL :: my_just_gradients
635 LOGICAL, ALLOCATABLE, DIMENSION(:) :: compute_charge, is_constraint
636 REAL(kind=dp) :: alpha, coef, eps_rho_rspace, exp_eval, &
637 prefactor, radius
638 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients
639 REAL(kind=dp), DIMENSION(3) :: dr_rs, r2, r_pbc, ra
640 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
641 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
642 TYPE(cdft_control_type), POINTER :: cdft_control
643 TYPE(cell_type), POINTER :: cell
644 TYPE(dft_control_type), POINTER :: dft_control
645 TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control
646 TYPE(hirshfeld_type), POINTER :: hirshfeld_env
647 TYPE(mp_para_env_type), POINTER :: para_env
648 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
649 TYPE(pw_env_type), POINTER :: pw_env
650 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
651 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: pw_single_dr
652 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
653 TYPE(qs_rho_type), POINTER :: rho
654 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
655 TYPE(realspace_grid_type) :: rs_rho_all, rs_rho_constr
656 TYPE(realspace_grid_type), ALLOCATABLE, &
657 DIMENSION(:) :: rs_single, rs_single_charge, rs_single_dr
658
659 NULLIFY (atom_list, atomic_kind_set, dft_control, &
660 hirshfeld_env, particle_set, pw_env, auxbas_pw_pool, para_env, &
661 auxbas_rs_desc, cdft_control, pab, &
662 hirshfeld_control, cell, rho_r, rho)
663
664 CALL timeset(routinen, handle)
665 CALL get_qs_env(qs_env, &
666 atomic_kind_set=atomic_kind_set, &
667 particle_set=particle_set, &
668 natom=natom, &
669 cell=cell, &
670 rho=rho, &
671 dft_control=dft_control, &
672 para_env=para_env, &
673 pw_env=pw_env)
674 CALL qs_rho_get(rho, rho_r=rho_r)
675
676 num_atoms = natom
677
678 cdft_control => dft_control%qs_control%cdft_control
679 hirshfeld_control => cdft_control%hirshfeld_control
680 hirshfeld_env => hirshfeld_control%hirshfeld_env
681
682 ! Check if only gradient should be calculated, if gradients should be precomputed
683 my_just_gradients = .false.
684 IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
685 IF (my_just_gradients) THEN
686 cdft_control%in_memory = .true.
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 lb_pw(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
698 ub_pw(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
699
700 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
701 auxbas_pw_pool=auxbas_pw_pool)
702 CALL rs_grid_create(rs_rho_all, auxbas_rs_desc)
703 CALL rs_grid_zero(rs_rho_all)
704
705 dr_rs(1) = rs_rho_all%desc%dh(1, 1)
706 dr_rs(2) = rs_rho_all%desc%dh(2, 2)
707 dr_rs(3) = rs_rho_all%desc%dh(3, 3)
708 lb_rs(1) = lbound(rs_rho_all%r(:, :, :), 1)
709 lb_rs(2) = lbound(rs_rho_all%r(:, :, :), 2)
710 lb_rs(3) = lbound(rs_rho_all%r(:, :, :), 3)
711 ub_rs(1) = ubound(rs_rho_all%r(:, :, :), 1)
712 ub_rs(2) = ubound(rs_rho_all%r(:, :, :), 2)
713 ub_rs(3) = ubound(rs_rho_all%r(:, :, :), 3)
714
715 ! For each CDFT group
716 DO igroup = 1, SIZE(cdft_control%group)
717
718 IF (igroup == 2 .AND. .NOT. cdft_control%in_memory) THEN
719 CALL rs_grid_zero(rs_rho_all)
720 END IF
721 bo = cdft_control%group(igroup)%weight%pw_grid%bounds_local
722
723 ! Coefficients
724 coefficients(:) = 0.0_dp
725 is_constraint = .false.
726 DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
727 coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
728 is_constraint(cdft_control%group(igroup)%atoms(i)) = .true.
729 END DO
730
731 ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group
732 CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc)
733 CALL rs_grid_zero(rs_rho_constr)
734
735 ! rs_single: Gaussian density over single atoms when required
736 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
737 ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic(cdft_control%natoms))
738 ALLOCATE (rs_single(cdft_control%natoms))
739 DO i = 1, cdft_control%natoms
740 CALL rs_grid_create(rs_single(i), auxbas_rs_desc)
741 CALL rs_grid_zero(rs_single(i))
742 END DO
743 END IF
744
745 ! Setup pw
746 CALL pw_zero(cdft_control%group(igroup)%weight)
747
748 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
749 CALL pw_set(cdft_control%group(igroup)%hw_rho_total_constraint, 1.0_dp)
750
751 IF (igroup == 1) THEN
752 CALL auxbas_pw_pool%create_pw(cdft_control%hw_rho_total)
753 CALL pw_set(cdft_control%hw_rho_total, 1.0_dp)
754
755 IF (hirshfeld_control%print_density) THEN
756 DO iatom = 1, cdft_control%natoms
757 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic(iatom))
758 CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic(iatom), 1.0_dp)
759 END DO
760 END IF
761 END IF
762
763 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
764 ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge(cdft_control%natoms))
765 ALLOCATE (rs_single_charge(cdft_control%natoms))
766 ALLOCATE (compute_charge(natom))
767 compute_charge = .false.
768
769 DO i = 1, cdft_control%natoms
770 CALL rs_grid_create(rs_single_charge(i), auxbas_rs_desc)
771 CALL rs_grid_zero(rs_single_charge(i))
772 compute_charge(cdft_control%atoms(i)) = .true.
773 END DO
774
775 DO iatom = 1, cdft_control%natoms
776 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom))
777 CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom), 1.0_dp)
778 END DO
779 END IF
780
781 ALLOCATE (pab(1, 1))
782 nthread = 1
783 ithread = 0
784
785 DO ikind = 1, SIZE(atomic_kind_set)
786 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
787 IF (numexp <= 0) cycle
788 CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
789 ALLOCATE (cores(num_species))
790
791 DO iex = 1, numexp
792 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
793 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
794 npme = 0
795 cores = 0
796 DO iatom = 1, num_species
797 atom_a = atom_list(iatom)
798 ra(:) = pbc(particle_set(atom_a)%r, cell)
799 IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
800 IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
801 npme = npme + 1
802 cores(npme) = iatom
803 END IF
804 ELSE
805 npme = npme + 1
806 cores(npme) = iatom
807 END IF
808 END DO
809 DO j = 1, npme
810 iatom = cores(j)
811 atom_a = atom_list(iatom)
812 pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
813 ra(:) = pbc(particle_set(atom_a)%r, cell)
814
815 IF (hirshfeld_control%use_atomic_cutoff) THEN
816 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
817 ra=ra, rb=ra, rp=ra, &
818 zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
819 pab=pab, o1=0, o2=0, & ! without map_consistent
820 prefactor=1.0_dp, cutoff=0.0_dp)
821 END IF
822
823 IF (igroup == 1) THEN
824 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
825 [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
826 rs_rho_all, radius=radius, &
827 ga_gb_function=grid_func_ab, use_subpatch=.true., &
828 subpatch_pattern=subpatch_pattern)
829 END IF
830
831 IF (is_constraint(atom_a)) 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], coefficients(atom_a), &
834 pab, 0, 0, rs_rho_constr, &
835 radius=radius, &
836 ga_gb_function=grid_func_ab, use_subpatch=.true., &
837 subpatch_pattern=subpatch_pattern)
838 END IF
839
840 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
841 IF (is_constraint(atom_a)) THEN
842 DO iatom = 1, cdft_control%natoms
843 IF (atom_a == cdft_control%atoms(iatom)) EXIT
844 END DO
845 cpassert(iatom <= cdft_control%natoms)
846 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
847 [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
848 rs_single(iatom), radius=radius, &
849 ga_gb_function=grid_func_ab, use_subpatch=.true., &
850 subpatch_pattern=subpatch_pattern)
851 END IF
852 END IF
853
854 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
855 IF (compute_charge(atom_a)) THEN
856 DO iatom = 1, cdft_control%natoms
857 IF (atom_a == cdft_control%atoms(iatom)) EXIT
858 END DO
859 cpassert(iatom <= cdft_control%natoms)
860 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
861 [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
862 rs_single_charge(iatom), radius=radius, &
863 ga_gb_function=grid_func_ab, use_subpatch=.true., &
864 subpatch_pattern=subpatch_pattern)
865 END IF
866 END IF
867
868 END DO
869 END DO
870 DEALLOCATE (cores)
871 END DO
872 DEALLOCATE (pab)
873
874 IF (igroup == 1) THEN
875 CALL transfer_rs2pw(rs_rho_all, cdft_control%hw_rho_total)
876 END IF
877
878 CALL transfer_rs2pw(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint)
879 CALL rs_grid_release(rs_rho_constr)
880
881 ! Calculate weight function
882 CALL hfun_scale(cdft_control%group(igroup)%weight%array, &
883 cdft_control%group(igroup)%hw_rho_total_constraint%array, &
884 cdft_control%hw_rho_total%array, divide=.true., &
885 small=hirshfeld_control%eps_cutoff)
886
887 ! Calculate charges
888 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
889 DO i = 1, cdft_control%natoms
890 CALL transfer_rs2pw(rs_single_charge(i), cdft_control%group(igroup)%hw_rho_atomic_charge(i))
891 CALL hfun_scale(cdft_control%charge(i)%array, &
892 cdft_control%group(igroup)%hw_rho_atomic_charge(i)%array, &
893 cdft_control%hw_rho_total%array, divide=.true., &
894 small=hirshfeld_control%eps_cutoff)
895 END DO
896 END IF
897
898 ! Print atomic densities if requested
899 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
900 DO i = 1, cdft_control%natoms
901 CALL transfer_rs2pw(rs_single(i), cdft_control%group(igroup)%hw_rho_atomic(i))
902 END DO
903 CALL cdft_print_hirshfeld_density(qs_env)
904 END IF
905
906 END DO
907
908 DO igroup = 1, SIZE(cdft_control%group)
909
910 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
911
912 IF (.NOT. cdft_control%in_memory .AND. igroup == 1) THEN
913 CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
914 END IF
915
916 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
917 DO i = 1, cdft_control%natoms
918 CALL rs_grid_release(rs_single(i))
919 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic(i))
920 END DO
921 DEALLOCATE (rs_single)
922 DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic)
923 END IF
924
925 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
926 DO i = 1, cdft_control%natoms
927 CALL rs_grid_release(rs_single_charge(i))
928 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(i))
929 END DO
930 DEALLOCATE (rs_single_charge)
931 DEALLOCATE (compute_charge)
932 DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge)
933 END IF
934
935 END DO
936
937 IF (cdft_control%in_memory) THEN
938 DO igroup = 1, SIZE(cdft_control%group)
939 ALLOCATE (cdft_control%group(igroup)%gradients_x(1*natom, lb_pw(1):ub_pw(1), &
940 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
941 cdft_control%group(igroup)%gradients_x(:, :, :, :) = 0.0_dp
942 END DO
943 END IF
944
945 IF (cdft_control%in_memory) THEN
946 DO igroup = 1, SIZE(cdft_control%group)
947
948 ALLOCATE (pab(1, 1))
949 nthread = 1
950 ithread = 0
951 atoms_memory = hirshfeld_control%atoms_memory
952
953 DO ikind = 1, SIZE(atomic_kind_set)
954 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
955 IF (numexp <= 0) cycle
956 CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
957
958 ALLOCATE (pw_single_dr(num_species))
959 ALLOCATE (rs_single_dr(num_species))
960
961 DO i = 1, num_species
962 CALL auxbas_pw_pool%create_pw(pw_single_dr(i))
963 CALL pw_zero(pw_single_dr(i))
964 END DO
965
966 atoms_memory_num = SIZE([(j, j=1, num_species, atoms_memory)])
967
968 ! Can't store all pw grids, therefore split into groups of size atom_memory
969 ! Ideally this code should be re-written to be more memory efficient
970 IF (num_species > atoms_memory) THEN
971 ALLOCATE (num_species_small(atoms_memory_num + 1))
972 num_species_small(1:atoms_memory_num) = [(j, j=1, num_species, atoms_memory)]
973 num_species_small(atoms_memory_num + 1) = num_species
974 ELSE
975 ALLOCATE (num_species_small(2))
976 num_species_small(:) = [1, num_species]
977 END IF
978
979 DO k = 1, SIZE(num_species_small) - 1
980 IF (num_species > atoms_memory) THEN
981 ALLOCATE (cores(num_species_small(k + 1) - (num_species_small(k) - 1)))
982 ELSE
983 ALLOCATE (cores(num_species))
984 END IF
985
986 DO i = num_species_small(k), num_species_small(k + 1)
987 CALL rs_grid_create(rs_single_dr(i), auxbas_rs_desc)
988 CALL rs_grid_zero(rs_single_dr(i))
989 END DO
990 DO iex = 1, numexp
991
992 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
993 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
994 prefactor = 2.0_dp*alpha
995 npme = 0
996 cores = 0
997
998 DO iatom = 1, SIZE(cores)
999 atom_a = atom_list(iatom + (num_species_small(k) - 1))
1000 ra(:) = pbc(particle_set(atom_a)%r, cell)
1001
1002 IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
1003 IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
1004 npme = npme + 1
1005 cores(npme) = iatom
1006 END IF
1007 ELSE
1008 npme = npme + 1
1009 cores(npme) = iatom
1010 END IF
1011 END DO
1012 DO j = 1, npme
1013 iatom = cores(j)
1014 atom_a = atom_list(iatom + (num_species_small(k) - 1))
1015 pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
1016 ra(:) = pbc(particle_set(atom_a)%r, cell)
1017 subpatch_pattern = 0
1018
1019 ! Calculate cutoff
1020 IF (hirshfeld_control%use_atomic_cutoff) THEN
1021 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
1022 ra=ra, rb=ra, rp=ra, &
1023 zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
1024 pab=pab, o1=0, o2=0, & ! without map_consistent
1025 prefactor=1.0_dp, cutoff=0.0_dp)
1026 END IF
1027
1028 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1029 [0.0_dp, 0.0_dp, 0.0_dp], prefactor, &
1030 pab, 0, 0, rs_single_dr(iatom + (num_species_small(k) - 1)), &
1031 radius=radius, &
1032 ga_gb_function=grid_func_ab, use_subpatch=.true., &
1033 subpatch_pattern=subpatch_pattern)
1034
1035 END DO
1036 END DO
1037
1038 DO iatom = num_species_small(k), num_species_small(k + 1)
1039 CALL transfer_rs2pw(rs_single_dr(iatom), pw_single_dr(iatom))
1040 CALL rs_grid_release(rs_single_dr(iatom))
1041 END DO
1042
1043 DEALLOCATE (cores)
1044 END DO
1045
1046 DO iatom = 1, num_species
1047 atom_a = atom_list(iatom)
1048 cdft_control%group(igroup)%gradients_x(atom_a, :, :, :) = pw_single_dr(iatom)%array(:, :, :)
1049 CALL auxbas_pw_pool%give_back_pw(pw_single_dr(iatom))
1050 END DO
1051
1052 DEALLOCATE (rs_single_dr)
1053 DEALLOCATE (num_species_small)
1054 DEALLOCATE (pw_single_dr)
1055 END DO
1056 DEALLOCATE (pab)
1057 END DO
1058 END IF
1059
1060 IF (cdft_control%in_memory) THEN
1061 DO igroup = 1, SIZE(cdft_control%group)
1062 ALLOCATE (cdft_control%group(igroup)%gradients_y(1*num_atoms, lb_pw(1):ub_pw(1), &
1063 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1064 ALLOCATE (cdft_control%group(igroup)%gradients_z(1*num_atoms, lb_pw(1):ub_pw(1), &
1065 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1066 cdft_control%group(igroup)%gradients_y(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1067 cdft_control%group(igroup)%gradients_z(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1068 END DO
1069 END IF
1070
1071 ! Calculate gradient if requested
1072 IF (cdft_control%in_memory) THEN
1073
1074 DO igroup = 1, SIZE(cdft_control%group)
1075
1076 ! Coefficients
1077 coefficients(:) = 0.0_dp
1078 is_constraint = .false.
1079 DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
1080 coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
1081 is_constraint(cdft_control%group(igroup)%atoms(i)) = .true.
1082 END DO
1083
1084 DO k = lb_pw(3), ub_pw(3)
1085 DO j = lb_pw(2), ub_pw(2)
1086 DO i = lb_pw(1), ub_pw(1)
1087 DO iatom = 1, natom
1088
1089 ra(:) = particle_set(iatom)%r
1090
1091 IF (cdft_control%hw_rho_total%array(i, j, k) > hirshfeld_control%eps_cutoff) THEN
1092
1093 exp_eval = (coefficients(iatom) - &
1094 cdft_control%group(igroup)%weight%array(i, j, k))/ &
1095 cdft_control%hw_rho_total%array(i, j, k)
1096
1097 r2 = real(i - rho_r(1)%pw_grid%bounds(1, 1), dp)* &
1098 rho_r(1)%pw_grid%dh(:, 1) + &
1099 REAL(j - rho_r(1)%pw_grid%bounds(1, 2), dp)* &
1100 rho_r(1)%pw_grid%dh(:, 2) + &
1101 REAL(k - rho_r(1)%pw_grid%bounds(1, 3), dp)* &
1102 rho_r(1)%pw_grid%dh(:, 3)
1103 r_pbc = pbc(ra, r2, cell)
1104
1105 ! Store gradient d/dR_x w, including term: (r_x - R_x)
1106 cdft_control%group(igroup)%gradients_x(iatom, i, j, k) = &
1107 cdft_control%group(igroup)%gradients_x(iatom, i, j, k)* &
1108 r_pbc(1)*exp_eval
1109
1110 ! Store gradient d/dR_y w, including term: (r_y - R_y)
1111 cdft_control%group(igroup)%gradients_y(iatom, i, j, k) = &
1112 cdft_control%group(igroup)%gradients_y(iatom, i, j, k)* &
1113 r_pbc(2)*exp_eval
1114
1115 ! Store gradient d/dR_z w, including term:(r_z - R_z)
1116 cdft_control%group(igroup)%gradients_z(iatom, i, j, k) = &
1117 cdft_control%group(igroup)%gradients_z(iatom, i, j, k)* &
1118 r_pbc(3)*exp_eval
1119
1120 END IF
1121 END DO
1122 END DO
1123 END DO
1124 END DO
1125 END DO
1126 CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
1127 END IF
1128
1129 CALL rs_grid_release(rs_rho_all)
1130
1131 IF (ALLOCATED(coefficients)) DEALLOCATE (coefficients)
1132 IF (ALLOCATED(is_constraint)) DEALLOCATE (is_constraint)
1133
1134 CALL timestop(handle)
1135
1136 END SUBROUTINE hirshfeld_constraint_low
1137
1138! **************************************************************************************************
1139!> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT
1140!> weight function and the realspace electron density
1141!> \param qs_env ...
1142! **************************************************************************************************
1143 SUBROUTINE cdft_constraint_integrate(qs_env)
1144 TYPE(qs_environment_type), POINTER :: qs_env
1145
1146 CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_integrate'
1147
1148 INTEGER :: handle, i, iatom, igroup, ivar, iw, nvar
1149 LOGICAL :: is_becke
1150 REAL(kind=dp) :: dvol, eps_cavity, sign
1151 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: de, strength, target_val
1152 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge
1153 TYPE(becke_constraint_type), POINTER :: becke_control
1154 TYPE(cdft_control_type), POINTER :: cdft_control
1155 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1156 TYPE(cp_logger_type), POINTER :: logger
1157 TYPE(dft_control_type), POINTER :: dft_control
1158 TYPE(mp_para_env_type), POINTER :: para_env
1159 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge, rho_r
1160 TYPE(qs_energy_type), POINTER :: energy
1161 TYPE(qs_rho_type), POINTER :: rho
1162 TYPE(section_vals_type), POINTER :: cdft_constraint_section
1163
1164 NULLIFY (para_env, dft_control, rho_r, energy, rho, &
1165 logger, cdft_constraint_section, group, charge)
1166 CALL timeset(routinen, handle)
1167 logger => cp_get_default_logger()
1168 CALL get_qs_env(qs_env, &
1169 rho=rho, &
1170 dft_control=dft_control, &
1171 para_env=para_env)
1172 CALL qs_rho_get(rho, rho_r=rho_r)
1173 cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1174 iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
1175 cdft_control => dft_control%qs_control%cdft_control
1176 is_becke = (cdft_control%type == outer_scf_becke_constraint)
1177 becke_control => cdft_control%becke_control
1178 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) THEN
1179 cpabort("Becke control has not been allocated.")
1180 END IF
1181 group => cdft_control%group
1182 ! Initialize
1183 nvar = SIZE(cdft_control%target)
1184 ALLOCATE (strength(nvar))
1185 ALLOCATE (target_val(nvar))
1186 ALLOCATE (de(nvar))
1187 strength(:) = cdft_control%strength(:)
1188 sign = 1.0_dp
1189 de = 0.0_dp
1190 dvol = group(1)%weight%pw_grid%dvol
1191 IF (cdft_control%atomic_charges) THEN
1192 charge => cdft_control%charge
1193 ALLOCATE (electronic_charge(cdft_control%natoms, dft_control%nspins))
1194 electronic_charge = 0.0_dp
1195 END IF
1196 ! Calculate value of constraint i.e. int ( rho(r) w(r) dr)
1197 DO i = 1, dft_control%nspins
1198 DO igroup = 1, SIZE(group)
1199 SELECT CASE (group(igroup)%constraint_type)
1200 CASE (cdft_charge_constraint)
1201 sign = 1.0_dp
1202 CASE (cdft_magnetization_constraint)
1203 IF (i == 1) THEN
1204 sign = 1.0_dp
1205 ELSE
1206 sign = -1.0_dp
1207 END IF
1208 CASE (cdft_alpha_constraint)
1209 sign = 1.0_dp
1210 IF (i == 2) cycle
1211 CASE (cdft_beta_constraint)
1212 sign = 1.0_dp
1213 IF (i == 1) cycle
1214 CASE DEFAULT
1215 cpabort("Unknown constraint type.")
1216 END SELECT
1217 IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1218 ! With external control, we can use cavity_mat as a mask to kahan sum
1219 eps_cavity = becke_control%eps_cavity
1220 IF (igroup /= 1) THEN
1221 CALL cp_abort(__location__, &
1222 "Multiple constraints not yet supported by parallel mixed calculations.")
1223 END IF
1224 de(igroup) = de(igroup) + sign*accurate_dot_product(group(igroup)%weight%array, rho_r(i)%array, &
1225 becke_control%cavity_mat, eps_cavity)*dvol
1226 ELSE
1227 de(igroup) = de(igroup) + sign*pw_integral_ab(group(igroup)%weight, rho_r(i), local_only=.true.)
1228 END IF
1229 END DO
1230 IF (cdft_control%atomic_charges) THEN
1231 DO iatom = 1, cdft_control%natoms
1232 electronic_charge(iatom, i) = pw_integral_ab(charge(iatom), rho_r(i), local_only=.true.)
1233 END DO
1234 END IF
1235 END DO
1236 CALL get_qs_env(qs_env, energy=energy)
1237 CALL para_env%sum(de)
1238 IF (cdft_control%atomic_charges) THEN
1239 CALL para_env%sum(electronic_charge)
1240 END IF
1241 ! Use fragment densities as reference value (= Becke deformation density)
1242 IF (cdft_control%fragment_density .AND. .NOT. cdft_control%fragments_integrated) THEN
1243 CALL prepare_fragment_constraint(qs_env)
1244 END IF
1245 ! Update constraint value and energy
1246 ! Fragment integration above can replace the input target.
1247 target_val(:) = cdft_control%target(:)
1248 cdft_control%value(:) = de(:)
1249 energy%cdft = 0.0_dp
1250 DO ivar = 1, nvar
1251 energy%cdft = energy%cdft + (de(ivar) - target_val(ivar))*strength(ivar)
1252 END DO
1253 ! Print constraint info and atomic CDFT charges
1254 IF (.NOT. dft_control%qs_control%gapw) THEN
1255 CALL cdft_constraint_print(qs_env, electronic_charge)
1256 END IF
1257 ! Deallocate tmp storage
1258 DEALLOCATE (de, strength, target_val)
1259 IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge)
1260 CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO")
1261 CALL timestop(handle)
1262
1263 END SUBROUTINE cdft_constraint_integrate
1264
1265! **************************************************************************************************
1266!> \brief Calculates atomic forces due to a CDFT constraint (Becke or Hirshfeld)
1267!> \param qs_env ...
1268! **************************************************************************************************
1269 SUBROUTINE cdft_constraint_force(qs_env)
1270 TYPE(qs_environment_type), POINTER :: qs_env
1271
1272 CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_force'
1273
1274 INTEGER :: handle, i, iatom, igroup, ikind, ispin, &
1275 j, k, natom, nvar
1276 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
1277 INTEGER, DIMENSION(2, 3) :: bo
1278 INTEGER, DIMENSION(3) :: lb, ub
1279 REAL(kind=dp) :: dvol, eps_cavity, sign
1280 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength
1281 REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
1282 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1283 TYPE(becke_constraint_type), POINTER :: becke_control
1284 TYPE(cdft_control_type), POINTER :: cdft_control
1285 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1286 TYPE(cell_type), POINTER :: cell
1287 TYPE(dft_control_type), POINTER :: dft_control
1288 TYPE(mp_para_env_type), POINTER :: para_env
1289 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1290 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1291 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1292 TYPE(qs_rho_type), POINTER :: rho
1293
1294 CALL timeset(routinen, handle)
1295 NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, &
1296 rho, rho_r, force, cutoffs, becke_control, group)
1297
1298 CALL get_qs_env(qs_env, &
1299 atomic_kind_set=atomic_kind_set, &
1300 natom=natom, &
1301 particle_set=particle_set, &
1302 cell=cell, &
1303 rho=rho, &
1304 force=force, &
1305 dft_control=dft_control, &
1306 para_env=para_env)
1307 CALL qs_rho_get(rho, rho_r=rho_r)
1308
1309 cdft_control => dft_control%qs_control%cdft_control
1310 becke_control => cdft_control%becke_control
1311 group => cdft_control%group
1312 nvar = SIZE(cdft_control%target)
1313 ALLOCATE (strength(nvar))
1314 strength(:) = cdft_control%strength(:)
1315 cutoffs => cdft_control%becke_control%cutoffs
1316 eps_cavity = cdft_control%becke_control%eps_cavity
1317
1318 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1319 atom_of_kind=atom_of_kind, &
1320 kind_of=kind_of)
1321 DO igroup = 1, SIZE(cdft_control%group)
1322 ALLOCATE (cdft_control%group(igroup)%integrated(3, natom))
1323 cdft_control%group(igroup)%integrated = 0.0_dp
1324 END DO
1325
1326 lb(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
1327 ub(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
1328 bo = cdft_control%group(1)%weight%pw_grid%bounds_local
1329 dvol = cdft_control%group(1)%weight%pw_grid%dvol
1330 sign = 1.0_dp
1331
1332 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1333 IF (.NOT. cdft_control%becke_control%in_memory) THEN
1334 CALL becke_constraint_low(qs_env, just_gradients=.true.)
1335 END IF
1336
1337 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1338 IF (.NOT. cdft_control%in_memory) THEN
1339 CALL hirshfeld_constraint_low(qs_env, just_gradients=.true.)
1340 END IF
1341 END IF
1342
1343 ! If no Becke Gaussian confinement
1344 IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN
1345 ! No external control
1346 DO k = bo(1, 1), bo(2, 1)
1347 DO j = bo(1, 2), bo(2, 2)
1348 DO i = bo(1, 3), bo(2, 3)
1349 ! First check if this grid point should be skipped
1350 IF (cdft_control%becke_control%cavity_confine) THEN
1351 IF (cdft_control%becke_control%cavity%array(k, j, i) < eps_cavity) cycle
1352 END IF
1353
1354 DO igroup = 1, SIZE(cdft_control%group)
1355 DO iatom = 1, natom
1356 DO ispin = 1, dft_control%nspins
1357
1358 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1359 CASE (cdft_charge_constraint)
1360 sign = 1.0_dp
1361 CASE (cdft_magnetization_constraint)
1362 IF (ispin == 1) THEN
1363 sign = 1.0_dp
1364 ELSE
1365 sign = -1.0_dp
1366 END IF
1367 CASE (cdft_alpha_constraint)
1368 sign = 1.0_dp
1369 IF (ispin == 2) cycle
1370 CASE (cdft_beta_constraint)
1371 sign = 1.0_dp
1372 IF (ispin == 1) cycle
1373 CASE DEFAULT
1374 cpabort("Unknown constraint type.")
1375 END SELECT
1376
1377 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1378
1379 cdft_control%group(igroup)%integrated(:, iatom) = &
1380 cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1381 cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1382 *rho_r(ispin)%array(k, j, i) &
1383 *dvol
1384
1385 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1386
1387 cdft_control%group(igroup)%integrated(1, iatom) = &
1388 cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1389 cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1390 *rho_r(ispin)%array(k, j, i) &
1391 *dvol
1392
1393 cdft_control%group(igroup)%integrated(2, iatom) = &
1394 cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1395 cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1396 *rho_r(ispin)%array(k, j, i) &
1397 *dvol
1398
1399 cdft_control%group(igroup)%integrated(3, iatom) = &
1400 cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1401 cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1402 *rho_r(ispin)%array(k, j, i) &
1403 *dvol
1404
1405 END IF
1406
1407 END DO
1408 END DO
1409 END DO
1410 END DO
1411 END DO
1412 END DO
1413
1414 ! If Becke Gaussian confinement
1415 ELSE
1416 DO k = lbound(cdft_control%becke_control%cavity_mat, 1), ubound(cdft_control%becke_control%cavity_mat, 1)
1417 DO j = lbound(cdft_control%becke_control%cavity_mat, 2), ubound(cdft_control%becke_control%cavity_mat, 2)
1418 DO i = lbound(cdft_control%becke_control%cavity_mat, 3), ubound(cdft_control%becke_control%cavity_mat, 3)
1419
1420 ! First check if this grid point should be skipped
1421 IF (cdft_control%becke_control%cavity_mat(k, j, i) < eps_cavity) cycle
1422
1423 DO igroup = 1, SIZE(group)
1424 DO iatom = 1, natom
1425 DO ispin = 1, dft_control%nspins
1426 SELECT CASE (group(igroup)%constraint_type)
1427 CASE (cdft_charge_constraint)
1428 sign = 1.0_dp
1429 CASE (cdft_magnetization_constraint)
1430 IF (ispin == 1) THEN
1431 sign = 1.0_dp
1432 ELSE
1433 sign = -1.0_dp
1434 END IF
1435 CASE (cdft_alpha_constraint)
1436 sign = 1.0_dp
1437 IF (ispin == 2) cycle
1438 CASE (cdft_beta_constraint)
1439 sign = 1.0_dp
1440 IF (ispin == 1) cycle
1441 CASE DEFAULT
1442 cpabort("Unknown constraint type.")
1443 END SELECT
1444
1445 ! Integrate gradient of weight function
1446 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1447
1448 cdft_control%group(igroup)%integrated(:, iatom) = &
1449 cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1450 cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1451 *rho_r(ispin)%array(k, j, i) &
1452 *dvol
1453
1454 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1455
1456 cdft_control%group(igroup)%integrated(1, iatom) = &
1457 cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1458 cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1459 *rho_r(ispin)%array(k, j, i) &
1460 *dvol
1461
1462 cdft_control%group(igroup)%integrated(2, iatom) = &
1463 cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1464 cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1465 *rho_r(ispin)%array(k, j, i) &
1466 *dvol
1467
1468 cdft_control%group(igroup)%integrated(3, iatom) = &
1469 cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1470 cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1471 *rho_r(ispin)%array(k, j, i) &
1472 *dvol
1473
1474 END IF
1475
1476 END DO
1477 END DO
1478 END DO
1479 END DO
1480 END DO
1481 END DO
1482 END IF
1483
1484 IF (.NOT. cdft_control%transfer_pot) THEN
1485 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1486 DO igroup = 1, SIZE(group)
1487 DEALLOCATE (cdft_control%group(igroup)%gradients)
1488 END DO
1489 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1490 DO igroup = 1, SIZE(group)
1491 DEALLOCATE (cdft_control%group(igroup)%gradients_x)
1492 DEALLOCATE (cdft_control%group(igroup)%gradients_y)
1493 DEALLOCATE (cdft_control%group(igroup)%gradients_z)
1494 END DO
1495 END IF
1496 END IF
1497
1498 DO igroup = 1, SIZE(group)
1499 CALL para_env%sum(group(igroup)%integrated)
1500 END DO
1501
1502 ! Update force only on master process. Otherwise force due to constraint becomes multiplied
1503 ! by the number of processes when the final force%rho_elec is constructed in qs_force
1504 ! by mp_summing [the final integrated(:,:) is distributed on all processors]
1505 IF (para_env%is_source()) THEN
1506 DO igroup = 1, SIZE(group)
1507 DO iatom = 1, natom
1508 ikind = kind_of(iatom)
1509 i = atom_of_kind(iatom)
1510 force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup)
1511 END DO
1512 END DO
1513 END IF
1514
1515 DEALLOCATE (strength)
1516 DO igroup = 1, SIZE(group)
1517 DEALLOCATE (group(igroup)%integrated)
1518 END DO
1519 NULLIFY (group)
1520
1521 CALL timestop(handle)
1522
1523 END SUBROUTINE cdft_constraint_force
1524
1525! **************************************************************************************************
1526!> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied
1527!> by the CDFT weight functions and integrated over the realspace grid.
1528!> \param qs_env ...
1529! **************************************************************************************************
1530 SUBROUTINE prepare_fragment_constraint(qs_env)
1531 TYPE(qs_environment_type), POINTER :: qs_env
1532
1533 CHARACTER(len=*), PARAMETER :: routinen = 'prepare_fragment_constraint'
1534
1535 INTEGER :: handle, i, iatom, igroup, natom, &
1536 nelectron_total, nfrag_spins
1537 LOGICAL :: is_becke, needs_spin_density
1538 REAL(kind=dp) :: dvol, multiplier(2), nelectron_frag
1539 TYPE(becke_constraint_type), POINTER :: becke_control
1540 TYPE(cdft_control_type), POINTER :: cdft_control
1541 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1542 TYPE(cp_logger_type), POINTER :: logger
1543 TYPE(dft_control_type), POINTER :: dft_control
1544 TYPE(mp_para_env_type), POINTER :: para_env
1545 TYPE(pw_env_type), POINTER :: pw_env
1546 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1547 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: rho_frag
1548 TYPE(qs_subsys_type), POINTER :: subsys
1549
1550 NULLIFY (para_env, dft_control, logger, subsys, pw_env, auxbas_pw_pool, group)
1551 CALL timeset(routinen, handle)
1552 logger => cp_get_default_logger()
1553 CALL get_qs_env(qs_env, &
1554 natom=natom, &
1555 dft_control=dft_control, &
1556 para_env=para_env)
1557
1558 cdft_control => dft_control%qs_control%cdft_control
1559 is_becke = (cdft_control%type == outer_scf_becke_constraint)
1560 becke_control => cdft_control%becke_control
1561 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) THEN
1562 cpabort("Becke control has not been allocated.")
1563 END IF
1564 group => cdft_control%group
1565 dvol = group(1)%weight%pw_grid%dvol
1566 ! Fragment densities are meaningful only for some calculation types
1567 IF (.NOT. qs_env%single_point_run) THEN
1568 CALL cp_abort(__location__, &
1569 "CDFT fragment constraints are only compatible with single "// &
1570 "point calculations (run_type ENERGY or ENERGY_FORCE).")
1571 END IF
1572 needs_spin_density = .false.
1573 multiplier = 1.0_dp
1574 nfrag_spins = 1
1575 DO igroup = 1, SIZE(group)
1576 SELECT CASE (group(igroup)%constraint_type)
1577 CASE (cdft_charge_constraint)
1578 ! Do nothing
1579 CASE (cdft_magnetization_constraint)
1580 needs_spin_density = .true.
1581 CASE (cdft_alpha_constraint, cdft_beta_constraint)
1582 CALL cp_abort(__location__, &
1583 "CDFT fragment constraint not yet compatible with "// &
1584 "spin specific constraints.")
1585 CASE DEFAULT
1586 cpabort("Unknown constraint type.")
1587 END SELECT
1588 END DO
1589 IF (needs_spin_density) THEN
1590 nfrag_spins = 2
1591 DO i = 1, 2
1592 IF (cdft_control%flip_fragment(i)) multiplier(i) = -1.0_dp
1593 END DO
1594 END IF
1595 ! Read fragment reference densities
1596 ALLOCATE (cdft_control%fragments(nfrag_spins, 2))
1597 ALLOCATE (rho_frag(nfrag_spins))
1598 CALL get_qs_env(qs_env, pw_env=pw_env)
1599 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1600 IF (dft_control%qs_control%gapw) THEN
1601 CALL verify_gapw_fragment_cube(cdft_control%fragment_a_fname, spin_density=.false.)
1602 CALL verify_gapw_fragment_cube(cdft_control%fragment_b_fname, spin_density=.false.)
1603 IF (needs_spin_density) THEN
1604 CALL verify_gapw_fragment_cube(cdft_control%fragment_a_spin_fname, spin_density=.true.)
1605 CALL verify_gapw_fragment_cube(cdft_control%fragment_b_spin_fname, spin_density=.true.)
1606 END IF
1607 END IF
1608 ! Total density (rho_alpha + rho_beta)
1609 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 1))
1610 CALL cp_cube_to_pw(cdft_control%fragments(1, 1), &
1611 cdft_control%fragment_a_fname, 1.0_dp)
1612 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 2))
1613 CALL cp_cube_to_pw(cdft_control%fragments(1, 2), &
1614 cdft_control%fragment_b_fname, 1.0_dp)
1615 ! Spin difference density (rho_alpha - rho_beta) if needed
1616 IF (needs_spin_density) THEN
1617 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 1))
1618 CALL cp_cube_to_pw(cdft_control%fragments(2, 1), &
1619 cdft_control%fragment_a_spin_fname, multiplier(1))
1620 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 2))
1621 CALL cp_cube_to_pw(cdft_control%fragments(2, 2), &
1622 cdft_control%fragment_b_spin_fname, multiplier(2))
1623 END IF
1624 ! Sum up fragments
1625 DO i = 1, nfrag_spins
1626 CALL auxbas_pw_pool%create_pw(rho_frag(i))
1627 CALL pw_copy(cdft_control%fragments(i, 1), rho_frag(i))
1628 CALL pw_axpy(cdft_control%fragments(i, 2), rho_frag(i), 1.0_dp)
1629 CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 1))
1630 CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 2))
1631 END DO
1632 DEALLOCATE (cdft_control%fragments)
1633 ! Check that the number of electrons is consistent
1634 CALL get_qs_env(qs_env, subsys=subsys)
1635 CALL qs_subsys_get(subsys, nelectron_total=nelectron_total)
1636 nelectron_frag = pw_integrate_function(rho_frag(1))
1637 IF (nint(nelectron_frag) /= nelectron_total) THEN
1638 CALL cp_abort(__location__, &
1639 "The number of electrons in the reference and interacting "// &
1640 "configurations does not match. Check your fragment cube files.")
1641 END IF
1642 ! Update constraint target value i.e. perform integration w_i*rho_frag_{tot/spin}*dr
1643 cdft_control%target = 0.0_dp
1644 DO igroup = 1, SIZE(group)
1645 IF (group(igroup)%constraint_type == cdft_charge_constraint) THEN
1646 i = 1
1647 ELSE
1648 i = 2
1649 END IF
1650 IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1651 cdft_control%target(igroup) = cdft_control%target(igroup) + &
1652 accurate_dot_product(group(igroup)%weight%array, rho_frag(i)%array, &
1653 becke_control%cavity_mat, becke_control%eps_cavity)*dvol
1654 ELSE
1655 cdft_control%target(igroup) = cdft_control%target(igroup) + &
1656 pw_integral_ab(group(igroup)%weight, rho_frag(i), local_only=.true.)
1657 END IF
1658 END DO
1659 CALL para_env%sum(cdft_control%target)
1660 ! Calculate reference atomic charges int( w_i * rho_frag * dr )
1661 IF (cdft_control%atomic_charges) THEN
1662 ALLOCATE (cdft_control%charges_fragment(cdft_control%natoms, nfrag_spins))
1663 DO i = 1, nfrag_spins
1664 DO iatom = 1, cdft_control%natoms
1665 cdft_control%charges_fragment(iatom, i) = &
1666 pw_integral_ab(cdft_control%charge(iatom), rho_frag(i), local_only=.true.)
1667 END DO
1668 END DO
1669 CALL para_env%sum(cdft_control%charges_fragment)
1670 END IF
1671 DO i = 1, nfrag_spins
1672 CALL auxbas_pw_pool%give_back_pw(rho_frag(i))
1673 END DO
1674 DEALLOCATE (rho_frag)
1675 cdft_control%fragments_integrated = .true.
1676
1677 CALL timestop(handle)
1678
1679 CONTAINS
1680
1681! **************************************************************************************************
1682!> \brief Require full-density cubes for GAPW fragment references.
1683!> \param filename cube file to inspect
1684!> \param spin_density whether a total spin-density title is expected
1685! **************************************************************************************************
1686 SUBROUTINE verify_gapw_fragment_cube(filename, spin_density)
1687 CHARACTER(LEN=*), INTENT(IN) :: filename
1688 LOGICAL, INTENT(IN) :: spin_density
1689
1690 CHARACTER(LEN=256) :: title
1691 INTEGER :: input_unit, io_status
1692 LOGICAL :: usable_cube
1693
1694 usable_cube = .false.
1695 IF (para_env%is_source()) THEN
1696 CALL open_file(file_name=trim(filename), file_status="OLD", &
1697 file_action="READ", unit_number=input_unit)
1698 READ (input_unit, '(A)', iostat=io_status) title
1699 IF (io_status == 0) THEN
1700 READ (input_unit, '(A)', iostat=io_status) title
1701 usable_cube = io_status == 0
1702 END IF
1703 CALL close_file(input_unit)
1704 IF (usable_cube) THEN
1705 IF (spin_density) THEN
1706 IF (index(title, "SPIN DENSITY") > 0 .AND. &
1707 index(title, "TOTAL SPIN DENSITY") == 0) usable_cube = .false.
1708 ELSE
1709 IF (index(title, "ELECTRON DENSITY") > 0 .AND. &
1710 index(title, "TOTAL ELECTRON DENSITY") == 0) usable_cube = .false.
1711 END IF
1712 END IF
1713 END IF
1714 CALL para_env%bcast(usable_cube, 0)
1715 IF (.NOT. usable_cube) THEN
1716 CALL cp_abort(__location__, &
1717 "GAPW fragment CDFT cannot use a regular CP2K GAPW density cube. "// &
1718 "Generate the reference with "// &
1719 "E_DENSITY_CUBE% DENSITY_INCLUDE TOTAL_DENSITY. "// &
1720 "Invalid reference: "//trim(filename))
1721 END IF
1722 END SUBROUTINE verify_gapw_fragment_cube
1723
1724 END SUBROUTINE prepare_fragment_constraint
1725
1726END MODULE qs_cdft_methods
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:126
cell_types::pbc
Definition cell_types.F:105
kahan_sum::accurate_dot_product
Definition kahan_sum.F:52
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_set
Definition pw_methods.F:127
pw_methods::pw_zero
Definition pw_methods.F:83
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:208
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_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
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:838
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:1058
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:119
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:27
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:116
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:18
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition input_constants.F:821
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition input_constants.F:821
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition input_constants.F:821
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition input_constants.F:773
input_constants::outer_scf_hirshfeld_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Definition input_constants.F:773
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition input_constants.F:821
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
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:27
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:102
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:589
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:1117
qs_cdft_utils::hirshfeld_constraint_init
subroutine, public hirshfeld_constraint_init(qs_env)
Initializes Gaussian Hirshfeld constraints.
Definition qs_cdft_utils.F:987
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:869
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:1329
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, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_force_types
Definition qs_force_types.F:13
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:137
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:478
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition realspace_grid_types.F:652
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:2020
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1913
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
hirshfeld_types::hirshfeld_type
quantities needed for a Hirshfeld based partitioning of real space
Definition hirshfeld_types.F:36
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_cdft_types::becke_constraint_type
Definition qs_cdft_types.F:76
qs_cdft_types::cdft_control_type
Definition qs_cdft_types.F:218
qs_cdft_types::cdft_group_type
control parameters for CDFT simulations
Definition qs_cdft_types.F:185
qs_cdft_types::hirshfeld_constraint_type
Definition qs_cdft_types.F:110
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56
realspace_grid_types::realspace_grid_desc_type
Definition realspace_grid_types.F:90
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136