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