(git:374b731)
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-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Subroutines for building CDFT constraints
10!> \par History
11!> separated from et_coupling [03.2017]
12!> \author Nico Holmberg [03.2017]
13! **************************************************************************************************
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 cdft_control%atomic_charges = .false.
687 hirshfeld_control%print_density = .false.
688 END IF
689
690 ALLOCATE (coefficients(natom))
691 ALLOCATE (is_constraint(natom))
692
693 subpatch_pattern = 0
694 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
695 radius = 100.0_dp
696
697 dr_pw(1) = rho_r(1)%pw_grid%dr(1)
698 dr_pw(2) = rho_r(1)%pw_grid%dr(2)
699 dr_pw(3) = rho_r(1)%pw_grid%dr(3)
700 lb_pw(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
701 ub_pw(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
702 npts = rho_r(1)%pw_grid%npts
703 origin(1) = (dr_pw(1)*npts(1))*0.5_dp
704 origin(2) = (dr_pw(2)*npts(2))*0.5_dp
705 origin(3) = (dr_pw(3)*npts(3))*0.5_dp
706
707 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
708 auxbas_pw_pool=auxbas_pw_pool)
709 CALL rs_grid_create(rs_rho_all, auxbas_rs_desc)
710 CALL rs_grid_zero(rs_rho_all)
711
712 dr_rs(1) = rs_rho_all%desc%dh(1, 1)
713 dr_rs(2) = rs_rho_all%desc%dh(2, 2)
714 dr_rs(3) = rs_rho_all%desc%dh(3, 3)
715 lb_rs(1) = lbound(rs_rho_all%r(:, :, :), 1)
716 lb_rs(2) = lbound(rs_rho_all%r(:, :, :), 2)
717 lb_rs(3) = lbound(rs_rho_all%r(:, :, :), 3)
718 ub_rs(1) = ubound(rs_rho_all%r(:, :, :), 1)
719 ub_rs(2) = ubound(rs_rho_all%r(:, :, :), 2)
720 ub_rs(3) = ubound(rs_rho_all%r(:, :, :), 3)
721
722 ! For each CDFT group
723 DO igroup = 1, SIZE(cdft_control%group)
724
725 IF (igroup == 2 .AND. .NOT. cdft_control%in_memory) THEN
726 CALL rs_grid_zero(rs_rho_all)
727 END IF
728 bo = cdft_control%group(igroup)%weight%pw_grid%bounds_local
729
730 ! Coefficients
731 coefficients(:) = 0.0_dp
732 is_constraint = .false.
733 DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
734 coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
735 is_constraint(cdft_control%group(igroup)%atoms(i)) = .true.
736 END DO
737
738 ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group
739 CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc)
740 CALL rs_grid_zero(rs_rho_constr)
741
742 ! rs_single: Gaussian density over single atoms when required
743 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
744 ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic(cdft_control%natoms))
745 ALLOCATE (rs_single(cdft_control%natoms))
746 DO i = 1, cdft_control%natoms
747 CALL rs_grid_create(rs_single(i), auxbas_rs_desc)
748 CALL rs_grid_zero(rs_single(i))
749 END DO
750 END IF
751
752 ! Setup pw
753 CALL pw_zero(cdft_control%group(igroup)%weight)
754
755 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
756 CALL pw_set(cdft_control%group(igroup)%hw_rho_total_constraint, 1.0_dp)
757
758 IF (igroup == 1) THEN
759 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total)
760 CALL pw_set(cdft_control%group(igroup)%hw_rho_total, 1.0_dp)
761
762 IF (hirshfeld_control%print_density) THEN
763 DO iatom = 1, cdft_control%natoms
764 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic(iatom))
765 CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic(iatom), 1.0_dp)
766 END DO
767 END IF
768 END IF
769
770 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
771 ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge(cdft_control%natoms))
772 ALLOCATE (rs_single_charge(cdft_control%natoms))
773 ALLOCATE (compute_charge(natom))
774 compute_charge = .false.
775
776 DO i = 1, cdft_control%natoms
777 CALL rs_grid_create(rs_single_charge(i), auxbas_rs_desc)
778 CALL rs_grid_zero(rs_single_charge(i))
779 compute_charge(cdft_control%atoms(i)) = .true.
780 END DO
781
782 DO iatom = 1, cdft_control%natoms
783 CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom))
784 CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom), 1.0_dp)
785 END DO
786 END IF
787
788 ALLOCATE (pab(1, 1))
789 nthread = 1
790 ithread = 0
791
792 DO ikind = 1, SIZE(atomic_kind_set)
793 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
794 IF (numexp <= 0) cycle
795 CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
796 ALLOCATE (cores(num_species))
797
798 DO iex = 1, numexp
799 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
800 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
801 npme = 0
802 cores = 0
803 DO iatom = 1, num_species
804 atom_a = atom_list(iatom)
805 ra(:) = pbc(particle_set(atom_a)%r, cell)
806 IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
807 IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
808 npme = npme + 1
809 cores(npme) = iatom
810 END IF
811 ELSE
812 npme = npme + 1
813 cores(npme) = iatom
814 END IF
815 END DO
816 DO j = 1, npme
817 iatom = cores(j)
818 atom_a = atom_list(iatom)
819 pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
820 ra(:) = pbc(particle_set(atom_a)%r, cell)
821
822 IF (hirshfeld_control%use_atomic_cutoff) THEN
823 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
824 ra=ra, rb=ra, rp=ra, &
825 zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
826 pab=pab, o1=0, o2=0, & ! without map_consistent
827 prefactor=1.0_dp, cutoff=0.0_dp)
828 ! IF (para_env%is_source()) PRINT *, "radius", radius
829 END IF
830
831 IF (igroup == 1) THEN
832 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
833 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
834 rs_rho_all, radius=radius, &
835 ga_gb_function=grid_func_ab, use_subpatch=.true., &
836 subpatch_pattern=subpatch_pattern)
837 END IF
838
839 IF (is_constraint(atom_a)) THEN
840 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
841 (/0.0_dp, 0.0_dp, 0.0_dp/), coefficients(atom_a), &
842 pab, 0, 0, rs_rho_constr, &
843 radius=radius, &
844 ga_gb_function=grid_func_ab, use_subpatch=.true., &
845 subpatch_pattern=subpatch_pattern)
846 END IF
847
848 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
849 IF (is_constraint(atom_a)) THEN
850 DO iatom = 1, cdft_control%natoms
851 IF (atom_a == cdft_control%atoms(iatom)) EXIT
852 END DO
853 cpassert(iatom <= cdft_control%natoms)
854 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
855 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
856 rs_single(iatom), radius=radius, &
857 ga_gb_function=grid_func_ab, use_subpatch=.true., &
858 subpatch_pattern=subpatch_pattern)
859 END IF
860 END IF
861
862 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
863 IF (compute_charge(atom_a)) THEN
864 DO iatom = 1, cdft_control%natoms
865 IF (atom_a == cdft_control%atoms(iatom)) EXIT
866 END DO
867 cpassert(iatom <= cdft_control%natoms)
868 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
869 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
870 rs_single_charge(iatom), radius=radius, &
871 ga_gb_function=grid_func_ab, use_subpatch=.true., &
872 subpatch_pattern=subpatch_pattern)
873 END IF
874 END IF
875
876 END DO
877 END DO
878 DEALLOCATE (cores)
879 END DO
880 DEALLOCATE (pab)
881
882 IF (igroup == 1) THEN
883 CALL transfer_rs2pw(rs_rho_all, cdft_control%group(igroup)%hw_rho_total)
884 END IF
885
886 CALL transfer_rs2pw(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint)
887 CALL rs_grid_release(rs_rho_constr)
888
889 ! Calculate weight function
890 CALL hfun_scale(cdft_control%group(igroup)%weight%array, &
891 cdft_control%group(igroup)%hw_rho_total_constraint%array, &
892 cdft_control%group(1)%hw_rho_total%array, divide=.true., &
893 small=hirshfeld_control%eps_cutoff)
894
895 ! Calculate charges
896 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
897 DO i = 1, cdft_control%natoms
898 CALL transfer_rs2pw(rs_single_charge(i), cdft_control%group(igroup)%hw_rho_atomic_charge(i))
899 CALL hfun_scale(cdft_control%charge(i)%array, &
900 cdft_control%group(igroup)%hw_rho_atomic_charge(i)%array, &
901 cdft_control%group(igroup)%hw_rho_total%array, divide=.true., &
902 small=hirshfeld_control%eps_cutoff)
903 END DO
904 END IF
905
906 ! Print atomic densities if requested
907 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
908 DO i = 1, cdft_control%natoms
909 CALL transfer_rs2pw(rs_single(i), cdft_control%group(igroup)%hw_rho_atomic(i))
910 END DO
911 CALL cdft_print_hirshfeld_density(qs_env)
912 END IF
913
914 END DO
915
916 DO igroup = 1, SIZE(cdft_control%group)
917
918 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
919
920 IF (.NOT. cdft_control%in_memory .AND. igroup == 1) THEN
921 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(1)%hw_rho_total)
922 END IF
923
924 IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
925 DO i = 1, cdft_control%natoms
926 CALL rs_grid_release(rs_single(i))
927 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic(i))
928 END DO
929 DEALLOCATE (rs_single)
930 DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic)
931 END IF
932
933 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
934 DO i = 1, cdft_control%natoms
935 CALL rs_grid_release(rs_single_charge(i))
936 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(i))
937 END DO
938 DEALLOCATE (rs_single_charge)
939 DEALLOCATE (compute_charge)
940 DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge)
941 END IF
942
943 END DO
944
945 IF (cdft_control%in_memory) THEN
946 DO igroup = 1, SIZE(cdft_control%group)
947 ALLOCATE (cdft_control%group(igroup)%gradients_x(1*natom, lb_pw(1):ub_pw(1), &
948 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
949 cdft_control%group(igroup)%gradients_x(:, :, :, :) = 0.0_dp
950 END DO
951 END IF
952
953 IF (cdft_control%in_memory) THEN
954 DO igroup = 1, SIZE(cdft_control%group)
955
956 ALLOCATE (pab(1, 1))
957 nthread = 1
958 ithread = 0
959 atoms_memory = hirshfeld_control%atoms_memory
960
961 DO ikind = 1, SIZE(atomic_kind_set)
962 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
963 IF (numexp <= 0) cycle
964 CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
965
966 ALLOCATE (pw_single_dr(num_species))
967 ALLOCATE (rs_single_dr(num_species))
968
969 DO i = 1, num_species
970 CALL auxbas_pw_pool%create_pw(pw_single_dr(i))
971 CALL pw_zero(pw_single_dr(i))
972 END DO
973
974 atoms_memory_num = SIZE((/(j, j=1, num_species, atoms_memory)/))
975
976 ! Can't store all pw grids, therefore split into groups of size atom_memory
977 ! Ideally this code should be re-written to be more memory efficient
978 IF (num_species > atoms_memory) THEN
979 ALLOCATE (num_species_small(atoms_memory_num + 1))
980 num_species_small(1:atoms_memory_num) = (/(j, j=1, num_species, atoms_memory)/)
981 num_species_small(atoms_memory_num + 1) = num_species
982 ELSE
983 ALLOCATE (num_species_small(2))
984 num_species_small(:) = (/1, num_species/)
985 END IF
986
987 DO k = 1, SIZE(num_species_small) - 1
988 IF (num_species > atoms_memory) THEN
989 ALLOCATE (cores(num_species_small(k + 1) - (num_species_small(k) - 1)))
990 ELSE
991 ALLOCATE (cores(num_species))
992 END IF
993
994 DO i = num_species_small(k), num_species_small(k + 1)
995 CALL rs_grid_create(rs_single_dr(i), auxbas_rs_desc)
996 CALL rs_grid_zero(rs_single_dr(i))
997 END DO
998 DO iex = 1, numexp
999
1000 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
1001 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
1002 prefactor = 2.0_dp*alpha
1003 npme = 0
1004 cores = 0
1005
1006 DO iatom = 1, SIZE(cores)
1007 atom_a = atom_list(iatom + (num_species_small(k) - 1))
1008 ra(:) = pbc(particle_set(atom_a)%r, cell)
1009
1010 IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
1011 IF (modulo(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
1012 npme = npme + 1
1013 cores(npme) = iatom
1014 END IF
1015 ELSE
1016 npme = npme + 1
1017 cores(npme) = iatom
1018 END IF
1019 END DO
1020 DO j = 1, npme
1021 iatom = cores(j)
1022 atom_a = atom_list(iatom + (num_species_small(k) - 1))
1023 pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
1024 ra(:) = pbc(particle_set(atom_a)%r, cell)
1025 subpatch_pattern = 0
1026
1027 ! Calculate cutoff
1028 IF (hirshfeld_control%use_atomic_cutoff) THEN
1029 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
1030 ra=ra, rb=ra, rp=ra, &
1031 zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
1032 pab=pab, o1=0, o2=0, & ! without map_consistent
1033 prefactor=1.0_dp, cutoff=0.0_dp)
1034 END IF
1035
1036 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1037 (/0.0_dp, 0.0_dp, 0.0_dp/), prefactor, &
1038 pab, 0, 0, rs_single_dr(iatom + (num_species_small(k) - 1)), &
1039 radius=radius, &
1040 ga_gb_function=grid_func_ab, use_subpatch=.true., &
1041 subpatch_pattern=subpatch_pattern)
1042
1043 END DO
1044 END DO
1045
1046 DO iatom = num_species_small(k), num_species_small(k + 1)
1047 CALL transfer_rs2pw(rs_single_dr(iatom), pw_single_dr(iatom))
1048 CALL rs_grid_release(rs_single_dr(iatom))
1049 END DO
1050
1051 DEALLOCATE (cores)
1052 END DO
1053
1054 DO iatom = 1, num_species
1055 atom_a = atom_list(iatom)
1056 cdft_control%group(igroup)%gradients_x(atom_a, :, :, :) = pw_single_dr(iatom)%array(:, :, :)
1057 CALL auxbas_pw_pool%give_back_pw(pw_single_dr(iatom))
1058 END DO
1059
1060 DEALLOCATE (rs_single_dr)
1061 DEALLOCATE (num_species_small)
1062 DEALLOCATE (pw_single_dr)
1063 END DO
1064 DEALLOCATE (pab)
1065 END DO
1066 END IF
1067
1068 IF (cdft_control%in_memory) THEN
1069 DO igroup = 1, SIZE(cdft_control%group)
1070 ALLOCATE (cdft_control%group(igroup)%gradients_y(1*num_atoms, lb_pw(1):ub_pw(1), &
1071 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1072 ALLOCATE (cdft_control%group(igroup)%gradients_z(1*num_atoms, lb_pw(1):ub_pw(1), &
1073 lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1074 cdft_control%group(igroup)%gradients_y(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1075 cdft_control%group(igroup)%gradients_z(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1076 END DO
1077 END IF
1078
1079 ! Calculate gradient if requested
1080 IF (cdft_control%in_memory) THEN
1081
1082 DO igroup = 1, SIZE(cdft_control%group)
1083 DO k = lb_pw(3), ub_pw(3)
1084 DO j = lb_pw(2), ub_pw(2)
1085 DO i = lb_pw(1), ub_pw(1)
1086 DO iatom = 1, natom
1087
1088 ra(:) = particle_set(iatom)%r
1089
1090 IF (cdft_control%group(igroup)%hw_rho_total%array(i, j, k) > hirshfeld_control%eps_cutoff) THEN
1091
1092 exp_eval = (coefficients(iatom) - &
1093 cdft_control%group(igroup)%weight%array(i, j, k))/ &
1094 cdft_control%group(igroup)%hw_rho_total%array(i, j, k)
1095
1096 r2 = (/i*dr_pw(1), j*dr_pw(2), k*dr_pw(3)/) + origin
1097 r_pbc = pbc(ra, r2, cell)
1098
1099 ! Store gradient d/dR_x w, including term: (r_x - R_x)
1100 cdft_control%group(igroup)%gradients_x(iatom, i, j, k) = &
1101 cdft_control%group(igroup)%gradients_x(iatom, i, j, k)* &
1102 r_pbc(1)*exp_eval
1103
1104 ! Store gradient d/dR_y w, including term: (r_y - R_y)
1105 cdft_control%group(igroup)%gradients_y(iatom, i, j, k) = &
1106 cdft_control%group(igroup)%gradients_y(iatom, i, j, k)* &
1107 r_pbc(2)*exp_eval
1108
1109 ! Store gradient d/dR_z w, including term:(r_z - R_z)
1110 cdft_control%group(igroup)%gradients_z(iatom, i, j, k) = &
1111 cdft_control%group(igroup)%gradients_z(iatom, i, j, k)* &
1112 r_pbc(3)*exp_eval
1113
1114 END IF
1115 END DO
1116 END DO
1117 END DO
1118 END DO
1119
1120 IF (igroup == 1) THEN
1121 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(1)%hw_rho_total)
1122 END IF
1123 END DO
1124 END IF
1125
1126 CALL rs_grid_release(rs_rho_all)
1127
1128 IF (ALLOCATED(coefficients)) DEALLOCATE (coefficients)
1129 IF (ALLOCATED(is_constraint)) DEALLOCATE (is_constraint)
1130
1131 CALL timestop(handle)
1132
1133 END SUBROUTINE hirshfeld_constraint_low
1134
1135! **************************************************************************************************
1136!> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT
1137!> weight function and the realspace electron density
1138!> \param qs_env ...
1139! **************************************************************************************************
1140 SUBROUTINE cdft_constraint_integrate(qs_env)
1141 TYPE(qs_environment_type), POINTER :: qs_env
1142
1143 CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_integrate'
1144
1145 INTEGER :: handle, i, iatom, igroup, ikind, ivar, &
1146 iw, jatom, natom, nvar
1147 LOGICAL :: is_becke, paw_atom
1148 REAL(kind=dp) :: dvol, eps_cavity, sign
1149 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: de, strength, target_val
1150 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge, gapw_offset
1151 TYPE(becke_constraint_type), POINTER :: becke_control
1152 TYPE(cdft_control_type), POINTER :: cdft_control
1153 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1154 TYPE(cp_logger_type), POINTER :: logger
1155 TYPE(dft_control_type), POINTER :: dft_control
1156 TYPE(mp_para_env_type), POINTER :: para_env
1157 TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
1158 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1159 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge, rho_r
1160 TYPE(qs_energy_type), POINTER :: energy
1161 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1162 TYPE(qs_rho_type), POINTER :: rho
1163 TYPE(rho0_mpole_type), POINTER :: rho0_mpole
1164 TYPE(section_vals_type), POINTER :: cdft_constraint_section
1165
1166 NULLIFY (para_env, dft_control, particle_set, rho_r, energy, rho, &
1167 logger, cdft_constraint_section, qs_kind_set, mp_rho, &
1168 rho0_mpole, group, charge)
1169 CALL timeset(routinen, handle)
1170 logger => cp_get_default_logger()
1171 CALL get_qs_env(qs_env, &
1172 particle_set=particle_set, &
1173 rho=rho, &
1174 natom=natom, &
1175 dft_control=dft_control, &
1176 para_env=para_env, &
1177 qs_kind_set=qs_kind_set)
1178 CALL qs_rho_get(rho, rho_r=rho_r)
1179 cpassert(ASSOCIATED(qs_kind_set))
1180 cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1181 iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
1182 cdft_control => dft_control%qs_control%cdft_control
1183 is_becke = (cdft_control%type == outer_scf_becke_constraint)
1184 becke_control => cdft_control%becke_control
1185 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
1186 cpabort("Becke control has not been allocated.")
1187 group => cdft_control%group
1188 ! Initialize
1189 nvar = SIZE(cdft_control%target)
1190 ALLOCATE (strength(nvar))
1191 ALLOCATE (target_val(nvar))
1192 ALLOCATE (de(nvar))
1193 strength(:) = cdft_control%strength(:)
1194 target_val(:) = cdft_control%target(:)
1195 sign = 1.0_dp
1196 de = 0.0_dp
1197 dvol = group(1)%weight%pw_grid%dvol
1198 IF (cdft_control%atomic_charges) THEN
1199 charge => cdft_control%charge
1200 ALLOCATE (electronic_charge(cdft_control%natoms, dft_control%nspins))
1201 electronic_charge = 0.0_dp
1202 END IF
1203 ! Calculate value of constraint i.e. int ( rho(r) w(r) dr)
1204 DO i = 1, dft_control%nspins
1205 DO igroup = 1, SIZE(group)
1206 SELECT CASE (group(igroup)%constraint_type)
1207 CASE (cdft_charge_constraint)
1208 sign = 1.0_dp
1209 CASE (cdft_magnetization_constraint)
1210 IF (i == 1) THEN
1211 sign = 1.0_dp
1212 ELSE
1213 sign = -1.0_dp
1214 END IF
1215 CASE (cdft_alpha_constraint)
1216 sign = 1.0_dp
1217 IF (i == 2) cycle
1218 CASE (cdft_beta_constraint)
1219 sign = 1.0_dp
1220 IF (i == 1) cycle
1221 CASE DEFAULT
1222 cpabort("Unknown constraint type.")
1223 END SELECT
1224 IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1225 ! With external control, we can use cavity_mat as a mask to kahan sum
1226 eps_cavity = becke_control%eps_cavity
1227 IF (igroup /= 1) &
1228 CALL cp_abort(__location__, &
1229 "Multiple constraints not yet supported by parallel mixed calculations.")
1230 de(igroup) = de(igroup) + sign*accurate_dot_product(group(igroup)%weight%array, rho_r(i)%array, &
1231 becke_control%cavity_mat, eps_cavity)*dvol
1232 ELSE
1233 de(igroup) = de(igroup) + sign*pw_integral_ab(group(igroup)%weight, rho_r(i), local_only=.true.)
1234 END IF
1235 END DO
1236 IF (cdft_control%atomic_charges) THEN
1237 DO iatom = 1, cdft_control%natoms
1238 electronic_charge(iatom, i) = pw_integral_ab(charge(iatom), rho_r(i), local_only=.true.)
1239 END DO
1240 END IF
1241 END DO
1242 CALL get_qs_env(qs_env, energy=energy)
1243 CALL para_env%sum(de)
1244 IF (cdft_control%atomic_charges) THEN
1245 CALL para_env%sum(electronic_charge)
1246 END IF
1247 ! Use fragment densities as reference value (= Becke deformation density)
1248 IF (cdft_control%fragment_density .AND. .NOT. cdft_control%fragments_integrated) THEN
1249 CALL prepare_fragment_constraint(qs_env)
1250 END IF
1251 IF (dft_control%qs_control%gapw) THEN
1252 ! GAPW: add core charges (rho_hard - rho_soft)
1253 IF (cdft_control%fragment_density) &
1254 CALL cp_abort(__location__, &
1255 "Fragment constraints not yet compatible with GAPW.")
1256 ALLOCATE (gapw_offset(nvar, dft_control%nspins))
1257 gapw_offset = 0.0_dp
1258 CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole)
1259 CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho)
1260 DO i = 1, dft_control%nspins
1261 DO igroup = 1, SIZE(group)
1262 DO iatom = 1, SIZE(group(igroup)%atoms)
1263 SELECT CASE (group(igroup)%constraint_type)
1264 CASE (cdft_charge_constraint)
1265 sign = 1.0_dp
1266 CASE (cdft_magnetization_constraint)
1267 IF (i == 1) THEN
1268 sign = 1.0_dp
1269 ELSE
1270 sign = -1.0_dp
1271 END IF
1272 CASE (cdft_alpha_constraint)
1273 sign = 1.0_dp
1274 IF (i == 2) cycle
1275 CASE (cdft_beta_constraint)
1276 sign = 1.0_dp
1277 IF (i == 1) cycle
1278 CASE DEFAULT
1279 cpabort("Unknown constraint type.")
1280 END SELECT
1281 jatom = group(igroup)%atoms(iatom)
1282 CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
1283 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
1284 IF (paw_atom) THEN
1285 gapw_offset(igroup, i) = gapw_offset(igroup, i) + sign*group(igroup)%coeff(iatom)*mp_rho(jatom)%q0(i)
1286 END IF
1287 END DO
1288 END DO
1289 END DO
1290 IF (cdft_control%atomic_charges) THEN
1291 DO iatom = 1, cdft_control%natoms
1292 jatom = cdft_control%atoms(iatom)
1293 CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
1294 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
1295 IF (paw_atom) THEN
1296 DO i = 1, dft_control%nspins
1297 electronic_charge(iatom, i) = electronic_charge(iatom, i) + mp_rho(jatom)%q0(i)
1298 END DO
1299 END IF
1300 END DO
1301 END IF
1302 DO i = 1, dft_control%nspins
1303 DO ivar = 1, nvar
1304 de(ivar) = de(ivar) + gapw_offset(ivar, i)
1305 END DO
1306 END DO
1307 DEALLOCATE (gapw_offset)
1308 END IF
1309 ! Update constraint value and energy
1310 cdft_control%value(:) = de(:)
1311 energy%cdft = 0.0_dp
1312 DO ivar = 1, nvar
1313 energy%cdft = energy%cdft + (de(ivar) - target_val(ivar))*strength(ivar)
1314 END DO
1315 ! Print constraint info and atomic CDFT charges
1316 CALL cdft_constraint_print(qs_env, electronic_charge)
1317 ! Deallocate tmp storage
1318 DEALLOCATE (de, strength, target_val)
1319 IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge)
1320 CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO")
1321 CALL timestop(handle)
1322
1323 END SUBROUTINE cdft_constraint_integrate
1324
1325! **************************************************************************************************
1326!> \brief Calculates atomic forces due to a CDFT constraint (Becke or Hirshfeld)
1327!> \param qs_env ...
1328! **************************************************************************************************
1329 SUBROUTINE cdft_constraint_force(qs_env)
1330 TYPE(qs_environment_type), POINTER :: qs_env
1331
1332 CHARACTER(len=*), PARAMETER :: routinen = 'cdft_constraint_force'
1333
1334 INTEGER :: handle, i, iatom, igroup, ikind, ispin, &
1335 j, k, natom, nvar
1336 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
1337 INTEGER, DIMENSION(2, 3) :: bo
1338 INTEGER, DIMENSION(3) :: lb, ub
1339 REAL(kind=dp) :: dvol, eps_cavity, sign
1340 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength
1341 REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
1342 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1343 TYPE(becke_constraint_type), POINTER :: becke_control
1344 TYPE(cdft_control_type), POINTER :: cdft_control
1345 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1346 TYPE(cell_type), POINTER :: cell
1347 TYPE(dft_control_type), POINTER :: dft_control
1348 TYPE(mp_para_env_type), POINTER :: para_env
1349 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1350 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1351 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1352 TYPE(qs_rho_type), POINTER :: rho
1353
1354 CALL timeset(routinen, handle)
1355 NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, &
1356 rho, rho_r, force, cutoffs, becke_control, group)
1357
1358 CALL get_qs_env(qs_env, &
1359 atomic_kind_set=atomic_kind_set, &
1360 natom=natom, &
1361 particle_set=particle_set, &
1362 cell=cell, &
1363 rho=rho, &
1364 force=force, &
1365 dft_control=dft_control, &
1366 para_env=para_env)
1367 CALL qs_rho_get(rho, rho_r=rho_r)
1368
1369 cdft_control => dft_control%qs_control%cdft_control
1370 becke_control => cdft_control%becke_control
1371 group => cdft_control%group
1372 nvar = SIZE(cdft_control%target)
1373 ALLOCATE (strength(nvar))
1374 strength(:) = cdft_control%strength(:)
1375 cutoffs => cdft_control%becke_control%cutoffs
1376 eps_cavity = cdft_control%becke_control%eps_cavity
1377
1378 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1379 atom_of_kind=atom_of_kind, &
1380 kind_of=kind_of)
1381 DO igroup = 1, SIZE(cdft_control%group)
1382 ALLOCATE (cdft_control%group(igroup)%integrated(3, natom))
1383 cdft_control%group(igroup)%integrated = 0.0_dp
1384 END DO
1385
1386 lb(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
1387 ub(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
1388 bo = cdft_control%group(1)%weight%pw_grid%bounds_local
1389 dvol = cdft_control%group(1)%weight%pw_grid%dvol
1390 sign = 1.0_dp
1391
1392 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1393 IF (.NOT. cdft_control%becke_control%in_memory) THEN
1394 CALL becke_constraint_low(qs_env, just_gradients=.true.)
1395 END IF
1396
1397 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1398 IF (.NOT. cdft_control%in_memory) THEN
1399 CALL hirshfeld_constraint_low(qs_env, just_gradients=.true.)
1400 END IF
1401 END IF
1402
1403 ! If no Becke Gaussian confinement
1404 IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN
1405 ! No external control
1406 DO k = bo(1, 1), bo(2, 1)
1407 DO j = bo(1, 2), bo(2, 2)
1408 DO i = bo(1, 3), bo(2, 3)
1409 ! First check if this grid point should be skipped
1410 IF (cdft_control%becke_control%cavity_confine) THEN
1411 IF (cdft_control%becke_control%cavity%array(k, j, i) < eps_cavity) cycle
1412 END IF
1413
1414 DO igroup = 1, SIZE(cdft_control%group)
1415 DO iatom = 1, natom
1416 DO ispin = 1, dft_control%nspins
1417
1418 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1419 CASE (cdft_charge_constraint)
1420 sign = 1.0_dp
1421 CASE (cdft_magnetization_constraint)
1422 IF (ispin == 1) THEN
1423 sign = 1.0_dp
1424 ELSE
1425 sign = -1.0_dp
1426 END IF
1427 CASE (cdft_alpha_constraint)
1428 sign = 1.0_dp
1429 IF (ispin == 2) cycle
1430 CASE (cdft_beta_constraint)
1431 sign = 1.0_dp
1432 IF (ispin == 1) cycle
1433 CASE DEFAULT
1434 cpabort("Unknown constraint type.")
1435 END SELECT
1436
1437 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1438
1439 cdft_control%group(igroup)%integrated(:, iatom) = &
1440 cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1441 cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1442 *rho_r(ispin)%array(k, j, i) &
1443 *dvol
1444
1445 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1446
1447 cdft_control%group(igroup)%integrated(1, iatom) = &
1448 cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1449 cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1450 *rho_r(ispin)%array(k, j, i) &
1451 *dvol
1452
1453 cdft_control%group(igroup)%integrated(2, iatom) = &
1454 cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1455 cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1456 *rho_r(ispin)%array(k, j, i) &
1457 *dvol
1458
1459 cdft_control%group(igroup)%integrated(3, iatom) = &
1460 cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1461 cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1462 *rho_r(ispin)%array(k, j, i) &
1463 *dvol
1464
1465 END IF
1466
1467 END DO
1468 END DO
1469 END DO
1470 END DO
1471 END DO
1472 END DO
1473
1474 ! If Becke Gaussian confinement
1475 ELSE
1476 DO k = lbound(cdft_control%becke_control%cavity_mat, 1), ubound(cdft_control%becke_control%cavity_mat, 1)
1477 DO j = lbound(cdft_control%becke_control%cavity_mat, 2), ubound(cdft_control%becke_control%cavity_mat, 2)
1478 DO i = lbound(cdft_control%becke_control%cavity_mat, 3), ubound(cdft_control%becke_control%cavity_mat, 3)
1479
1480 ! First check if this grid point should be skipped
1481 IF (cdft_control%becke_control%cavity_mat(k, j, i) < eps_cavity) cycle
1482
1483 DO igroup = 1, SIZE(group)
1484 DO iatom = 1, natom
1485 DO ispin = 1, dft_control%nspins
1486 SELECT CASE (group(igroup)%constraint_type)
1487 CASE (cdft_charge_constraint)
1488 sign = 1.0_dp
1489 CASE (cdft_magnetization_constraint)
1490 IF (ispin == 1) THEN
1491 sign = 1.0_dp
1492 ELSE
1493 sign = -1.0_dp
1494 END IF
1495 CASE (cdft_alpha_constraint)
1496 sign = 1.0_dp
1497 IF (ispin == 2) cycle
1498 CASE (cdft_beta_constraint)
1499 sign = 1.0_dp
1500 IF (ispin == 1) cycle
1501 CASE DEFAULT
1502 cpabort("Unknown constraint type.")
1503 END SELECT
1504
1505 ! Integrate gradient of weight function
1506 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1507
1508 cdft_control%group(igroup)%integrated(:, iatom) = &
1509 cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1510 cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1511 *rho_r(ispin)%array(k, j, i) &
1512 *dvol
1513
1514 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1515
1516 cdft_control%group(igroup)%integrated(1, iatom) = &
1517 cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1518 cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1519 *rho_r(ispin)%array(k, j, i) &
1520 *dvol
1521
1522 cdft_control%group(igroup)%integrated(2, iatom) = &
1523 cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1524 cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1525 *rho_r(ispin)%array(k, j, i) &
1526 *dvol
1527
1528 cdft_control%group(igroup)%integrated(3, iatom) = &
1529 cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1530 cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1531 *rho_r(ispin)%array(k, j, i) &
1532 *dvol
1533
1534 END IF
1535
1536 END DO
1537 END DO
1538 END DO
1539 END DO
1540 END DO
1541 END DO
1542 END IF
1543
1544 IF (.NOT. cdft_control%transfer_pot) THEN
1545 IF (cdft_control%type == outer_scf_becke_constraint) THEN
1546 DO igroup = 1, SIZE(group)
1547 DEALLOCATE (cdft_control%group(igroup)%gradients)
1548 END DO
1549 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1550 DO igroup = 1, SIZE(group)
1551 DEALLOCATE (cdft_control%group(igroup)%gradients_x)
1552 DEALLOCATE (cdft_control%group(igroup)%gradients_y)
1553 DEALLOCATE (cdft_control%group(igroup)%gradients_z)
1554 END DO
1555 END IF
1556 END IF
1557
1558 DO igroup = 1, SIZE(group)
1559 CALL para_env%sum(group(igroup)%integrated)
1560 END DO
1561
1562 ! Update force only on master process. Otherwise force due to constraint becomes multiplied
1563 ! by the number of processes when the final force%rho_elec is constructed in qs_force
1564 ! by mp_summing [the final integrated(:,:) is distributed on all processors]
1565 IF (para_env%is_source()) THEN
1566 DO igroup = 1, SIZE(group)
1567 DO iatom = 1, natom
1568 ikind = kind_of(iatom)
1569 i = atom_of_kind(iatom)
1570 force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup)
1571 END DO
1572 END DO
1573 END IF
1574
1575 DEALLOCATE (strength)
1576 DO igroup = 1, SIZE(group)
1577 DEALLOCATE (group(igroup)%integrated)
1578 END DO
1579 NULLIFY (group)
1580
1581 CALL timestop(handle)
1582
1583 END SUBROUTINE cdft_constraint_force
1584
1585! **************************************************************************************************
1586!> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied
1587!> by the CDFT weight functions and integrated over the realspace grid.
1588!> \param qs_env ...
1589! **************************************************************************************************
1590 SUBROUTINE prepare_fragment_constraint(qs_env)
1591 TYPE(qs_environment_type), POINTER :: qs_env
1592
1593 CHARACTER(len=*), PARAMETER :: routinen = 'prepare_fragment_constraint'
1594
1595 INTEGER :: handle, i, iatom, igroup, natom, &
1596 nelectron_total, nfrag_spins
1597 LOGICAL :: is_becke, needs_spin_density
1598 REAL(kind=dp) :: dvol, multiplier(2), nelectron_frag
1599 TYPE(becke_constraint_type), POINTER :: becke_control
1600 TYPE(cdft_control_type), POINTER :: cdft_control
1601 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1602 TYPE(cp_logger_type), POINTER :: logger
1603 TYPE(dft_control_type), POINTER :: dft_control
1604 TYPE(mp_para_env_type), POINTER :: para_env
1605 TYPE(pw_env_type), POINTER :: pw_env
1606 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1607 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: rho_frag
1608 TYPE(qs_subsys_type), POINTER :: subsys
1609
1610 NULLIFY (para_env, dft_control, logger, subsys, pw_env, auxbas_pw_pool, group)
1611 CALL timeset(routinen, handle)
1612 logger => cp_get_default_logger()
1613 CALL get_qs_env(qs_env, &
1614 natom=natom, &
1615 dft_control=dft_control, &
1616 para_env=para_env)
1617
1618 cdft_control => dft_control%qs_control%cdft_control
1619 is_becke = (cdft_control%type == outer_scf_becke_constraint)
1620 becke_control => cdft_control%becke_control
1621 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
1622 cpabort("Becke control has not been allocated.")
1623 group => cdft_control%group
1624 dvol = group(1)%weight%pw_grid%dvol
1625 ! Fragment densities are meaningful only for some calculation types
1626 IF (.NOT. qs_env%single_point_run) &
1627 CALL cp_abort(__location__, &
1628 "CDFT fragment constraints are only compatible with single "// &
1629 "point calculations (run_type ENERGY or ENERGY_FORCE).")
1630 IF (dft_control%qs_control%gapw) &
1631 CALL cp_abort(__location__, &
1632 "CDFT fragment constraint not compatible with GAPW.")
1633 needs_spin_density = .false.
1634 multiplier = 1.0_dp
1635 nfrag_spins = 1
1636 DO igroup = 1, SIZE(group)
1637 SELECT CASE (group(igroup)%constraint_type)
1638 CASE (cdft_charge_constraint)
1639 ! Do nothing
1640 CASE (cdft_magnetization_constraint)
1641 needs_spin_density = .true.
1642 CASE (cdft_alpha_constraint, cdft_beta_constraint)
1643 CALL cp_abort(__location__, &
1644 "CDFT fragment constraint not yet compatible with "// &
1645 "spin specific constraints.")
1646 CASE DEFAULT
1647 cpabort("Unknown constraint type.")
1648 END SELECT
1649 END DO
1650 IF (needs_spin_density) THEN
1651 nfrag_spins = 2
1652 DO i = 1, 2
1653 IF (cdft_control%flip_fragment(i)) multiplier(i) = -1.0_dp
1654 END DO
1655 END IF
1656 ! Read fragment reference densities
1657 ALLOCATE (cdft_control%fragments(nfrag_spins, 2))
1658 ALLOCATE (rho_frag(nfrag_spins))
1659 CALL get_qs_env(qs_env, pw_env=pw_env)
1660 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1661 ! Total density (rho_alpha + rho_beta)
1662 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 1))
1663 CALL cp_cube_to_pw(cdft_control%fragments(1, 1), &
1664 cdft_control%fragment_a_fname, 1.0_dp)
1665 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 2))
1666 CALL cp_cube_to_pw(cdft_control%fragments(1, 2), &
1667 cdft_control%fragment_b_fname, 1.0_dp)
1668 ! Spin difference density (rho_alpha - rho_beta) if needed
1669 IF (needs_spin_density) THEN
1670 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 1))
1671 CALL cp_cube_to_pw(cdft_control%fragments(2, 1), &
1672 cdft_control%fragment_a_spin_fname, multiplier(1))
1673 CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 2))
1674 CALL cp_cube_to_pw(cdft_control%fragments(2, 2), &
1675 cdft_control%fragment_b_spin_fname, multiplier(2))
1676 END IF
1677 ! Sum up fragments
1678 DO i = 1, nfrag_spins
1679 CALL auxbas_pw_pool%create_pw(rho_frag(i))
1680 CALL pw_copy(cdft_control%fragments(i, 1), rho_frag(i))
1681 CALL pw_axpy(cdft_control%fragments(i, 2), rho_frag(i), 1.0_dp)
1682 CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 1))
1683 CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 2))
1684 END DO
1685 DEALLOCATE (cdft_control%fragments)
1686 ! Check that the number of electrons is consistent
1687 CALL get_qs_env(qs_env, subsys=subsys)
1688 CALL qs_subsys_get(subsys, nelectron_total=nelectron_total)
1689 nelectron_frag = pw_integrate_function(rho_frag(1))
1690 IF (nint(nelectron_frag) /= nelectron_total) &
1691 CALL cp_abort(__location__, &
1692 "The number of electrons in the reference and interacting "// &
1693 "configurations does not match. Check your fragment cube files.")
1694 ! Update constraint target value i.e. perform integration w_i*rho_frag_{tot/spin}*dr
1695 cdft_control%target = 0.0_dp
1696 DO igroup = 1, SIZE(group)
1697 IF (group(igroup)%constraint_type == cdft_charge_constraint) THEN
1698 i = 1
1699 ELSE
1700 i = 2
1701 END IF
1702 IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1703 cdft_control%target(igroup) = cdft_control%target(igroup) + &
1704 accurate_dot_product(group(igroup)%weight%array, rho_frag(i)%array, &
1705 becke_control%cavity_mat, becke_control%eps_cavity)*dvol
1706 ELSE
1707 cdft_control%target(igroup) = cdft_control%target(igroup) + &
1708 pw_integral_ab(group(igroup)%weight, rho_frag(i), local_only=.true.)
1709 END IF
1710 END DO
1711 CALL para_env%sum(cdft_control%target)
1712 ! Calculate reference atomic charges int( w_i * rho_frag * dr )
1713 IF (cdft_control%atomic_charges) THEN
1714 ALLOCATE (cdft_control%charges_fragment(cdft_control%natoms, nfrag_spins))
1715 DO i = 1, nfrag_spins
1716 DO iatom = 1, cdft_control%natoms
1717 cdft_control%charges_fragment(iatom, i) = &
1718 pw_integral_ab(cdft_control%charge(iatom), rho_frag(i), local_only=.true.)
1719 END DO
1720 END DO
1721 CALL para_env%sum(cdft_control%charges_fragment)
1722 END IF
1723 DO i = 1, nfrag_spins
1724 CALL auxbas_pw_pool%give_back_pw(rho_frag(i))
1725 END DO
1726 DEALLOCATE (rho_frag)
1727 cdft_control%fragments_integrated = .true.
1728
1729 CALL timestop(handle)
1730
1731 END SUBROUTINE prepare_fragment_constraint
1732
1733END 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:117
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:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_cube_to_pw
subroutine, public cp_cube_to_pw(grid, filename, scaling, silent)
Thin wrapper around routine cube_to_pw.
Definition cp_realspace_grid_cube.F:113
grid_api
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
grid_api::grid_func_ab
integer, parameter, public grid_func_ab
Definition grid_api.F:29
grid_api::collocate_pgf_product
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
Definition grid_api.F:119
hirshfeld_types
The types needed for the calculation of Hirshfeld charges and related functions.
Definition hirshfeld_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition input_constants.F:766
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition input_constants.F:766
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition input_constants.F:766
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition input_constants.F:718
input_constants::outer_scf_hirshfeld_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Definition input_constants.F:718
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition input_constants.F:766
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_cdft_methods
Subroutines for building CDFT constraints.
Definition qs_cdft_methods.F:14
qs_cdft_methods::becke_constraint
subroutine, public becke_constraint(qs_env, calc_pot, calculate_forces)
Driver routine for calculating a Becke constraint.
Definition qs_cdft_methods.F:105
qs_cdft_methods::hirshfeld_constraint
subroutine, public hirshfeld_constraint(qs_env, calc_pot, calculate_forces)
Driver routine for calculating a Hirshfeld constraint.
Definition qs_cdft_methods.F:588
qs_cdft_types
Defines CDFT control structures.
Definition qs_cdft_types.F:14
qs_cdft_utils
Utility subroutines for CDFT calculations.
Definition qs_cdft_utils.F:14
qs_cdft_utils::cdft_constraint_print
subroutine, public cdft_constraint_print(qs_env, electronic_charge)
Prints information about CDFT constraints.
Definition qs_cdft_utils.F:1103
qs_cdft_utils::hirshfeld_constraint_init
subroutine, public hirshfeld_constraint_init(qs_env)
Initializes Gaussian Hirshfeld constraints.
Definition qs_cdft_utils.F:974
qs_cdft_utils::becke_constraint_init
subroutine, public becke_constraint_init(qs_env)
Initializes the Becke constraint environment.
Definition qs_cdft_utils.F:97
qs_cdft_utils::hfun_scale
subroutine, public hfun_scale(fout, fun1, fun2, divide, small)
Calculate fout = fun1/fun2 or fout = fun1*fun2.
Definition qs_cdft_utils.F:856
qs_cdft_utils::cdft_print_hirshfeld_density
subroutine, public cdft_print_hirshfeld_density(qs_env)
Prints Hirshfeld weight function and promolecule density.
Definition qs_cdft_utils.F:1313
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
qs_rho0_types
Definition qs_rho0_types.F:8
qs_rho0_types::get_rho0_mpole
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, qlm_gg, qlm_car, qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs)
...
Definition qs_rho0_types.F:425
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:134
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition realspace_grid_types.F:482
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition realspace_grid_types.F:656
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition realspace_grid_types.F:2016
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1909
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:559
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:126
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:220
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:215
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:137