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