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