(git:5e7fe52)
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
143 INTEGER :: handle, i, iatom, igroup, ind(3), ip, j, &
144 jatom, jp, k, natom, np(3), nskipped
145 INTEGER, ALLOCATABLE, DIMENSION(:) :: catom
146 INTEGER, DIMENSION(2, 3) :: bo, bo_conf
147 LOGICAL :: in_memory, my_just_gradients
148 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint, skip_me
149 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: atom_in_group
150 REAL(kind=dp) :: dist1, dist2, dmyexp, dvol, eps_cavity, &
151 my1, my1_homo, myexp, sum_cell_f_all, &
152 th, tmp_const
153 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cell_functions, ds_dr_i, ds_dr_j, &
154 sum_cell_f_group
155 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_sum_pm_dr, dp_i_dri
156 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dp_i_drj
157 REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dmy_dr_i, dmy_dr_j, &
158 dr, dr1_r2, dr_i_dr, dr_ij_dr, &
159 dr_j_dr, grid_p, r, r1, shift
160 REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
161 TYPE(becke_constraint_type), POINTER :: becke_control
162 TYPE(cdft_control_type), POINTER :: cdft_control
163 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
164 TYPE(cell_type), POINTER :: cell
165 TYPE(dft_control_type), POINTER :: dft_control
166 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
167 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge
168
169 NULLIFY (cutoffs, cell, dft_control, particle_set, group, charge, cdft_control)
170 CALL timeset(routinen, handle)
171 ! Get simulation environment
172 CALL get_qs_env(qs_env, &
173 cell=cell, &
174 particle_set=particle_set, &
175 natom=natom, &
176 dft_control=dft_control)
177 cdft_control => dft_control%qs_control%cdft_control
178 becke_control => cdft_control%becke_control
179 group => cdft_control%group
180 cutoffs => becke_control%cutoffs
181 IF (cdft_control%atomic_charges) THEN
182 charge => cdft_control%charge
183 END IF
184 in_memory = .false.
185 IF (cdft_control%save_pot) THEN
186 in_memory = becke_control%in_memory
187 END IF
188 eps_cavity = becke_control%eps_cavity
189 ! Decide if only gradients need to be calculated
190 my_just_gradients = .false.
191 IF (PRESENT(just_gradients)) my_just_gradients = just_gradients
192 IF (my_just_gradients) THEN
193 in_memory = .true.
194 ! Pairwise distances need to be recalculated
195 IF (becke_control%vector_buffer%store_vectors) THEN
196 ALLOCATE (becke_control%vector_buffer%distances(natom))
197 ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom))
198 IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom))
199 ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom))
200 END IF
201 ALLOCATE (becke_control%vector_buffer%R12(natom, natom))
202 DO i = 1, 3
203 cell_v(i) = cell%hmat(i, i)
204 END DO
205 DO iatom = 1, natom - 1
206 DO jatom = iatom + 1, natom
207 r = particle_set(iatom)%r
208 r1 = particle_set(jatom)%r
209 DO i = 1, 3
210 r(i) = modulo(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
211 r1(i) = modulo(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
212 END DO
213 dist_vec = (r - r1) - anint((r - r1)/cell_v)*cell_v
214 IF (becke_control%vector_buffer%store_vectors) THEN
215 becke_control%vector_buffer%position_vecs(:, iatom) = r(:)
216 IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:)
217 IF (in_memory) THEN
218 becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
219 becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
220 END IF
221 END IF
222 becke_control%vector_buffer%R12(iatom, jatom) = norm2(dist_vec)
223 becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom)
224 END DO
225 END DO
226 END IF
227 ALLOCATE (catom(cdft_control%natoms))
228 IF (cdft_control%save_pot .OR. &
229 becke_control%cavity_confine .OR. &
230 becke_control%should_skip) THEN
231 ALLOCATE (is_constraint(natom))
232 is_constraint = .false.
233 END IF
234 ! This boolean is needed to prevent calculation of atom pairs ji when the pair ij has
235 ! already been calculated (data for pair ji is set using symmetry)
236 ! With gradient precomputation, symmetry exploited for both weight function and gradients
237 ALLOCATE (skip_me(natom))
238 DO i = 1, cdft_control%natoms
239 catom(i) = cdft_control%atoms(i)
240 ! Notice that here is_constraint=.TRUE. also for dummy atoms to properly compute their Becke charges
241 ! A subsequent check (atom_in_group) ensures that the gradients of these dummy atoms are correct
242 IF (cdft_control%save_pot .OR. &
243 becke_control%cavity_confine .OR. &
244 becke_control%should_skip) THEN
245 is_constraint(catom(i)) = .true.
246 END IF
247 END DO
248 bo = group(1)%weight%pw_grid%bounds_local
249 dvol = group(1)%weight%pw_grid%dvol
250 dr = group(1)%weight%pw_grid%dr
251 np = group(1)%weight%pw_grid%npts
252 shift = -real(modulo(np, 2), dp)*dr/2.0_dp
253 DO i = 1, 3
254 cell_v(i) = cell%hmat(i, i)
255 END DO
256 ! If requested, allocate storage for gradients
257 IF (in_memory) THEN
258 bo_conf = bo
259 ! With confinement active, we dont need to store gradients outside
260 ! the confinement bounds since they vanish for all particles
261 IF (becke_control%cavity_confine) THEN
262 bo_conf(1, 3) = becke_control%confine_bounds(1)
263 bo_conf(2, 3) = becke_control%confine_bounds(2)
264 END IF
265 ALLOCATE (atom_in_group(SIZE(group), natom))
266 atom_in_group = .false.
267 DO igroup = 1, SIZE(group)
268 ALLOCATE (group(igroup)%gradients(3*natom, bo_conf(1, 1):bo_conf(2, 1), &
269 bo_conf(1, 2):bo_conf(2, 2), &
270 bo_conf(1, 3):bo_conf(2, 3)))
271 group(igroup)%gradients = 0.0_dp
272 ALLOCATE (group(igroup)%d_sum_const_dR(3, natom))
273 group(igroup)%d_sum_const_dR = 0.0_dp
274 DO ip = 1, SIZE(group(igroup)%atoms)
275 atom_in_group(igroup, group(igroup)%atoms(ip)) = .true.
276 END DO
277 END DO
278 END IF
279 ! Allocate remaining work
280 ALLOCATE (sum_cell_f_group(SIZE(group)))
281 ALLOCATE (cell_functions(natom))
282 IF (in_memory) THEN
283 ALLOCATE (ds_dr_j(3))
284 ALLOCATE (ds_dr_i(3))
285 ALLOCATE (d_sum_pm_dr(3, natom))
286 ALLOCATE (dp_i_drj(3, natom, natom))
287 ALLOCATE (dp_i_dri(3, natom))
288 th = 1.0e-8_dp
289 END IF
290 ! Build constraint
291 DO k = bo(1, 1), bo(2, 1)
292 DO j = bo(1, 2), bo(2, 2)
293 DO i = bo(1, 3), bo(2, 3)
294 ! If the grid point is too far from all constraint atoms and cavity confinement is active,
295 ! we can skip this grid point as it does not contribute to the weight or gradients
296 IF (becke_control%cavity_confine) THEN
297 IF (becke_control%cavity%array(k, j, i) < eps_cavity) cycle
298 END IF
299 ind = [k, j, i]
300 grid_p(1) = k*dr(1) + shift(1)
301 grid_p(2) = j*dr(2) + shift(2)
302 grid_p(3) = i*dr(3) + shift(3)
303 nskipped = 0
304 cell_functions = 1.0_dp
305 skip_me = .false.
306 IF (becke_control%vector_buffer%store_vectors) becke_control%vector_buffer%distances = 0.0_dp
307 IF (in_memory) THEN
308 d_sum_pm_dr = 0.0_dp
309 DO igroup = 1, SIZE(group)
310 group(igroup)%d_sum_const_dR = 0.0_dp
311 END DO
312 dp_i_dri = 0.0_dp
313 END IF
314 ! Iterate over all atoms in the system
315 DO iatom = 1, natom
316 IF (skip_me(iatom)) THEN
317 cell_functions(iatom) = 0.0_dp
318 IF (becke_control%should_skip) THEN
319 IF (is_constraint(iatom)) nskipped = nskipped + 1
320 IF (nskipped == cdft_control%natoms) THEN
321 IF (in_memory) THEN
322 IF (becke_control%cavity_confine) THEN
323 becke_control%cavity%array(k, j, i) = 0.0_dp
324 END IF
325 END IF
326 EXIT
327 END IF
328 END IF
329 cycle
330 END IF
331 IF (becke_control%vector_buffer%store_vectors) THEN
332 IF (becke_control%vector_buffer%distances(iatom) == 0.0_dp) THEN
333 r = becke_control%vector_buffer%position_vecs(:, iatom)
334 dist_vec = (r - grid_p) - anint((r - grid_p)/cell_v)*cell_v
335 dist1 = norm2(dist_vec)
336 becke_control%vector_buffer%distance_vecs(:, iatom) = dist_vec
337 becke_control%vector_buffer%distances(iatom) = dist1
338 ELSE
339 dist_vec = becke_control%vector_buffer%distance_vecs(:, iatom)
340 dist1 = becke_control%vector_buffer%distances(iatom)
341 END IF
342 ELSE
343 r = particle_set(iatom)%r
344 DO ip = 1, 3
345 r(ip) = modulo(r(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
346 END DO
347 dist_vec = (r - grid_p) - anint((r - grid_p)/cell_v)*cell_v
348 dist1 = norm2(dist_vec)
349 END IF
350 IF (dist1 <= cutoffs(iatom)) THEN
351 IF (in_memory) THEN
352 IF (dist1 <= th) dist1 = th
353 dr_i_dr(:) = dist_vec(:)/dist1
354 END IF
355 DO jatom = 1, natom
356 IF (jatom /= iatom) THEN
357 ! Using pairwise symmetry, execute block only for such j<i
358 ! that have previously not been looped over
359 ! Note that if skip_me(jatom) = .TRUE., this means that the outer
360 ! loop over iatom skipped this index when iatom=jatom, but we still
361 ! need to compute the pair for iatom>jatom
362 IF (jatom < iatom) THEN
363 IF (.NOT. skip_me(jatom)) cycle
364 END IF
365 IF (becke_control%vector_buffer%store_vectors) THEN
366 IF (becke_control%vector_buffer%distances(jatom) == 0.0_dp) THEN
367 r1 = becke_control%vector_buffer%position_vecs(:, jatom)
368 dist_vec = (r1 - grid_p) - anint((r1 - grid_p)/cell_v)*cell_v
369 dist2 = norm2(dist_vec)
370 becke_control%vector_buffer%distance_vecs(:, jatom) = dist_vec
371 becke_control%vector_buffer%distances(jatom) = dist2
372 ELSE
373 dist_vec = becke_control%vector_buffer%distance_vecs(:, jatom)
374 dist2 = becke_control%vector_buffer%distances(jatom)
375 END IF
376 ELSE
377 r1 = particle_set(jatom)%r
378 DO ip = 1, 3
379 r1(ip) = modulo(r1(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
380 END DO
381 dist_vec = (r1 - grid_p) - anint((r1 - grid_p)/cell_v)*cell_v
382 dist2 = norm2(dist_vec)
383 END IF
384 IF (in_memory) THEN
385 IF (becke_control%vector_buffer%store_vectors) THEN
386 dr1_r2 = becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom)
387 ELSE
388 dr1_r2 = (r - r1) - anint((r - r1)/cell_v)*cell_v
389 END IF
390 IF (dist2 <= th) dist2 = th
391 tmp_const = (becke_control%vector_buffer%R12(iatom, jatom)**3)
392 dr_ij_dr(:) = dr1_r2(:)/tmp_const
393 !derivative w.r.t. Rj
394 dr_j_dr = dist_vec(:)/dist2
395 dmy_dr_j(:) = -(dr_j_dr(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dr(:))
396 !derivative w.r.t. Ri
397 dmy_dr_i(:) = dr_i_dr(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dr(:)
398 END IF
399 ! myij
400 my1 = (dist1 - dist2)/becke_control%vector_buffer%R12(iatom, jatom)
401 IF (becke_control%adjust) THEN
402 my1_homo = my1 ! Homonuclear quantity needed for gradient
403 my1 = my1 + becke_control%aij(iatom, jatom)*(1.0_dp - my1**2)
404 END IF
405 ! f(myij)
406 myexp = 1.5_dp*my1 - 0.5_dp*my1**3
407 IF (in_memory) THEN
408 dmyexp = 1.5_dp - 1.5_dp*my1**2
409 tmp_const = (1.5_dp**2)*dmyexp*(1 - myexp**2)* &
410 (1.0_dp - ((1.5_dp*myexp - 0.5_dp*(myexp**3))**2))
411 ! d s(myij)/d R_i
412 ds_dr_i(:) = -0.5_dp*tmp_const*dmy_dr_i(:)
413 ! d s(myij)/d R_j
414 ds_dr_j(:) = -0.5_dp*tmp_const*dmy_dr_j(:)
415 IF (becke_control%adjust) THEN
416 tmp_const = 1.0_dp - 2.0_dp*my1_homo* &
417 becke_control%aij(iatom, jatom)
418 ds_dr_i(:) = ds_dr_i(:)*tmp_const
419 ! tmp_const is same for both since aij=-aji and myij=-myji
420 ds_dr_j(:) = ds_dr_j(:)*tmp_const
421 END IF
422 END IF
423 ! s(myij) = f[f(f{myij})]
424 myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
425 myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
426 tmp_const = 0.5_dp*(1.0_dp - myexp)
427 cell_functions(iatom) = cell_functions(iatom)*tmp_const
428 IF (in_memory) THEN
429 IF (abs(tmp_const) <= th) tmp_const = tmp_const + th
430 ! P_i independent part of dP_i/dR_i
431 dp_i_dri(:, iatom) = dp_i_dri(:, iatom) + ds_dr_i(:)/tmp_const
432 ! P_i independent part of dP_i/dR_j
433 dp_i_drj(:, iatom, jatom) = ds_dr_j(:)/tmp_const
434 END IF
435
436 IF (dist2 <= cutoffs(jatom)) THEN
437 tmp_const = 0.5_dp*(1.0_dp + myexp) ! s(myji)
438 cell_functions(jatom) = cell_functions(jatom)*tmp_const
439 IF (in_memory) THEN
440 IF (abs(tmp_const) <= th) tmp_const = tmp_const + th
441 ! P_j independent part of dP_j/dR_i
442 ! d s(myji)/d R_i = -d s(myij)/d R_i
443 dp_i_drj(:, jatom, iatom) = -ds_dr_i(:)/tmp_const
444 ! P_j independent part of dP_j/dR_j
445 ! d s(myji)/d R_j = -d s(myij)/d R_j
446 dp_i_dri(:, jatom) = dp_i_dri(:, jatom) - ds_dr_j(:)/tmp_const
447 END IF
448 ELSE
449 skip_me(jatom) = .true.
450 END IF
451 END IF
452 END DO ! jatom
453 IF (in_memory) THEN
454 ! Final value of dP_i_dRi
455 dp_i_dri(:, iatom) = cell_functions(iatom)*dp_i_dri(:, iatom)
456 ! Update relevant sums with value
457 d_sum_pm_dr(:, iatom) = d_sum_pm_dr(:, iatom) + dp_i_dri(:, iatom)
458 IF (is_constraint(iatom)) THEN
459 DO igroup = 1, SIZE(group)
460 IF (.NOT. atom_in_group(igroup, iatom)) cycle
461 DO jp = 1, SIZE(group(igroup)%atoms)
462 IF (iatom == group(igroup)%atoms(jp)) THEN
463 ip = jp
464 EXIT
465 END IF
466 END DO
467 group(igroup)%d_sum_const_dR(1:3, iatom) = group(igroup)%d_sum_const_dR(1:3, iatom) + &
468 group(igroup)%coeff(ip)*dp_i_dri(:, iatom)
469 END DO
470 END IF
471 DO jatom = 1, natom
472 IF (jatom /= iatom) THEN
473 ! Final value of dP_i_dRj
474 dp_i_drj(:, iatom, jatom) = cell_functions(iatom)*dp_i_drj(:, iatom, jatom)
475 ! Update where needed
476 d_sum_pm_dr(:, jatom) = d_sum_pm_dr(:, jatom) + dp_i_drj(:, iatom, jatom)
477 IF (is_constraint(iatom)) THEN
478 DO igroup = 1, SIZE(group)
479 IF (.NOT. atom_in_group(igroup, iatom)) cycle
480 ip = -1
481 DO jp = 1, SIZE(group(igroup)%atoms)
482 IF (iatom == group(igroup)%atoms(jp)) THEN
483 ip = jp
484 EXIT
485 END IF
486 END DO
487 group(igroup)%d_sum_const_dR(1:3, jatom) = group(igroup)%d_sum_const_dR(1:3, jatom) + &
488 group(igroup)%coeff(ip)* &
489 dp_i_drj(:, iatom, jatom)
490 END DO
491 END IF
492 END IF
493 END DO
494 END IF
495 ELSE
496 cell_functions(iatom) = 0.0_dp
497 skip_me(iatom) = .true.
498 IF (becke_control%should_skip) THEN
499 IF (is_constraint(iatom)) nskipped = nskipped + 1
500 IF (nskipped == cdft_control%natoms) THEN
501 IF (in_memory) THEN
502 IF (becke_control%cavity_confine) THEN
503 becke_control%cavity%array(k, j, i) = 0.0_dp
504 END IF
505 END IF
506 EXIT
507 END IF
508 END IF
509 END IF
510 END DO !iatom
511 IF (nskipped == cdft_control%natoms) cycle
512 ! Sum up cell functions
513 sum_cell_f_group = 0.0_dp
514 DO igroup = 1, SIZE(group)
515 DO ip = 1, SIZE(group(igroup)%atoms)
516 sum_cell_f_group(igroup) = sum_cell_f_group(igroup) + group(igroup)%coeff(ip)* &
517 cell_functions(group(igroup)%atoms(ip))
518 END DO
519 END DO
520 sum_cell_f_all = 0.0_dp
521 DO ip = 1, natom
522 sum_cell_f_all = sum_cell_f_all + cell_functions(ip)
523 END DO
524 ! Gradients at (k,j,i)
525 IF (in_memory .AND. abs(sum_cell_f_all) > 0.0_dp) THEN
526 DO igroup = 1, SIZE(group)
527 DO iatom = 1, natom
528 group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) = &
529 group(igroup)%d_sum_const_dR(1:3, iatom)/sum_cell_f_all - sum_cell_f_group(igroup)* &
530 d_sum_pm_dr(1:3, iatom)/(sum_cell_f_all**2)
531 END DO
532 END DO
533 END IF
534 ! Weight function(s) at (k,j,i)
535 IF (.NOT. my_just_gradients .AND. abs(sum_cell_f_all) > 0.000001) THEN
536 DO igroup = 1, SIZE(group)
537 group(igroup)%weight%array(k, j, i) = sum_cell_f_group(igroup)/sum_cell_f_all
538 END DO
539 IF (cdft_control%atomic_charges) THEN
540 DO iatom = 1, cdft_control%natoms
541 charge(iatom)%array(k, j, i) = cell_functions(catom(iatom))/sum_cell_f_all
542 END DO
543 END IF
544 END IF
545 END DO
546 END DO
547 END DO
548 ! Release storage
549 IF (in_memory) THEN
550 DEALLOCATE (ds_dr_j)
551 DEALLOCATE (ds_dr_i)
552 DEALLOCATE (d_sum_pm_dr)
553 DEALLOCATE (dp_i_drj)
554 DEALLOCATE (dp_i_dri)
555 DO igroup = 1, SIZE(group)
556 DEALLOCATE (group(igroup)%d_sum_const_dR)
557 END DO
558 DEALLOCATE (atom_in_group)
559 IF (becke_control%vector_buffer%store_vectors) THEN
560 DEALLOCATE (becke_control%vector_buffer%pair_dist_vecs)
561 END IF
562 END IF
563 NULLIFY (cutoffs)
564 IF (ALLOCATED(is_constraint)) THEN
565 DEALLOCATE (is_constraint)
566 END IF
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 smooth density contribution. GAPW one-center terms are added
608 ! while the Kohn-Sham matrix is assembled.
609 CALL cdft_constraint_integrate(qs_env)
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, 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_rs, 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 hirshfeld_control%print_density = .false.
687 END IF
688
689 ALLOCATE (coefficients(natom))
690 ALLOCATE (is_constraint(natom))
691
692 subpatch_pattern = 0
693 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
694 radius = 100.0_dp
695
696 lb_pw(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
697 ub_pw(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
698
699 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
700 auxbas_pw_pool=auxbas_pw_pool)
701 CALL rs_grid_create(rs_rho_all, auxbas_rs_desc)
702 CALL rs_grid_zero(rs_rho_all)
703
704 dr_rs(1) = rs_rho_all%desc%dh(1, 1)
705 dr_rs(2) = rs_rho_all%desc%dh(2, 2)
706 dr_rs(3) = rs_rho_all%desc%dh(3, 3)
707 lb_rs(1) = lbound(rs_rho_all%r(:, :, :), 1)
708 lb_rs(2) = lbound(rs_rho_all%r(:, :, :), 2)
709 lb_rs(3) = lbound(rs_rho_all%r(:, :, :), 3)
710 ub_rs(1) = ubound(rs_rho_all%r(:, :, :), 1)
711 ub_rs(2) = ubound(rs_rho_all%r(:, :, :), 2)
712 ub_rs(3) = ubound(rs_rho_all%r(:, :, :), 3)
713
714 ! For each CDFT group
715 DO igroup = 1, SIZE(cdft_control%group)
716
717 IF (igroup == 2 .AND. .NOT. cdft_control%in_memory) THEN
718 CALL rs_grid_zero(rs_rho_all)
719 END IF
720 bo = cdft_control%group(igroup)%weight%pw_grid%bounds_local
721
722 ! Coefficients
723 coefficients(:) = 0.0_dp
724 is_constraint = .false.
725 DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
726 coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
727 is_constraint(cdft_control%group(igroup)%atoms(i)) = .true.
728 END DO
729
730 ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group
731 CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc)
732 CALL rs_grid_zero(rs_rho_constr)
733
734 ! rs_single: Gaussian density over single atoms when required
735 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
736 ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic(cdft_control%natoms))
737 ALLOCATE (rs_single(cdft_control%natoms))
738 DO i = 1, cdft_control%natoms
739 CALL rs_grid_create(rs_single(i), auxbas_rs_desc)
740 CALL rs_grid_zero(rs_single(i))
741 END DO
742 END IF
743
744 ! Setup pw
745 CALL pw_zero(cdft_control%group(igroup)%weight)
746
747 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
748 CALL pw_set(cdft_control%group(igroup)%hw_rho_total_constraint, 1.0_dp)
749
750 IF (igroup == 1) THEN
751 CALL auxbas_pw_pool%create_pw(cdft_control%hw_rho_total)
752 CALL pw_set(cdft_control%hw_rho_total, 1.0_dp)
753
754 IF (hirshfeld_control%print_density) THEN
755 DO iatom = 1, cdft_control%natoms
756 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic(iatom))
757 CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic(iatom), 1.0_dp)
758 END DO
759 END IF
760 END IF
761
762 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
763 ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge(cdft_control%natoms))
764 ALLOCATE (rs_single_charge(cdft_control%natoms))
765 ALLOCATE (compute_charge(natom))
766 compute_charge = .false.
767
768 DO i = 1, cdft_control%natoms
769 CALL rs_grid_create(rs_single_charge(i), auxbas_rs_desc)
770 CALL rs_grid_zero(rs_single_charge(i))
771 compute_charge(cdft_control%atoms(i)) = .true.
772 END DO
773
774 DO iatom = 1, cdft_control%natoms
775 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom))
776 CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom), 1.0_dp)
777 END DO
778 END IF
779
780 ALLOCATE (pab(1, 1))
781 nthread = 1
782 ithread = 0
783
784 DO ikind = 1, SIZE(atomic_kind_set)
785 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
786 IF (numexp <= 0) cycle
787 CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
788 ALLOCATE (cores(num_species))
789
790 DO iex = 1, numexp
791 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
792 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
793 npme = 0
794 cores = 0
795 DO iatom = 1, num_species
796 atom_a = atom_list(iatom)
797 ra(:) = pbc(particle_set(atom_a)%r, cell)
798 IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
799 IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
800 npme = npme + 1
801 cores(npme) = iatom
802 END IF
803 ELSE
804 npme = npme + 1
805 cores(npme) = iatom
806 END IF
807 END DO
808 DO j = 1, npme
809 iatom = cores(j)
810 atom_a = atom_list(iatom)
811 pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
812 ra(:) = pbc(particle_set(atom_a)%r, cell)
813
814 IF (hirshfeld_control%use_atomic_cutoff) THEN
815 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
816 ra=ra, rb=ra, rp=ra, &
817 zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
818 pab=pab, o1=0, o2=0, & ! without map_consistent
819 prefactor=1.0_dp, cutoff=0.0_dp)
820 END IF
821
822 IF (igroup == 1) THEN
823 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
824 [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
825 rs_rho_all, radius=radius, &
826 ga_gb_function=grid_func_ab, use_subpatch=.true., &
827 subpatch_pattern=subpatch_pattern)
828 END IF
829
830 IF (is_constraint(atom_a)) THEN
831 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
832 [0.0_dp, 0.0_dp, 0.0_dp], coefficients(atom_a), &
833 pab, 0, 0, rs_rho_constr, &
834 radius=radius, &
835 ga_gb_function=grid_func_ab, use_subpatch=.true., &
836 subpatch_pattern=subpatch_pattern)
837 END IF
838
839 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
840 IF (is_constraint(atom_a)) THEN
841 DO iatom = 1, cdft_control%natoms
842 IF (atom_a == cdft_control%atoms(iatom)) EXIT
843 END DO
844 cpassert(iatom <= cdft_control%natoms)
845 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
846 [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
847 rs_single(iatom), radius=radius, &
848 ga_gb_function=grid_func_ab, use_subpatch=.true., &
849 subpatch_pattern=subpatch_pattern)
850 END IF
851 END IF
852
853 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
854 IF (compute_charge(atom_a)) THEN
855 DO iatom = 1, cdft_control%natoms
856 IF (atom_a == cdft_control%atoms(iatom)) EXIT
857 END DO
858 cpassert(iatom <= cdft_control%natoms)
859 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
860 [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
861 rs_single_charge(iatom), radius=radius, &
862 ga_gb_function=grid_func_ab, use_subpatch=.true., &
863 subpatch_pattern=subpatch_pattern)
864 END IF
865 END IF
866
867 END DO
868 END DO
869 DEALLOCATE (cores)
870 END DO
871 DEALLOCATE (pab)
872
873 IF (igroup == 1) THEN
874 CALL transfer_rs2pw(rs_rho_all, cdft_control%hw_rho_total)
875 END IF
876
877 CALL transfer_rs2pw(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint)
878 CALL rs_grid_release(rs_rho_constr)
879
880 ! Calculate weight function
881 CALL hfun_scale(cdft_control%group(igroup)%weight%array, &
882 cdft_control%group(igroup)%hw_rho_total_constraint%array, &
883 cdft_control%hw_rho_total%array, divide=.true., &
884 small=hirshfeld_control%eps_cutoff)
885
886 ! Calculate charges
887 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
888 DO i = 1, cdft_control%natoms
889 CALL transfer_rs2pw(rs_single_charge(i), cdft_control%group(igroup)%hw_rho_atomic_charge(i))
890 CALL hfun_scale(cdft_control%charge(i)%array, &
891 cdft_control%group(igroup)%hw_rho_atomic_charge(i)%array, &
892 cdft_control%hw_rho_total%array, divide=.true., &
893 small=hirshfeld_control%eps_cutoff)
894 END DO
895 END IF
896
897 ! Print atomic densities if requested
898 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
899 DO i = 1, cdft_control%natoms
900 CALL transfer_rs2pw(rs_single(i), cdft_control%group(igroup)%hw_rho_atomic(i))
901 END DO
902 CALL cdft_print_hirshfeld_density(qs_env)
903 END IF
904
905 END DO
906
907 DO igroup = 1, SIZE(cdft_control%group)
908
909 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
910
911 IF (.NOT. cdft_control%in_memory .AND. igroup == 1) THEN
912 CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
913 END IF
914
915 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
916 DO i = 1, cdft_control%natoms
917 CALL rs_grid_release(rs_single(i))
918 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic(i))
919 END DO
920 DEALLOCATE (rs_single)
921 DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic)
922 END IF
923
924 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
925 DO i = 1, cdft_control%natoms
926 CALL rs_grid_release(rs_single_charge(i))
927 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(i))
928 END DO
929 DEALLOCATE (rs_single_charge)
930 DEALLOCATE (compute_charge)
931 DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge)
932 END IF
933
934 END DO
935
936 IF (cdft_control%in_memory) THEN
937 DO igroup = 1, SIZE(cdft_control%group)
938 ALLOCATE (cdft_control%group(igroup)%gradients_x(1*natom, lb_pw(1):ub_pw(1), &
939 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
940 cdft_control%group(igroup)%gradients_x(:, :, :, :) = 0.0_dp
941 END DO
942 END IF
943
944 IF (cdft_control%in_memory) THEN
945 DO igroup = 1, SIZE(cdft_control%group)
946
947 ALLOCATE (pab(1, 1))
948 nthread = 1
949 ithread = 0
950 atoms_memory = hirshfeld_control%atoms_memory
951
952 DO ikind = 1, SIZE(atomic_kind_set)
953 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
954 IF (numexp <= 0) cycle
955 CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
956
957 ALLOCATE (pw_single_dr(num_species))
958 ALLOCATE (rs_single_dr(num_species))
959
960 DO i = 1, num_species
961 CALL auxbas_pw_pool%create_pw(pw_single_dr(i))
962 CALL pw_zero(pw_single_dr(i))
963 END DO
964
965 atoms_memory_num = SIZE([(j, j=1, num_species, atoms_memory)])
966
967 ! Can't store all pw grids, therefore split into groups of size atom_memory
968 ! Ideally this code should be re-written to be more memory efficient
969 IF (num_species > atoms_memory) THEN
970 ALLOCATE (num_species_small(atoms_memory_num + 1))
971 num_species_small(1:atoms_memory_num) = [(j, j=1, num_species, atoms_memory)]
972 num_species_small(atoms_memory_num + 1) = num_species
973 ELSE
974 ALLOCATE (num_species_small(2))
975 num_species_small(:) = [1, num_species]
976 END IF
977
978 DO k = 1, SIZE(num_species_small) - 1
979 IF (num_species > atoms_memory) THEN
980 ALLOCATE (cores(num_species_small(k + 1) - (num_species_small(k) - 1)))
981 ELSE
982 ALLOCATE (cores(num_species))
983 END IF
984
985 DO i = num_species_small(k), num_species_small(k + 1)
986 CALL rs_grid_create(rs_single_dr(i), auxbas_rs_desc)
987 CALL rs_grid_zero(rs_single_dr(i))
988 END DO
989 DO iex = 1, numexp
990
991 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
992 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
993 prefactor = 2.0_dp*alpha
994 npme = 0
995 cores = 0
996
997 DO iatom = 1, SIZE(cores)
998 atom_a = atom_list(iatom + (num_species_small(k) - 1))
999 ra(:) = pbc(particle_set(atom_a)%r, cell)
1000
1001 IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
1002 IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
1003 npme = npme + 1
1004 cores(npme) = iatom
1005 END IF
1006 ELSE
1007 npme = npme + 1
1008 cores(npme) = iatom
1009 END IF
1010 END DO
1011 DO j = 1, npme
1012 iatom = cores(j)
1013 atom_a = atom_list(iatom + (num_species_small(k) - 1))
1014 pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
1015 ra(:) = pbc(particle_set(atom_a)%r, cell)
1016 subpatch_pattern = 0
1017
1018 ! Calculate cutoff
1019 IF (hirshfeld_control%use_atomic_cutoff) THEN
1020 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
1021 ra=ra, rb=ra, rp=ra, &
1022 zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
1023 pab=pab, o1=0, o2=0, & ! without map_consistent
1024 prefactor=1.0_dp, cutoff=0.0_dp)
1025 END IF
1026
1027 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1028 [0.0_dp, 0.0_dp, 0.0_dp], prefactor, &
1029 pab, 0, 0, rs_single_dr(iatom + (num_species_small(k) - 1)), &
1030 radius=radius, &
1031 ga_gb_function=grid_func_ab, use_subpatch=.true., &
1032 subpatch_pattern=subpatch_pattern)
1033
1034 END DO
1035 END DO
1036
1037 DO iatom = num_species_small(k), num_species_small(k + 1)
1038 CALL transfer_rs2pw(rs_single_dr(iatom), pw_single_dr(iatom))
1039 CALL rs_grid_release(rs_single_dr(iatom))
1040 END DO
1041
1042 DEALLOCATE (cores)
1043 END DO
1044
1045 DO iatom = 1, num_species
1046 atom_a = atom_list(iatom)
1047 cdft_control%group(igroup)%gradients_x(atom_a, :, :, :) = pw_single_dr(iatom)%array(:, :, :)
1048 CALL auxbas_pw_pool%give_back_pw(pw_single_dr(iatom))
1049 END DO
1050
1051 DEALLOCATE (rs_single_dr)
1052 DEALLOCATE (num_species_small)
1053 DEALLOCATE (pw_single_dr)
1054 END DO
1055 DEALLOCATE (pab)
1056 END DO
1057 END IF
1058
1059 IF (cdft_control%in_memory) THEN
1060 DO igroup = 1, SIZE(cdft_control%group)
1061 ALLOCATE (cdft_control%group(igroup)%gradients_y(1*num_atoms, lb_pw(1):ub_pw(1), &
1062 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1063 ALLOCATE (cdft_control%group(igroup)%gradients_z(1*num_atoms, lb_pw(1):ub_pw(1), &
1064 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1065 cdft_control%group(igroup)%gradients_y(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1066 cdft_control%group(igroup)%gradients_z(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1067 END DO
1068 END IF
1069
1070 ! Calculate gradient if requested
1071 IF (cdft_control%in_memory) THEN
1072
1073 DO igroup = 1, SIZE(cdft_control%group)
1074
1075 ! Coefficients
1076 coefficients(:) = 0.0_dp
1077 is_constraint = .false.
1078 DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
1079 coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
1080 is_constraint(cdft_control%group(igroup)%atoms(i)) = .true.
1081 END DO
1082
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%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%hw_rho_total%array(i, j, k)
1095
1096 r2 = real(i - rho_r(1)%pw_grid%bounds(1, 1), dp)* &
1097 rho_r(1)%pw_grid%dh(:, 1) + &
1098 REAL(j - rho_r(1)%pw_grid%bounds(1, 2), dp)* &
1099 rho_r(1)%pw_grid%dh(:, 2) + &
1100 REAL(k - rho_r(1)%pw_grid%bounds(1, 3), dp)* &
1101 rho_r(1)%pw_grid%dh(:, 3)
1102 r_pbc = pbc(ra, r2, cell)
1103
1104 ! Store gradient d/dR_x w, including term: (r_x - R_x)
1105 cdft_control%group(igroup)%gradients_x(iatom, i, j, k) = &
1106 cdft_control%group(igroup)%gradients_x(iatom, i, j, k)* &
1107 r_pbc(1)*exp_eval
1108
1109 ! Store gradient d/dR_y w, including term: (r_y - R_y)
1110 cdft_control%group(igroup)%gradients_y(iatom, i, j, k) = &
1111 cdft_control%group(igroup)%gradients_y(iatom, i, j, k)* &
1112 r_pbc(2)*exp_eval
1113
1114 ! Store gradient d/dR_z w, including term:(r_z - R_z)
1115 cdft_control%group(igroup)%gradients_z(iatom, i, j, k) = &
1116 cdft_control%group(igroup)%gradients_z(iatom, i, j, k)* &
1117 r_pbc(3)*exp_eval
1118
1119 END IF
1120 END DO
1121 END DO
1122 END DO
1123 END DO
1124 END DO
1125 CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
1126 END IF
1127
1128 CALL rs_grid_release(rs_rho_all)
1129
1130 IF (ALLOCATED(coefficients)) DEALLOCATE (coefficients)
1131 IF (ALLOCATED(is_constraint)) DEALLOCATE (is_constraint)
1132
1133 CALL timestop(handle)
1134
1135 END SUBROUTINE hirshfeld_constraint_low
1136
1137! **************************************************************************************************
1138!> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT
1139!> weight function and the realspace electron density
1140!> \param qs_env ...
1141! **************************************************************************************************
1142 SUBROUTINE cdft_constraint_integrate(qs_env)
1143 TYPE(qs_environment_type), POINTER :: qs_env
1144
1145 CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_integrate'
1146
1147 INTEGER :: handle, i, iatom, igroup, ivar, iw, nvar
1148 LOGICAL :: is_becke
1149 REAL(kind=dp) :: dvol, eps_cavity, sign
1150 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: de, strength, target_val
1151 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge
1152 TYPE(becke_constraint_type), POINTER :: becke_control
1153 TYPE(cdft_control_type), POINTER :: cdft_control
1154 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1155 TYPE(cp_logger_type), POINTER :: logger
1156 TYPE(dft_control_type), POINTER :: dft_control
1157 TYPE(mp_para_env_type), POINTER :: para_env
1158 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge, rho_r
1159 TYPE(qs_energy_type), POINTER :: energy
1160 TYPE(qs_rho_type), POINTER :: rho
1161 TYPE(section_vals_type), POINTER :: cdft_constraint_section
1162
1163 NULLIFY (para_env, dft_control, rho_r, energy, rho, &
1164 logger, cdft_constraint_section, group, charge)
1165 CALL timeset(routinen, handle)
1166 logger => cp_get_default_logger()
1167 CALL get_qs_env(qs_env, &
1168 rho=rho, &
1169 dft_control=dft_control, &
1170 para_env=para_env)
1171 CALL qs_rho_get(rho, rho_r=rho_r)
1172 cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1173 iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
1174 cdft_control => dft_control%qs_control%cdft_control
1175 is_becke = (cdft_control%type == outer_scf_becke_constraint)
1176 becke_control => cdft_control%becke_control
1177 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) THEN
1178 cpabort("Becke control has not been allocated.")
1179 END IF
1180 group => cdft_control%group
1181 ! Initialize
1182 nvar = SIZE(cdft_control%target)
1183 ALLOCATE (strength(nvar))
1184 ALLOCATE (target_val(nvar))
1185 ALLOCATE (de(nvar))
1186 strength(:) = cdft_control%strength(:)
1187 target_val(:) = cdft_control%target(:)
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 cdft_control%value(:) = de(:)
1247 energy%cdft = 0.0_dp
1248 DO ivar = 1, nvar
1249 energy%cdft = energy%cdft + (de(ivar) - target_val(ivar))*strength(ivar)
1250 END DO
1251 ! Print constraint info and atomic CDFT charges
1252 IF (.NOT. dft_control%qs_control%gapw) THEN
1253 CALL cdft_constraint_print(qs_env, electronic_charge)
1254 END IF
1255 ! Deallocate tmp storage
1256 DEALLOCATE (de, strength, target_val)
1257 IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge)
1258 CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO")
1259 CALL timestop(handle)
1260
1261 END SUBROUTINE cdft_constraint_integrate
1262
1263! **************************************************************************************************
1264!> \brief Calculates atomic forces due to a CDFT constraint (Becke or Hirshfeld)
1265!> \param qs_env ...
1266! **************************************************************************************************
1267 SUBROUTINE cdft_constraint_force(qs_env)
1268 TYPE(qs_environment_type), POINTER :: qs_env
1269
1270 CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_force'
1271
1272 INTEGER :: handle, i, iatom, igroup, ikind, ispin, &
1273 j, k, natom, nvar
1274 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
1275 INTEGER, DIMENSION(2, 3) :: bo
1276 INTEGER, DIMENSION(3) :: lb, ub
1277 REAL(kind=dp) :: dvol, eps_cavity, sign
1278 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength
1279 REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
1280 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1281 TYPE(becke_constraint_type), POINTER :: becke_control
1282 TYPE(cdft_control_type), POINTER :: cdft_control
1283 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1284 TYPE(cell_type), POINTER :: cell
1285 TYPE(dft_control_type), POINTER :: dft_control
1286 TYPE(mp_para_env_type), POINTER :: para_env
1287 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1288 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1289 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1290 TYPE(qs_rho_type), POINTER :: rho
1291
1292 CALL timeset(routinen, handle)
1293 NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, &
1294 rho, rho_r, force, cutoffs, becke_control, group)
1295
1296 CALL get_qs_env(qs_env, &
1297 atomic_kind_set=atomic_kind_set, &
1298 natom=natom, &
1299 particle_set=particle_set, &
1300 cell=cell, &
1301 rho=rho, &
1302 force=force, &
1303 dft_control=dft_control, &
1304 para_env=para_env)
1305 CALL qs_rho_get(rho, rho_r=rho_r)
1306
1307 cdft_control => dft_control%qs_control%cdft_control
1308 becke_control => cdft_control%becke_control
1309 group => cdft_control%group
1310 nvar = SIZE(cdft_control%target)
1311 ALLOCATE (strength(nvar))
1312 strength(:) = cdft_control%strength(:)
1313 cutoffs => cdft_control%becke_control%cutoffs
1314 eps_cavity = cdft_control%becke_control%eps_cavity
1315
1316 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1317 atom_of_kind=atom_of_kind, &
1318 kind_of=kind_of)
1319 DO igroup = 1, SIZE(cdft_control%group)
1320 ALLOCATE (cdft_control%group(igroup)%integrated(3, natom))
1321 cdft_control%group(igroup)%integrated = 0.0_dp
1322 END DO
1323
1324 lb(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
1325 ub(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
1326 bo = cdft_control%group(1)%weight%pw_grid%bounds_local
1327 dvol = cdft_control%group(1)%weight%pw_grid%dvol
1328 sign = 1.0_dp
1329
1330 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1331 IF (.NOT. cdft_control%becke_control%in_memory) THEN
1332 CALL becke_constraint_low(qs_env, just_gradients=.true.)
1333 END IF
1334
1335 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1336 IF (.NOT. cdft_control%in_memory) THEN
1337 CALL hirshfeld_constraint_low(qs_env, just_gradients=.true.)
1338 END IF
1339 END IF
1340
1341 ! If no Becke Gaussian confinement
1342 IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN
1343 ! No external control
1344 DO k = bo(1, 1), bo(2, 1)
1345 DO j = bo(1, 2), bo(2, 2)
1346 DO i = bo(1, 3), bo(2, 3)
1347 ! First check if this grid point should be skipped
1348 IF (cdft_control%becke_control%cavity_confine) THEN
1349 IF (cdft_control%becke_control%cavity%array(k, j, i) < eps_cavity) cycle
1350 END IF
1351
1352 DO igroup = 1, SIZE(cdft_control%group)
1353 DO iatom = 1, natom
1354 DO ispin = 1, dft_control%nspins
1355
1356 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1357 CASE (cdft_charge_constraint)
1358 sign = 1.0_dp
1359 CASE (cdft_magnetization_constraint)
1360 IF (ispin == 1) THEN
1361 sign = 1.0_dp
1362 ELSE
1363 sign = -1.0_dp
1364 END IF
1365 CASE (cdft_alpha_constraint)
1366 sign = 1.0_dp
1367 IF (ispin == 2) cycle
1368 CASE (cdft_beta_constraint)
1369 sign = 1.0_dp
1370 IF (ispin == 1) cycle
1371 CASE DEFAULT
1372 cpabort("Unknown constraint type.")
1373 END SELECT
1374
1375 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1376
1377 cdft_control%group(igroup)%integrated(:, iatom) = &
1378 cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1379 cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1380 *rho_r(ispin)%array(k, j, i) &
1381 *dvol
1382
1383 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1384
1385 cdft_control%group(igroup)%integrated(1, iatom) = &
1386 cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1387 cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1388 *rho_r(ispin)%array(k, j, i) &
1389 *dvol
1390
1391 cdft_control%group(igroup)%integrated(2, iatom) = &
1392 cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1393 cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1394 *rho_r(ispin)%array(k, j, i) &
1395 *dvol
1396
1397 cdft_control%group(igroup)%integrated(3, iatom) = &
1398 cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1399 cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1400 *rho_r(ispin)%array(k, j, i) &
1401 *dvol
1402
1403 END IF
1404
1405 END DO
1406 END DO
1407 END DO
1408 END DO
1409 END DO
1410 END DO
1411
1412 ! If Becke Gaussian confinement
1413 ELSE
1414 DO k = lbound(cdft_control%becke_control%cavity_mat, 1), ubound(cdft_control%becke_control%cavity_mat, 1)
1415 DO j = lbound(cdft_control%becke_control%cavity_mat, 2), ubound(cdft_control%becke_control%cavity_mat, 2)
1416 DO i = lbound(cdft_control%becke_control%cavity_mat, 3), ubound(cdft_control%becke_control%cavity_mat, 3)
1417
1418 ! First check if this grid point should be skipped
1419 IF (cdft_control%becke_control%cavity_mat(k, j, i) < eps_cavity) cycle
1420
1421 DO igroup = 1, SIZE(group)
1422 DO iatom = 1, natom
1423 DO ispin = 1, dft_control%nspins
1424 SELECT CASE (group(igroup)%constraint_type)
1425 CASE (cdft_charge_constraint)
1426 sign = 1.0_dp
1427 CASE (cdft_magnetization_constraint)
1428 IF (ispin == 1) THEN
1429 sign = 1.0_dp
1430 ELSE
1431 sign = -1.0_dp
1432 END IF
1433 CASE (cdft_alpha_constraint)
1434 sign = 1.0_dp
1435 IF (ispin == 2) cycle
1436 CASE (cdft_beta_constraint)
1437 sign = 1.0_dp
1438 IF (ispin == 1) cycle
1439 CASE DEFAULT
1440 cpabort("Unknown constraint type.")
1441 END SELECT
1442
1443 ! Integrate gradient of weight function
1444 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1445
1446 cdft_control%group(igroup)%integrated(:, iatom) = &
1447 cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1448 cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1449 *rho_r(ispin)%array(k, j, i) &
1450 *dvol
1451
1452 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1453
1454 cdft_control%group(igroup)%integrated(1, iatom) = &
1455 cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1456 cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1457 *rho_r(ispin)%array(k, j, i) &
1458 *dvol
1459
1460 cdft_control%group(igroup)%integrated(2, iatom) = &
1461 cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1462 cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1463 *rho_r(ispin)%array(k, j, i) &
1464 *dvol
1465
1466 cdft_control%group(igroup)%integrated(3, iatom) = &
1467 cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1468 cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1469 *rho_r(ispin)%array(k, j, i) &
1470 *dvol
1471
1472 END IF
1473
1474 END DO
1475 END DO
1476 END DO
1477 END DO
1478 END DO
1479 END DO
1480 END IF
1481
1482 IF (.NOT. cdft_control%transfer_pot) THEN
1483 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1484 DO igroup = 1, SIZE(group)
1485 DEALLOCATE (cdft_control%group(igroup)%gradients)
1486 END DO
1487 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1488 DO igroup = 1, SIZE(group)
1489 DEALLOCATE (cdft_control%group(igroup)%gradients_x)
1490 DEALLOCATE (cdft_control%group(igroup)%gradients_y)
1491 DEALLOCATE (cdft_control%group(igroup)%gradients_z)
1492 END DO
1493 END IF
1494 END IF
1495
1496 DO igroup = 1, SIZE(group)
1497 CALL para_env%sum(group(igroup)%integrated)
1498 END DO
1499
1500 ! Update force only on master process. Otherwise force due to constraint becomes multiplied
1501 ! by the number of processes when the final force%rho_elec is constructed in qs_force
1502 ! by mp_summing [the final integrated(:,:) is distributed on all processors]
1503 IF (para_env%is_source()) THEN
1504 DO igroup = 1, SIZE(group)
1505 DO iatom = 1, natom
1506 ikind = kind_of(iatom)
1507 i = atom_of_kind(iatom)
1508 force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup)
1509 END DO
1510 END DO
1511 END IF
1512
1513 DEALLOCATE (strength)
1514 DO igroup = 1, SIZE(group)
1515 DEALLOCATE (group(igroup)%integrated)
1516 END DO
1517 NULLIFY (group)
1518
1519 CALL timestop(handle)
1520
1521 END SUBROUTINE cdft_constraint_force
1522
1523! **************************************************************************************************
1524!> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied
1525!> by the CDFT weight functions and integrated over the realspace grid.
1526!> \param qs_env ...
1527! **************************************************************************************************
1528 SUBROUTINE prepare_fragment_constraint(qs_env)
1529 TYPE(qs_environment_type), POINTER :: qs_env
1530
1531 CHARACTER(len=*), PARAMETER :: routinen = 'prepare_fragment_constraint'
1532
1533 INTEGER :: handle, i, iatom, igroup, natom, &
1534 nelectron_total, nfrag_spins
1535 LOGICAL :: is_becke, needs_spin_density
1536 REAL(kind=dp) :: dvol, multiplier(2), nelectron_frag
1537 TYPE(becke_constraint_type), POINTER :: becke_control
1538 TYPE(cdft_control_type), POINTER :: cdft_control
1539 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1540 TYPE(cp_logger_type), POINTER :: logger
1541 TYPE(dft_control_type), POINTER :: dft_control
1542 TYPE(mp_para_env_type), POINTER :: para_env
1543 TYPE(pw_env_type), POINTER :: pw_env
1544 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1545 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: rho_frag
1546 TYPE(qs_subsys_type), POINTER :: subsys
1547
1548 NULLIFY (para_env, dft_control, logger, subsys, pw_env, auxbas_pw_pool, group)
1549 CALL timeset(routinen, handle)
1550 logger => cp_get_default_logger()
1551 CALL get_qs_env(qs_env, &
1552 natom=natom, &
1553 dft_control=dft_control, &
1554 para_env=para_env)
1555
1556 cdft_control => dft_control%qs_control%cdft_control
1557 is_becke = (cdft_control%type == outer_scf_becke_constraint)
1558 becke_control => cdft_control%becke_control
1559 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) THEN
1560 cpabort("Becke control has not been allocated.")
1561 END IF
1562 group => cdft_control%group
1563 dvol = group(1)%weight%pw_grid%dvol
1564 ! Fragment densities are meaningful only for some calculation types
1565 IF (.NOT. qs_env%single_point_run) THEN
1566 CALL cp_abort(__location__, &
1567 "CDFT fragment constraints are only compatible with single "// &
1568 "point calculations (run_type ENERGY or ENERGY_FORCE).")
1569 END IF
1570 needs_spin_density = .false.
1571 multiplier = 1.0_dp
1572 nfrag_spins = 1
1573 DO igroup = 1, SIZE(group)
1574 SELECT CASE (group(igroup)%constraint_type)
1575 CASE (cdft_charge_constraint)
1576 ! Do nothing
1577 CASE (cdft_magnetization_constraint)
1578 needs_spin_density = .true.
1579 CASE (cdft_alpha_constraint, cdft_beta_constraint)
1580 CALL cp_abort(__location__, &
1581 "CDFT fragment constraint not yet compatible with "// &
1582 "spin specific constraints.")
1583 CASE DEFAULT
1584 cpabort("Unknown constraint type.")
1585 END SELECT
1586 END DO
1587 IF (needs_spin_density) THEN
1588 nfrag_spins = 2
1589 DO i = 1, 2
1590 IF (cdft_control%flip_fragment(i)) multiplier(i) = -1.0_dp
1591 END DO
1592 END IF
1593 ! Read fragment reference densities
1594 ALLOCATE (cdft_control%fragments(nfrag_spins, 2))
1595 ALLOCATE (rho_frag(nfrag_spins))
1596 CALL get_qs_env(qs_env, pw_env=pw_env)
1597 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1598 IF (dft_control%qs_control%gapw) THEN
1599 CALL verify_gapw_fragment_cube(cdft_control%fragment_a_fname, spin_density=.false.)
1600 CALL verify_gapw_fragment_cube(cdft_control%fragment_b_fname, spin_density=.false.)
1601 IF (needs_spin_density) THEN
1602 CALL verify_gapw_fragment_cube(cdft_control%fragment_a_spin_fname, spin_density=.true.)
1603 CALL verify_gapw_fragment_cube(cdft_control%fragment_b_spin_fname, spin_density=.true.)
1604 END IF
1605 END IF
1606 ! Total density (rho_alpha + rho_beta)
1607 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 1))
1608 CALL cp_cube_to_pw(cdft_control%fragments(1, 1), &
1609 cdft_control%fragment_a_fname, 1.0_dp)
1610 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 2))
1611 CALL cp_cube_to_pw(cdft_control%fragments(1, 2), &
1612 cdft_control%fragment_b_fname, 1.0_dp)
1613 ! Spin difference density (rho_alpha - rho_beta) if needed
1614 IF (needs_spin_density) THEN
1615 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 1))
1616 CALL cp_cube_to_pw(cdft_control%fragments(2, 1), &
1617 cdft_control%fragment_a_spin_fname, multiplier(1))
1618 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 2))
1619 CALL cp_cube_to_pw(cdft_control%fragments(2, 2), &
1620 cdft_control%fragment_b_spin_fname, multiplier(2))
1621 END IF
1622 ! Sum up fragments
1623 DO i = 1, nfrag_spins
1624 CALL auxbas_pw_pool%create_pw(rho_frag(i))
1625 CALL pw_copy(cdft_control%fragments(i, 1), rho_frag(i))
1626 CALL pw_axpy(cdft_control%fragments(i, 2), rho_frag(i), 1.0_dp)
1627 CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 1))
1628 CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 2))
1629 END DO
1630 DEALLOCATE (cdft_control%fragments)
1631 ! Check that the number of electrons is consistent
1632 CALL get_qs_env(qs_env, subsys=subsys)
1633 CALL qs_subsys_get(subsys, nelectron_total=nelectron_total)
1634 nelectron_frag = pw_integrate_function(rho_frag(1))
1635 IF (nint(nelectron_frag) /= nelectron_total) THEN
1636 CALL cp_abort(__location__, &
1637 "The number of electrons in the reference and interacting "// &
1638 "configurations does not match. Check your fragment cube files.")
1639 END IF
1640 ! Update constraint target value i.e. perform integration w_i*rho_frag_{tot/spin}*dr
1641 cdft_control%target = 0.0_dp
1642 DO igroup = 1, SIZE(group)
1643 IF (group(igroup)%constraint_type == cdft_charge_constraint) THEN
1644 i = 1
1645 ELSE
1646 i = 2
1647 END IF
1648 IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1649 cdft_control%target(igroup) = cdft_control%target(igroup) + &
1650 accurate_dot_product(group(igroup)%weight%array, rho_frag(i)%array, &
1651 becke_control%cavity_mat, becke_control%eps_cavity)*dvol
1652 ELSE
1653 cdft_control%target(igroup) = cdft_control%target(igroup) + &
1654 pw_integral_ab(group(igroup)%weight, rho_frag(i), local_only=.true.)
1655 END IF
1656 END DO
1657 CALL para_env%sum(cdft_control%target)
1658 ! Calculate reference atomic charges int( w_i * rho_frag * dr )
1659 IF (cdft_control%atomic_charges) THEN
1660 ALLOCATE (cdft_control%charges_fragment(cdft_control%natoms, nfrag_spins))
1661 DO i = 1, nfrag_spins
1662 DO iatom = 1, cdft_control%natoms
1663 cdft_control%charges_fragment(iatom, i) = &
1664 pw_integral_ab(cdft_control%charge(iatom), rho_frag(i), local_only=.true.)
1665 END DO
1666 END DO
1667 CALL para_env%sum(cdft_control%charges_fragment)
1668 END IF
1669 DO i = 1, nfrag_spins
1670 CALL auxbas_pw_pool%give_back_pw(rho_frag(i))
1671 END DO
1672 DEALLOCATE (rho_frag)
1673 cdft_control%fragments_integrated = .true.
1674
1675 CALL timestop(handle)
1676
1677 CONTAINS
1678
1679! **************************************************************************************************
1680!> \brief Require full-density cubes for GAPW fragment references.
1681!> \param filename cube file to inspect
1682!> \param spin_density whether a total spin-density title is expected
1683! **************************************************************************************************
1684 SUBROUTINE verify_gapw_fragment_cube(filename, spin_density)
1685 CHARACTER(LEN=*), INTENT(IN) :: filename
1686 LOGICAL, INTENT(IN) :: spin_density
1687
1688 CHARACTER(LEN=256) :: title
1689 INTEGER :: input_unit, io_status
1690 LOGICAL :: usable_cube
1691
1692 usable_cube = .false.
1693 IF (para_env%is_source()) THEN
1694 CALL open_file(file_name=trim(filename), file_status="OLD", &
1695 file_action="READ", unit_number=input_unit)
1696 READ (input_unit, '(A)', iostat=io_status) title
1697 IF (io_status == 0) THEN
1698 READ (input_unit, '(A)', iostat=io_status) title
1699 usable_cube = io_status == 0
1700 END IF
1701 CALL close_file(input_unit)
1702 IF (usable_cube) THEN
1703 IF (spin_density) THEN
1704 IF (index(title, "SPIN DENSITY") > 0 .AND. &
1705 index(title, "TOTAL SPIN DENSITY") == 0) usable_cube = .false.
1706 ELSE
1707 IF (index(title, "ELECTRON DENSITY") > 0 .AND. &
1708 index(title, "TOTAL ELECTRON DENSITY") == 0) usable_cube = .false.
1709 END IF
1710 END IF
1711 END IF
1712 CALL para_env%bcast(usable_cube, 0)
1713 IF (.NOT. usable_cube) THEN
1714 CALL cp_abort(__location__, &
1715 "GAPW fragment CDFT cannot use a regular CP2K GAPW density cube. "// &
1716 "Generate the reference with "// &
1717 "E_DENSITY_CUBE% DENSITY_INCLUDE TOTAL_DENSITY. "// &
1718 "Invalid reference: "//trim(filename))
1719 END IF
1720 END SUBROUTINE verify_gapw_fragment_cube
1721
1722 END SUBROUTINE prepare_fragment_constraint
1723
1724END 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:125
cell_types::pbc
Definition cell_types.F:104
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: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_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
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:812
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition input_constants.F:812
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition input_constants.F:812
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition input_constants.F:764
input_constants::outer_scf_hirshfeld_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Definition input_constants.F:764
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition input_constants.F:812
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: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: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:539
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:2014
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1907
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:758
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:767
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:219
qs_cdft_types::cdft_group_type
control parameters for CDFT simulations
Definition qs_cdft_types.F:186
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:222
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