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