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