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