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 20 : 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 616 : 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 144 : 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 : hirshfeld_control%print_density = .FALSE.
687 : END IF
688 :
689 72 : ALLOCATE (coefficients(natom))
690 72 : ALLOCATE (is_constraint(natom))
691 :
692 24 : subpatch_pattern = 0
693 24 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
694 24 : radius = 100.0_dp
695 :
696 24 : dr_pw(1) = rho_r(1)%pw_grid%dr(1)
697 24 : dr_pw(2) = rho_r(1)%pw_grid%dr(2)
698 24 : dr_pw(3) = rho_r(1)%pw_grid%dr(3)
699 96 : lb_pw(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
700 96 : ub_pw(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
701 96 : npts = rho_r(1)%pw_grid%npts
702 24 : origin(1) = (dr_pw(1)*npts(1))*0.5_dp
703 24 : origin(2) = (dr_pw(2)*npts(2))*0.5_dp
704 24 : origin(3) = (dr_pw(3)*npts(3))*0.5_dp
705 :
706 : CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
707 24 : auxbas_pw_pool=auxbas_pw_pool)
708 24 : CALL rs_grid_create(rs_rho_all, auxbas_rs_desc)
709 24 : CALL rs_grid_zero(rs_rho_all)
710 :
711 24 : dr_rs(1) = rs_rho_all%desc%dh(1, 1)
712 24 : dr_rs(2) = rs_rho_all%desc%dh(2, 2)
713 24 : dr_rs(3) = rs_rho_all%desc%dh(3, 3)
714 24 : lb_rs(1) = LBOUND(rs_rho_all%r(:, :, :), 1)
715 24 : lb_rs(2) = LBOUND(rs_rho_all%r(:, :, :), 2)
716 24 : lb_rs(3) = LBOUND(rs_rho_all%r(:, :, :), 3)
717 24 : ub_rs(1) = UBOUND(rs_rho_all%r(:, :, :), 1)
718 24 : ub_rs(2) = UBOUND(rs_rho_all%r(:, :, :), 2)
719 24 : ub_rs(3) = UBOUND(rs_rho_all%r(:, :, :), 3)
720 :
721 : ! For each CDFT group
722 48 : DO igroup = 1, SIZE(cdft_control%group)
723 :
724 24 : IF (igroup == 2 .AND. .NOT. cdft_control%in_memory) THEN
725 0 : CALL rs_grid_zero(rs_rho_all)
726 : END IF
727 240 : bo = cdft_control%group(igroup)%weight%pw_grid%bounds_local
728 :
729 : ! Coefficients
730 72 : coefficients(:) = 0.0_dp
731 72 : is_constraint = .FALSE.
732 70 : DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
733 46 : coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
734 70 : is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE.
735 : END DO
736 :
737 : ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group
738 24 : CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc)
739 24 : CALL rs_grid_zero(rs_rho_constr)
740 :
741 : ! rs_single: Gaussian density over single atoms when required
742 24 : IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
743 0 : ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic(cdft_control%natoms))
744 0 : ALLOCATE (rs_single(cdft_control%natoms))
745 0 : DO i = 1, cdft_control%natoms
746 0 : CALL rs_grid_create(rs_single(i), auxbas_rs_desc)
747 0 : CALL rs_grid_zero(rs_single(i))
748 : END DO
749 : END IF
750 :
751 : ! Setup pw
752 24 : CALL pw_zero(cdft_control%group(igroup)%weight)
753 :
754 24 : CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
755 24 : CALL pw_set(cdft_control%group(igroup)%hw_rho_total_constraint, 1.0_dp)
756 :
757 24 : IF (igroup == 1) THEN
758 24 : CALL auxbas_pw_pool%create_pw(cdft_control%hw_rho_total)
759 24 : CALL pw_set(cdft_control%hw_rho_total, 1.0_dp)
760 :
761 24 : IF (hirshfeld_control%print_density) THEN
762 0 : DO iatom = 1, cdft_control%natoms
763 0 : CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic(iatom))
764 0 : CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic(iatom), 1.0_dp)
765 : END DO
766 : END IF
767 : END IF
768 :
769 24 : IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
770 40 : ALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge(cdft_control%natoms))
771 192 : ALLOCATE (rs_single_charge(cdft_control%natoms))
772 24 : ALLOCATE (compute_charge(natom))
773 24 : compute_charge = .FALSE.
774 :
775 24 : DO i = 1, cdft_control%natoms
776 16 : CALL rs_grid_create(rs_single_charge(i), auxbas_rs_desc)
777 16 : CALL rs_grid_zero(rs_single_charge(i))
778 24 : compute_charge(cdft_control%atoms(i)) = .TRUE.
779 : END DO
780 :
781 24 : DO iatom = 1, cdft_control%natoms
782 16 : CALL auxbas_pw_pool%create_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom))
783 24 : CALL pw_set(cdft_control%group(igroup)%hw_rho_atomic_charge(iatom), 1.0_dp)
784 : END DO
785 : END IF
786 :
787 24 : ALLOCATE (pab(1, 1))
788 24 : nthread = 1
789 24 : ithread = 0
790 :
791 72 : DO ikind = 1, SIZE(atomic_kind_set)
792 48 : numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
793 48 : IF (numexp <= 0) CYCLE
794 48 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
795 144 : ALLOCATE (cores(num_species))
796 :
797 180 : DO iex = 1, numexp
798 132 : alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
799 132 : coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
800 132 : npme = 0
801 264 : cores = 0
802 264 : DO iatom = 1, num_species
803 132 : atom_a = atom_list(iatom)
804 132 : ra(:) = pbc(particle_set(atom_a)%r, cell)
805 264 : IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
806 128 : IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
807 64 : npme = npme + 1
808 64 : cores(npme) = iatom
809 : END IF
810 : ELSE
811 4 : npme = npme + 1
812 4 : cores(npme) = iatom
813 : END IF
814 : END DO
815 248 : DO j = 1, npme
816 68 : iatom = cores(j)
817 68 : atom_a = atom_list(iatom)
818 68 : pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
819 68 : ra(:) = pbc(particle_set(atom_a)%r, cell)
820 :
821 68 : IF (hirshfeld_control%use_atomic_cutoff) THEN
822 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
823 : ra=ra, rb=ra, rp=ra, &
824 : zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
825 : pab=pab, o1=0, o2=0, & ! without map_consistent
826 68 : prefactor=1.0_dp, cutoff=0.0_dp)
827 : END IF
828 :
829 68 : IF (igroup == 1) THEN
830 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
831 : (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
832 : rs_rho_all, radius=radius, &
833 : ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
834 68 : subpatch_pattern=subpatch_pattern)
835 : END IF
836 :
837 68 : IF (is_constraint(atom_a)) THEN
838 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
839 : (/0.0_dp, 0.0_dp, 0.0_dp/), coefficients(atom_a), &
840 : pab, 0, 0, rs_rho_constr, &
841 : radius=radius, &
842 : ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
843 67 : subpatch_pattern=subpatch_pattern)
844 : END IF
845 :
846 68 : IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
847 0 : IF (is_constraint(atom_a)) THEN
848 0 : DO iatom = 1, cdft_control%natoms
849 0 : IF (atom_a == cdft_control%atoms(iatom)) EXIT
850 : END DO
851 0 : CPASSERT(iatom <= cdft_control%natoms)
852 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
853 : (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
854 : rs_single(iatom), radius=radius, &
855 : ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
856 0 : subpatch_pattern=subpatch_pattern)
857 : END IF
858 : END IF
859 :
860 200 : IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
861 22 : IF (compute_charge(atom_a)) THEN
862 33 : DO iatom = 1, cdft_control%natoms
863 33 : IF (atom_a == cdft_control%atoms(iatom)) EXIT
864 : END DO
865 22 : CPASSERT(iatom <= cdft_control%natoms)
866 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
867 : (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
868 : rs_single_charge(iatom), radius=radius, &
869 : ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
870 22 : subpatch_pattern=subpatch_pattern)
871 : END IF
872 : END IF
873 :
874 : END DO
875 : END DO
876 120 : DEALLOCATE (cores)
877 : END DO
878 24 : DEALLOCATE (pab)
879 :
880 24 : IF (igroup == 1) THEN
881 24 : CALL transfer_rs2pw(rs_rho_all, cdft_control%hw_rho_total)
882 : END IF
883 :
884 24 : CALL transfer_rs2pw(rs_rho_constr, cdft_control%group(igroup)%hw_rho_total_constraint)
885 24 : CALL rs_grid_release(rs_rho_constr)
886 :
887 : ! Calculate weight function
888 : CALL hfun_scale(cdft_control%group(igroup)%weight%array, &
889 : cdft_control%group(igroup)%hw_rho_total_constraint%array, &
890 : cdft_control%hw_rho_total%array, divide=.TRUE., &
891 24 : small=hirshfeld_control%eps_cutoff)
892 :
893 : ! Calculate charges
894 24 : IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
895 24 : DO i = 1, cdft_control%natoms
896 16 : CALL transfer_rs2pw(rs_single_charge(i), cdft_control%group(igroup)%hw_rho_atomic_charge(i))
897 : CALL hfun_scale(cdft_control%charge(i)%array, &
898 : cdft_control%group(igroup)%hw_rho_atomic_charge(i)%array, &
899 : cdft_control%hw_rho_total%array, divide=.TRUE., &
900 24 : small=hirshfeld_control%eps_cutoff)
901 : END DO
902 : END IF
903 :
904 : ! Print atomic densities if requested
905 48 : IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
906 0 : DO i = 1, cdft_control%natoms
907 0 : CALL transfer_rs2pw(rs_single(i), cdft_control%group(igroup)%hw_rho_atomic(i))
908 : END DO
909 0 : CALL cdft_print_hirshfeld_density(qs_env)
910 : END IF
911 :
912 : END DO
913 :
914 48 : DO igroup = 1, SIZE(cdft_control%group)
915 :
916 24 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_total_constraint)
917 :
918 24 : IF (.NOT. cdft_control%in_memory .AND. igroup == 1) THEN
919 22 : CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
920 : END IF
921 :
922 24 : IF (hirshfeld_control%print_density .AND. igroup == 1) THEN
923 0 : DO i = 1, cdft_control%natoms
924 0 : CALL rs_grid_release(rs_single(i))
925 0 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic(i))
926 : END DO
927 0 : DEALLOCATE (rs_single)
928 0 : DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic)
929 : END IF
930 :
931 48 : IF (cdft_control%atomic_charges .AND. igroup == 1) THEN
932 24 : DO i = 1, cdft_control%natoms
933 16 : CALL rs_grid_release(rs_single_charge(i))
934 24 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(igroup)%hw_rho_atomic_charge(i))
935 : END DO
936 24 : DEALLOCATE (rs_single_charge)
937 8 : DEALLOCATE (compute_charge)
938 8 : DEALLOCATE (cdft_control%group(igroup)%hw_rho_atomic_charge)
939 : END IF
940 :
941 : END DO
942 :
943 24 : IF (cdft_control%in_memory) THEN
944 4 : DO igroup = 1, SIZE(cdft_control%group)
945 : ALLOCATE (cdft_control%group(igroup)%gradients_x(1*natom, lb_pw(1):ub_pw(1), &
946 12 : lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
947 195284 : cdft_control%group(igroup)%gradients_x(:, :, :, :) = 0.0_dp
948 : END DO
949 : END IF
950 :
951 24 : IF (cdft_control%in_memory) THEN
952 4 : DO igroup = 1, SIZE(cdft_control%group)
953 :
954 2 : ALLOCATE (pab(1, 1))
955 2 : nthread = 1
956 2 : ithread = 0
957 2 : atoms_memory = hirshfeld_control%atoms_memory
958 :
959 6 : DO ikind = 1, SIZE(atomic_kind_set)
960 4 : numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
961 4 : IF (numexp <= 0) CYCLE
962 4 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=num_species, atom_list=atom_list)
963 :
964 16 : ALLOCATE (pw_single_dr(num_species))
965 96 : ALLOCATE (rs_single_dr(num_species))
966 :
967 8 : DO i = 1, num_species
968 4 : CALL auxbas_pw_pool%create_pw(pw_single_dr(i))
969 8 : CALL pw_zero(pw_single_dr(i))
970 : END DO
971 :
972 16 : atoms_memory_num = SIZE((/(j, j=1, num_species, atoms_memory)/))
973 :
974 : ! Can't store all pw grids, therefore split into groups of size atom_memory
975 : ! Ideally this code should be re-written to be more memory efficient
976 4 : IF (num_species > atoms_memory) THEN
977 0 : ALLOCATE (num_species_small(atoms_memory_num + 1))
978 0 : num_species_small(1:atoms_memory_num) = (/(j, j=1, num_species, atoms_memory)/)
979 0 : num_species_small(atoms_memory_num + 1) = num_species
980 : ELSE
981 4 : ALLOCATE (num_species_small(2))
982 12 : num_species_small(:) = (/1, num_species/)
983 : END IF
984 :
985 8 : DO k = 1, SIZE(num_species_small) - 1
986 4 : IF (num_species > atoms_memory) THEN
987 0 : ALLOCATE (cores(num_species_small(k + 1) - (num_species_small(k) - 1)))
988 : ELSE
989 12 : ALLOCATE (cores(num_species))
990 : END IF
991 :
992 8 : DO i = num_species_small(k), num_species_small(k + 1)
993 4 : CALL rs_grid_create(rs_single_dr(i), auxbas_rs_desc)
994 8 : CALL rs_grid_zero(rs_single_dr(i))
995 : END DO
996 36 : DO iex = 1, numexp
997 :
998 32 : alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
999 32 : coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
1000 32 : prefactor = 2.0_dp*alpha
1001 32 : npme = 0
1002 64 : cores = 0
1003 :
1004 64 : DO iatom = 1, SIZE(cores)
1005 32 : atom_a = atom_list(iatom + (num_species_small(k) - 1))
1006 32 : ra(:) = pbc(particle_set(atom_a)%r, cell)
1007 :
1008 64 : IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN
1009 32 : IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN
1010 16 : npme = npme + 1
1011 16 : cores(npme) = iatom
1012 : END IF
1013 : ELSE
1014 0 : npme = npme + 1
1015 0 : cores(npme) = iatom
1016 : END IF
1017 : END DO
1018 52 : DO j = 1, npme
1019 16 : iatom = cores(j)
1020 16 : atom_a = atom_list(iatom + (num_species_small(k) - 1))
1021 16 : pab(1, 1) = coef*hirshfeld_env%charges(atom_a)
1022 16 : ra(:) = pbc(particle_set(atom_a)%r, cell)
1023 : subpatch_pattern = 0
1024 :
1025 : ! Calculate cutoff
1026 16 : IF (hirshfeld_control%use_atomic_cutoff) THEN
1027 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
1028 : ra=ra, rb=ra, rp=ra, &
1029 : zetp=alpha, eps=hirshfeld_control%atomic_cutoff, &
1030 : pab=pab, o1=0, o2=0, & ! without map_consistent
1031 16 : prefactor=1.0_dp, cutoff=0.0_dp)
1032 : END IF
1033 :
1034 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
1035 : (/0.0_dp, 0.0_dp, 0.0_dp/), prefactor, &
1036 : pab, 0, 0, rs_single_dr(iatom + (num_species_small(k) - 1)), &
1037 : radius=radius, &
1038 : ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., &
1039 48 : subpatch_pattern=subpatch_pattern)
1040 :
1041 : END DO
1042 : END DO
1043 :
1044 8 : DO iatom = num_species_small(k), num_species_small(k + 1)
1045 4 : CALL transfer_rs2pw(rs_single_dr(iatom), pw_single_dr(iatom))
1046 8 : CALL rs_grid_release(rs_single_dr(iatom))
1047 : END DO
1048 :
1049 8 : DEALLOCATE (cores)
1050 : END DO
1051 :
1052 8 : DO iatom = 1, num_species
1053 4 : atom_a = atom_list(iatom)
1054 134564 : cdft_control%group(igroup)%gradients_x(atom_a, :, :, :) = pw_single_dr(iatom)%array(:, :, :)
1055 8 : CALL auxbas_pw_pool%give_back_pw(pw_single_dr(iatom))
1056 : END DO
1057 :
1058 8 : DEALLOCATE (rs_single_dr)
1059 4 : DEALLOCATE (num_species_small)
1060 10 : DEALLOCATE (pw_single_dr)
1061 : END DO
1062 4 : DEALLOCATE (pab)
1063 : END DO
1064 : END IF
1065 :
1066 24 : IF (cdft_control%in_memory) THEN
1067 4 : DO igroup = 1, SIZE(cdft_control%group)
1068 : ALLOCATE (cdft_control%group(igroup)%gradients_y(1*num_atoms, lb_pw(1):ub_pw(1), &
1069 12 : lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1070 : ALLOCATE (cdft_control%group(igroup)%gradients_z(1*num_atoms, lb_pw(1):ub_pw(1), &
1071 10 : lb_pw(2):ub_pw(2), lb_pw(3):ub_pw(3)))
1072 195282 : cdft_control%group(igroup)%gradients_y(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1073 195284 : cdft_control%group(igroup)%gradients_z(:, :, :, :) = cdft_control%group(igroup)%gradients_x(:, :, :, :)
1074 : END DO
1075 : END IF
1076 :
1077 : ! Calculate gradient if requested
1078 24 : IF (cdft_control%in_memory) THEN
1079 :
1080 4 : DO igroup = 1, SIZE(cdft_control%group)
1081 :
1082 : ! Coefficients
1083 6 : coefficients(:) = 0.0_dp
1084 6 : is_constraint = .FALSE.
1085 6 : DO i = 1, SIZE(cdft_control%group(igroup)%atoms)
1086 4 : coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i)
1087 6 : is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE.
1088 : END DO
1089 :
1090 84 : DO k = lb_pw(3), ub_pw(3)
1091 3282 : DO j = lb_pw(2), ub_pw(2)
1092 67280 : DO i = lb_pw(1), ub_pw(1)
1093 195200 : DO iatom = 1, natom
1094 :
1095 512000 : ra(:) = particle_set(iatom)%r
1096 :
1097 192000 : IF (cdft_control%hw_rho_total%array(i, j, k) > hirshfeld_control%eps_cutoff) THEN
1098 :
1099 : exp_eval = (coefficients(iatom) - &
1100 : cdft_control%group(igroup)%weight%array(i, j, k))/ &
1101 128000 : cdft_control%hw_rho_total%array(i, j, k)
1102 :
1103 512000 : r2 = (/i*dr_pw(1), j*dr_pw(2), k*dr_pw(3)/) + origin
1104 128000 : r_pbc = pbc(ra, r2, cell)
1105 :
1106 : ! Store gradient d/dR_x w, including term: (r_x - R_x)
1107 : cdft_control%group(igroup)%gradients_x(iatom, i, j, k) = &
1108 : cdft_control%group(igroup)%gradients_x(iatom, i, j, k)* &
1109 128000 : r_pbc(1)*exp_eval
1110 :
1111 : ! Store gradient d/dR_y w, including term: (r_y - R_y)
1112 : cdft_control%group(igroup)%gradients_y(iatom, i, j, k) = &
1113 : cdft_control%group(igroup)%gradients_y(iatom, i, j, k)* &
1114 128000 : r_pbc(2)*exp_eval
1115 :
1116 : ! Store gradient d/dR_z w, including term:(r_z - R_z)
1117 : cdft_control%group(igroup)%gradients_z(iatom, i, j, k) = &
1118 : cdft_control%group(igroup)%gradients_z(iatom, i, j, k)* &
1119 128000 : r_pbc(3)*exp_eval
1120 :
1121 : END IF
1122 : END DO
1123 : END DO
1124 : END DO
1125 : END DO
1126 : END DO
1127 2 : CALL auxbas_pw_pool%give_back_pw(cdft_control%hw_rho_total)
1128 : END IF
1129 :
1130 24 : CALL rs_grid_release(rs_rho_all)
1131 :
1132 24 : IF (ALLOCATED(coefficients)) DEALLOCATE (coefficients)
1133 24 : IF (ALLOCATED(is_constraint)) DEALLOCATE (is_constraint)
1134 :
1135 24 : CALL timestop(handle)
1136 :
1137 72 : END SUBROUTINE hirshfeld_constraint_low
1138 :
1139 : ! **************************************************************************************************
1140 : !> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT
1141 : !> weight function and the realspace electron density
1142 : !> \param qs_env ...
1143 : ! **************************************************************************************************
1144 3034 : SUBROUTINE cdft_constraint_integrate(qs_env)
1145 : TYPE(qs_environment_type), POINTER :: qs_env
1146 :
1147 : CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_integrate'
1148 :
1149 : INTEGER :: handle, i, iatom, igroup, ikind, ivar, &
1150 : iw, jatom, natom, nvar
1151 : LOGICAL :: is_becke, paw_atom
1152 : REAL(kind=dp) :: dvol, eps_cavity, sign
1153 3034 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dE, strength, target_val
1154 3034 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge, gapw_offset
1155 : TYPE(becke_constraint_type), POINTER :: becke_control
1156 : TYPE(cdft_control_type), POINTER :: cdft_control
1157 3034 : TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1158 : TYPE(cp_logger_type), POINTER :: logger
1159 : TYPE(dft_control_type), POINTER :: dft_control
1160 : TYPE(mp_para_env_type), POINTER :: para_env
1161 3034 : TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho
1162 3034 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1163 3034 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: charge, rho_r
1164 : TYPE(qs_energy_type), POINTER :: energy
1165 3034 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1166 : TYPE(qs_rho_type), POINTER :: rho
1167 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
1168 : TYPE(section_vals_type), POINTER :: cdft_constraint_section
1169 :
1170 3034 : NULLIFY (para_env, dft_control, particle_set, rho_r, energy, rho, &
1171 3034 : logger, cdft_constraint_section, qs_kind_set, mp_rho, &
1172 3034 : rho0_mpole, group, charge)
1173 3034 : CALL timeset(routineN, handle)
1174 3034 : logger => cp_get_default_logger()
1175 : CALL get_qs_env(qs_env, &
1176 : particle_set=particle_set, &
1177 : rho=rho, &
1178 : natom=natom, &
1179 : dft_control=dft_control, &
1180 : para_env=para_env, &
1181 3034 : qs_kind_set=qs_kind_set)
1182 3034 : CALL qs_rho_get(rho, rho_r=rho_r)
1183 3034 : CPASSERT(ASSOCIATED(qs_kind_set))
1184 3034 : cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1185 3034 : iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
1186 3034 : cdft_control => dft_control%qs_control%cdft_control
1187 3034 : is_becke = (cdft_control%type == outer_scf_becke_constraint)
1188 3034 : becke_control => cdft_control%becke_control
1189 3034 : IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
1190 0 : CPABORT("Becke control has not been allocated.")
1191 3034 : group => cdft_control%group
1192 : ! Initialize
1193 3034 : nvar = SIZE(cdft_control%target)
1194 9102 : ALLOCATE (strength(nvar))
1195 6068 : ALLOCATE (target_val(nvar))
1196 6068 : ALLOCATE (dE(nvar))
1197 7044 : strength(:) = cdft_control%strength(:)
1198 7044 : target_val(:) = cdft_control%target(:)
1199 7044 : sign = 1.0_dp
1200 7044 : dE = 0.0_dp
1201 3034 : dvol = group(1)%weight%pw_grid%dvol
1202 3034 : IF (cdft_control%atomic_charges) THEN
1203 1598 : charge => cdft_control%charge
1204 6392 : ALLOCATE (electronic_charge(cdft_control%natoms, dft_control%nspins))
1205 11018 : electronic_charge = 0.0_dp
1206 : END IF
1207 : ! Calculate value of constraint i.e. int ( rho(r) w(r) dr)
1208 9102 : DO i = 1, dft_control%nspins
1209 14088 : DO igroup = 1, SIZE(group)
1210 8020 : SELECT CASE (group(igroup)%constraint_type)
1211 : CASE (cdft_charge_constraint)
1212 16 : sign = 1.0_dp
1213 : CASE (cdft_magnetization_constraint)
1214 16 : IF (i == 1) THEN
1215 : sign = 1.0_dp
1216 : ELSE
1217 8 : sign = -1.0_dp
1218 : END IF
1219 : CASE (cdft_alpha_constraint)
1220 1944 : sign = 1.0_dp
1221 1944 : IF (i == 2) CYCLE
1222 : CASE (cdft_beta_constraint)
1223 1944 : sign = 1.0_dp
1224 1944 : IF (i == 1) CYCLE
1225 : CASE DEFAULT
1226 8020 : CPABORT("Unknown constraint type.")
1227 : END SELECT
1228 12144 : IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1229 : ! With external control, we can use cavity_mat as a mask to kahan sum
1230 180 : eps_cavity = becke_control%eps_cavity
1231 180 : IF (igroup /= 1) &
1232 : CALL cp_abort(__LOCATION__, &
1233 0 : "Multiple constraints not yet supported by parallel mixed calculations.")
1234 : dE(igroup) = dE(igroup) + sign*accurate_dot_product(group(igroup)%weight%array, rho_r(i)%array, &
1235 180 : becke_control%cavity_mat, eps_cavity)*dvol
1236 : ELSE
1237 5896 : dE(igroup) = dE(igroup) + sign*pw_integral_ab(group(igroup)%weight, rho_r(i), local_only=.TRUE.)
1238 : END IF
1239 : END DO
1240 9102 : IF (cdft_control%atomic_charges) THEN
1241 9420 : DO iatom = 1, cdft_control%natoms
1242 9420 : electronic_charge(iatom, i) = pw_integral_ab(charge(iatom), rho_r(i), local_only=.TRUE.)
1243 : END DO
1244 : END IF
1245 : END DO
1246 3034 : CALL get_qs_env(qs_env, energy=energy)
1247 3034 : CALL para_env%sum(dE)
1248 3034 : IF (cdft_control%atomic_charges) THEN
1249 1598 : CALL para_env%sum(electronic_charge)
1250 : END IF
1251 : ! Use fragment densities as reference value (= Becke deformation density)
1252 3034 : IF (cdft_control%fragment_density .AND. .NOT. cdft_control%fragments_integrated) THEN
1253 10 : CALL prepare_fragment_constraint(qs_env)
1254 : END IF
1255 3034 : IF (dft_control%qs_control%gapw) THEN
1256 : ! GAPW: add core charges (rho_hard - rho_soft)
1257 46 : IF (cdft_control%fragment_density) &
1258 : CALL cp_abort(__LOCATION__, &
1259 0 : "Fragment constraints not yet compatible with GAPW.")
1260 184 : ALLOCATE (gapw_offset(nvar, dft_control%nspins))
1261 230 : gapw_offset = 0.0_dp
1262 46 : CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole)
1263 46 : CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho)
1264 138 : DO i = 1, dft_control%nspins
1265 230 : DO igroup = 1, SIZE(group)
1266 460 : DO iatom = 1, SIZE(group(igroup)%atoms)
1267 276 : SELECT CASE (group(igroup)%constraint_type)
1268 : CASE (cdft_charge_constraint)
1269 0 : sign = 1.0_dp
1270 : CASE (cdft_magnetization_constraint)
1271 0 : IF (i == 1) THEN
1272 : sign = 1.0_dp
1273 : ELSE
1274 0 : sign = -1.0_dp
1275 : END IF
1276 : CASE (cdft_alpha_constraint)
1277 0 : sign = 1.0_dp
1278 0 : IF (i == 2) CYCLE
1279 : CASE (cdft_beta_constraint)
1280 0 : sign = 1.0_dp
1281 0 : IF (i == 1) CYCLE
1282 : CASE DEFAULT
1283 276 : CPABORT("Unknown constraint type.")
1284 : END SELECT
1285 276 : jatom = group(igroup)%atoms(iatom)
1286 276 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
1287 276 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
1288 368 : IF (paw_atom) THEN
1289 276 : gapw_offset(igroup, i) = gapw_offset(igroup, i) + sign*group(igroup)%coeff(iatom)*mp_rho(jatom)%q0(i)
1290 : END IF
1291 : END DO
1292 : END DO
1293 : END DO
1294 46 : IF (cdft_control%atomic_charges) THEN
1295 184 : DO iatom = 1, cdft_control%natoms
1296 138 : jatom = cdft_control%atoms(iatom)
1297 138 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
1298 138 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
1299 184 : IF (paw_atom) THEN
1300 414 : DO i = 1, dft_control%nspins
1301 414 : electronic_charge(iatom, i) = electronic_charge(iatom, i) + mp_rho(jatom)%q0(i)
1302 : END DO
1303 : END IF
1304 : END DO
1305 : END IF
1306 138 : DO i = 1, dft_control%nspins
1307 230 : DO ivar = 1, nvar
1308 184 : dE(ivar) = dE(ivar) + gapw_offset(ivar, i)
1309 : END DO
1310 : END DO
1311 46 : DEALLOCATE (gapw_offset)
1312 : END IF
1313 : ! Update constraint value and energy
1314 7044 : cdft_control%value(:) = dE(:)
1315 3034 : energy%cdft = 0.0_dp
1316 7044 : DO ivar = 1, nvar
1317 7044 : energy%cdft = energy%cdft + (dE(ivar) - target_val(ivar))*strength(ivar)
1318 : END DO
1319 : ! Print constraint info and atomic CDFT charges
1320 3034 : CALL cdft_constraint_print(qs_env, electronic_charge)
1321 : ! Deallocate tmp storage
1322 3034 : DEALLOCATE (dE, strength, target_val)
1323 3034 : IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge)
1324 3034 : CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO")
1325 3034 : CALL timestop(handle)
1326 :
1327 6068 : END SUBROUTINE cdft_constraint_integrate
1328 :
1329 : ! **************************************************************************************************
1330 : !> \brief Calculates atomic forces due to a CDFT constraint (Becke or Hirshfeld)
1331 : !> \param qs_env ...
1332 : ! **************************************************************************************************
1333 118 : SUBROUTINE cdft_constraint_force(qs_env)
1334 : TYPE(qs_environment_type), POINTER :: qs_env
1335 :
1336 : CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_force'
1337 :
1338 : INTEGER :: handle, i, iatom, igroup, ikind, ispin, &
1339 : j, k, natom, nvar
1340 118 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
1341 : INTEGER, DIMENSION(2, 3) :: bo
1342 : INTEGER, DIMENSION(3) :: lb, ub
1343 : REAL(kind=dp) :: dvol, eps_cavity, sign
1344 118 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength
1345 118 : REAL(KIND=dp), DIMENSION(:), POINTER :: cutoffs
1346 118 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1347 : TYPE(becke_constraint_type), POINTER :: becke_control
1348 : TYPE(cdft_control_type), POINTER :: cdft_control
1349 118 : TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1350 : TYPE(cell_type), POINTER :: cell
1351 : TYPE(dft_control_type), POINTER :: dft_control
1352 : TYPE(mp_para_env_type), POINTER :: para_env
1353 118 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1354 118 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1355 118 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1356 : TYPE(qs_rho_type), POINTER :: rho
1357 :
1358 118 : CALL timeset(routineN, handle)
1359 118 : NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, &
1360 118 : rho, rho_r, force, cutoffs, becke_control, group)
1361 :
1362 : CALL get_qs_env(qs_env, &
1363 : atomic_kind_set=atomic_kind_set, &
1364 : natom=natom, &
1365 : particle_set=particle_set, &
1366 : cell=cell, &
1367 : rho=rho, &
1368 : force=force, &
1369 : dft_control=dft_control, &
1370 118 : para_env=para_env)
1371 118 : CALL qs_rho_get(rho, rho_r=rho_r)
1372 :
1373 118 : cdft_control => dft_control%qs_control%cdft_control
1374 118 : becke_control => cdft_control%becke_control
1375 118 : group => cdft_control%group
1376 118 : nvar = SIZE(cdft_control%target)
1377 354 : ALLOCATE (strength(nvar))
1378 252 : strength(:) = cdft_control%strength(:)
1379 118 : cutoffs => cdft_control%becke_control%cutoffs
1380 118 : eps_cavity = cdft_control%becke_control%eps_cavity
1381 :
1382 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
1383 : atom_of_kind=atom_of_kind, &
1384 118 : kind_of=kind_of)
1385 252 : DO igroup = 1, SIZE(cdft_control%group)
1386 402 : ALLOCATE (cdft_control%group(igroup)%integrated(3, natom))
1387 1324 : cdft_control%group(igroup)%integrated = 0.0_dp
1388 : END DO
1389 :
1390 472 : lb(1:3) = rho_r(1)%pw_grid%bounds_local(1, 1:3)
1391 472 : ub(1:3) = rho_r(1)%pw_grid%bounds_local(2, 1:3)
1392 1180 : bo = cdft_control%group(1)%weight%pw_grid%bounds_local
1393 118 : dvol = cdft_control%group(1)%weight%pw_grid%dvol
1394 118 : sign = 1.0_dp
1395 :
1396 118 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
1397 116 : IF (.NOT. cdft_control%becke_control%in_memory) THEN
1398 10 : CALL becke_constraint_low(qs_env, just_gradients=.TRUE.)
1399 : END IF
1400 :
1401 2 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1402 2 : IF (.NOT. cdft_control%in_memory) THEN
1403 2 : CALL hirshfeld_constraint_low(qs_env, just_gradients=.TRUE.)
1404 : END IF
1405 : END IF
1406 :
1407 : ! If no Becke Gaussian confinement
1408 118 : IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN
1409 : ! No external control
1410 2106 : DO k = bo(1, 1), bo(2, 1)
1411 93674 : DO j = bo(1, 2), bo(2, 2)
1412 4518732 : DO i = bo(1, 3), bo(2, 3)
1413 : ! First check if this grid point should be skipped
1414 4425152 : IF (cdft_control%becke_control%cavity_confine) THEN
1415 4273152 : IF (cdft_control%becke_control%cavity%array(k, j, i) < eps_cavity) CYCLE
1416 : END IF
1417 :
1418 2713222 : DO igroup = 1, SIZE(cdft_control%group)
1419 8512463 : DO iatom = 1, natom
1420 9537059 : DO ispin = 1, dft_control%nspins
1421 :
1422 5449748 : SELECT CASE (cdft_control%group(igroup)%constraint_type)
1423 : CASE (cdft_charge_constraint)
1424 0 : sign = 1.0_dp
1425 : CASE (cdft_magnetization_constraint)
1426 0 : IF (ispin == 1) THEN
1427 : sign = 1.0_dp
1428 : ELSE
1429 0 : sign = -1.0_dp
1430 : END IF
1431 : CASE (cdft_alpha_constraint)
1432 412880 : sign = 1.0_dp
1433 412880 : IF (ispin == 2) CYCLE
1434 : CASE (cdft_beta_constraint)
1435 412880 : sign = 1.0_dp
1436 412880 : IF (ispin == 1) CYCLE
1437 : CASE DEFAULT
1438 5449748 : CPABORT("Unknown constraint type.")
1439 : END SELECT
1440 :
1441 7761742 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
1442 :
1443 : cdft_control%group(igroup)%integrated(:, iatom) = &
1444 : cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1445 : cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1446 : *rho_r(ispin)%array(k, j, i) &
1447 19123472 : *dvol
1448 :
1449 256000 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1450 :
1451 : cdft_control%group(igroup)%integrated(1, iatom) = &
1452 : cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1453 : cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1454 : *rho_r(ispin)%array(k, j, i) &
1455 256000 : *dvol
1456 :
1457 : cdft_control%group(igroup)%integrated(2, iatom) = &
1458 : cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1459 : cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1460 : *rho_r(ispin)%array(k, j, i) &
1461 256000 : *dvol
1462 :
1463 : cdft_control%group(igroup)%integrated(3, iatom) = &
1464 : cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1465 : cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1466 : *rho_r(ispin)%array(k, j, i) &
1467 256000 : *dvol
1468 :
1469 : END IF
1470 :
1471 : END DO
1472 : END DO
1473 : END DO
1474 : END DO
1475 : END DO
1476 : END DO
1477 :
1478 : ! If Becke Gaussian confinement
1479 : ELSE
1480 1224 : DO k = LBOUND(cdft_control%becke_control%cavity_mat, 1), UBOUND(cdft_control%becke_control%cavity_mat, 1)
1481 61848 : DO j = LBOUND(cdft_control%becke_control%cavity_mat, 2), UBOUND(cdft_control%becke_control%cavity_mat, 2)
1482 2969728 : DO i = LBOUND(cdft_control%becke_control%cavity_mat, 3), UBOUND(cdft_control%becke_control%cavity_mat, 3)
1483 :
1484 : ! First check if this grid point should be skipped
1485 2793472 : IF (cdft_control%becke_control%cavity_mat(k, j, i) < eps_cavity) CYCLE
1486 :
1487 1747632 : DO igroup = 1, SIZE(group)
1488 5327368 : DO iatom = 1, natom
1489 5912424 : DO ispin = 1, dft_control%nspins
1490 3378528 : SELECT CASE (group(igroup)%constraint_type)
1491 : CASE (cdft_charge_constraint)
1492 0 : sign = 1.0_dp
1493 : CASE (cdft_magnetization_constraint)
1494 0 : IF (ispin == 1) THEN
1495 : sign = 1.0_dp
1496 : ELSE
1497 0 : sign = -1.0_dp
1498 : END IF
1499 : CASE (cdft_alpha_constraint)
1500 0 : sign = 1.0_dp
1501 0 : IF (ispin == 2) CYCLE
1502 : CASE (cdft_beta_constraint)
1503 0 : sign = 1.0_dp
1504 0 : IF (ispin == 1) CYCLE
1505 : CASE DEFAULT
1506 3378528 : CPABORT("Unknown constraint type.")
1507 : END SELECT
1508 :
1509 : ! Integrate gradient of weight function
1510 5067792 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
1511 :
1512 : cdft_control%group(igroup)%integrated(:, iatom) = &
1513 : cdft_control%group(igroup)%integrated(:, iatom) + sign* &
1514 : cdft_control%group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) &
1515 : *rho_r(ispin)%array(k, j, i) &
1516 13514112 : *dvol
1517 :
1518 0 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1519 :
1520 : cdft_control%group(igroup)%integrated(1, iatom) = &
1521 : cdft_control%group(igroup)%integrated(1, iatom) + sign* &
1522 : cdft_control%group(igroup)%gradients_x(iatom, k, j, i) &
1523 : *rho_r(ispin)%array(k, j, i) &
1524 0 : *dvol
1525 :
1526 : cdft_control%group(igroup)%integrated(2, iatom) = &
1527 : cdft_control%group(igroup)%integrated(2, iatom) + sign* &
1528 : cdft_control%group(igroup)%gradients_y(iatom, k, j, i) &
1529 : *rho_r(ispin)%array(k, j, i) &
1530 0 : *dvol
1531 :
1532 : cdft_control%group(igroup)%integrated(3, iatom) = &
1533 : cdft_control%group(igroup)%integrated(3, iatom) + sign* &
1534 : cdft_control%group(igroup)%gradients_z(iatom, k, j, i) &
1535 : *rho_r(ispin)%array(k, j, i) &
1536 0 : *dvol
1537 :
1538 : END IF
1539 :
1540 : END DO
1541 : END DO
1542 : END DO
1543 : END DO
1544 : END DO
1545 : END DO
1546 : END IF
1547 :
1548 118 : IF (.NOT. cdft_control%transfer_pot) THEN
1549 98 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
1550 208 : DO igroup = 1, SIZE(group)
1551 208 : DEALLOCATE (cdft_control%group(igroup)%gradients)
1552 : END DO
1553 2 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1554 4 : DO igroup = 1, SIZE(group)
1555 2 : DEALLOCATE (cdft_control%group(igroup)%gradients_x)
1556 2 : DEALLOCATE (cdft_control%group(igroup)%gradients_y)
1557 4 : DEALLOCATE (cdft_control%group(igroup)%gradients_z)
1558 : END DO
1559 : END IF
1560 : END IF
1561 :
1562 252 : DO igroup = 1, SIZE(group)
1563 2396 : CALL para_env%sum(group(igroup)%integrated)
1564 : END DO
1565 :
1566 : ! Update force only on master process. Otherwise force due to constraint becomes multiplied
1567 : ! by the number of processes when the final force%rho_elec is constructed in qs_force
1568 : ! by mp_summing [the final integrated(:,:) is distributed on all processors]
1569 118 : IF (para_env%is_source()) THEN
1570 150 : DO igroup = 1, SIZE(group)
1571 308 : DO iatom = 1, natom
1572 158 : ikind = kind_of(iatom)
1573 158 : i = atom_of_kind(iatom)
1574 1185 : force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup)
1575 : END DO
1576 : END DO
1577 : END IF
1578 :
1579 118 : DEALLOCATE (strength)
1580 252 : DO igroup = 1, SIZE(group)
1581 252 : DEALLOCATE (group(igroup)%integrated)
1582 : END DO
1583 118 : NULLIFY (group)
1584 :
1585 118 : CALL timestop(handle)
1586 :
1587 236 : END SUBROUTINE cdft_constraint_force
1588 :
1589 : ! **************************************************************************************************
1590 : !> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied
1591 : !> by the CDFT weight functions and integrated over the realspace grid.
1592 : !> \param qs_env ...
1593 : ! **************************************************************************************************
1594 10 : SUBROUTINE prepare_fragment_constraint(qs_env)
1595 : TYPE(qs_environment_type), POINTER :: qs_env
1596 :
1597 : CHARACTER(len=*), PARAMETER :: routineN = 'prepare_fragment_constraint'
1598 :
1599 : INTEGER :: handle, i, iatom, igroup, natom, &
1600 : nelectron_total, nfrag_spins
1601 : LOGICAL :: is_becke, needs_spin_density
1602 : REAL(kind=dp) :: dvol, multiplier(2), nelectron_frag
1603 : TYPE(becke_constraint_type), POINTER :: becke_control
1604 : TYPE(cdft_control_type), POINTER :: cdft_control
1605 10 : TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1606 : TYPE(cp_logger_type), POINTER :: logger
1607 : TYPE(dft_control_type), POINTER :: dft_control
1608 : TYPE(mp_para_env_type), POINTER :: para_env
1609 : TYPE(pw_env_type), POINTER :: pw_env
1610 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1611 10 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: rho_frag
1612 : TYPE(qs_subsys_type), POINTER :: subsys
1613 :
1614 10 : NULLIFY (para_env, dft_control, logger, subsys, pw_env, auxbas_pw_pool, group)
1615 10 : CALL timeset(routineN, handle)
1616 10 : logger => cp_get_default_logger()
1617 : CALL get_qs_env(qs_env, &
1618 : natom=natom, &
1619 : dft_control=dft_control, &
1620 10 : para_env=para_env)
1621 :
1622 10 : cdft_control => dft_control%qs_control%cdft_control
1623 10 : is_becke = (cdft_control%type == outer_scf_becke_constraint)
1624 10 : becke_control => cdft_control%becke_control
1625 10 : IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) &
1626 0 : CPABORT("Becke control has not been allocated.")
1627 10 : group => cdft_control%group
1628 10 : dvol = group(1)%weight%pw_grid%dvol
1629 : ! Fragment densities are meaningful only for some calculation types
1630 10 : IF (.NOT. qs_env%single_point_run) &
1631 : CALL cp_abort(__LOCATION__, &
1632 : "CDFT fragment constraints are only compatible with single "// &
1633 0 : "point calculations (run_type ENERGY or ENERGY_FORCE).")
1634 10 : IF (dft_control%qs_control%gapw) &
1635 : CALL cp_abort(__LOCATION__, &
1636 0 : "CDFT fragment constraint not compatible with GAPW.")
1637 30 : needs_spin_density = .FALSE.
1638 30 : multiplier = 1.0_dp
1639 10 : nfrag_spins = 1
1640 22 : DO igroup = 1, SIZE(group)
1641 10 : SELECT CASE (group(igroup)%constraint_type)
1642 : CASE (cdft_charge_constraint)
1643 : ! Do nothing
1644 : CASE (cdft_magnetization_constraint)
1645 6 : needs_spin_density = .TRUE.
1646 : CASE (cdft_alpha_constraint, cdft_beta_constraint)
1647 : CALL cp_abort(__LOCATION__, &
1648 : "CDFT fragment constraint not yet compatible with "// &
1649 0 : "spin specific constraints.")
1650 : CASE DEFAULT
1651 12 : CPABORT("Unknown constraint type.")
1652 : END SELECT
1653 : END DO
1654 10 : IF (needs_spin_density) THEN
1655 12 : nfrag_spins = 2
1656 12 : DO i = 1, 2
1657 12 : IF (cdft_control%flip_fragment(i)) multiplier(i) = -1.0_dp
1658 : END DO
1659 : END IF
1660 : ! Read fragment reference densities
1661 68 : ALLOCATE (cdft_control%fragments(nfrag_spins, 2))
1662 34 : ALLOCATE (rho_frag(nfrag_spins))
1663 10 : CALL get_qs_env(qs_env, pw_env=pw_env)
1664 10 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1665 : ! Total density (rho_alpha + rho_beta)
1666 10 : CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 1))
1667 : CALL cp_cube_to_pw(cdft_control%fragments(1, 1), &
1668 10 : cdft_control%fragment_a_fname, 1.0_dp)
1669 10 : CALL auxbas_pw_pool%create_pw(cdft_control%fragments(1, 2))
1670 : CALL cp_cube_to_pw(cdft_control%fragments(1, 2), &
1671 10 : cdft_control%fragment_b_fname, 1.0_dp)
1672 : ! Spin difference density (rho_alpha - rho_beta) if needed
1673 10 : IF (needs_spin_density) THEN
1674 4 : CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 1))
1675 : CALL cp_cube_to_pw(cdft_control%fragments(2, 1), &
1676 4 : cdft_control%fragment_a_spin_fname, multiplier(1))
1677 4 : CALL auxbas_pw_pool%create_pw(cdft_control%fragments(2, 2))
1678 : CALL cp_cube_to_pw(cdft_control%fragments(2, 2), &
1679 4 : cdft_control%fragment_b_spin_fname, multiplier(2))
1680 : END IF
1681 : ! Sum up fragments
1682 24 : DO i = 1, nfrag_spins
1683 14 : CALL auxbas_pw_pool%create_pw(rho_frag(i))
1684 14 : CALL pw_copy(cdft_control%fragments(i, 1), rho_frag(i))
1685 14 : CALL pw_axpy(cdft_control%fragments(i, 2), rho_frag(i), 1.0_dp)
1686 14 : CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 1))
1687 24 : CALL auxbas_pw_pool%give_back_pw(cdft_control%fragments(i, 2))
1688 : END DO
1689 10 : DEALLOCATE (cdft_control%fragments)
1690 : ! Check that the number of electrons is consistent
1691 10 : CALL get_qs_env(qs_env, subsys=subsys)
1692 10 : CALL qs_subsys_get(subsys, nelectron_total=nelectron_total)
1693 10 : nelectron_frag = pw_integrate_function(rho_frag(1))
1694 10 : IF (NINT(nelectron_frag) /= nelectron_total) &
1695 : CALL cp_abort(__LOCATION__, &
1696 : "The number of electrons in the reference and interacting "// &
1697 0 : "configurations does not match. Check your fragment cube files.")
1698 : ! Update constraint target value i.e. perform integration w_i*rho_frag_{tot/spin}*dr
1699 22 : cdft_control%target = 0.0_dp
1700 22 : DO igroup = 1, SIZE(group)
1701 12 : IF (group(igroup)%constraint_type == cdft_charge_constraint) THEN
1702 : i = 1
1703 : ELSE
1704 6 : i = 2
1705 : END IF
1706 22 : IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN
1707 : cdft_control%target(igroup) = cdft_control%target(igroup) + &
1708 : accurate_dot_product(group(igroup)%weight%array, rho_frag(i)%array, &
1709 0 : becke_control%cavity_mat, becke_control%eps_cavity)*dvol
1710 : ELSE
1711 : cdft_control%target(igroup) = cdft_control%target(igroup) + &
1712 12 : pw_integral_ab(group(igroup)%weight, rho_frag(i), local_only=.TRUE.)
1713 : END IF
1714 : END DO
1715 34 : CALL para_env%sum(cdft_control%target)
1716 : ! Calculate reference atomic charges int( w_i * rho_frag * dr )
1717 10 : IF (cdft_control%atomic_charges) THEN
1718 40 : ALLOCATE (cdft_control%charges_fragment(cdft_control%natoms, nfrag_spins))
1719 24 : DO i = 1, nfrag_spins
1720 44 : DO iatom = 1, cdft_control%natoms
1721 : cdft_control%charges_fragment(iatom, i) = &
1722 34 : pw_integral_ab(cdft_control%charge(iatom), rho_frag(i), local_only=.TRUE.)
1723 : END DO
1724 : END DO
1725 78 : CALL para_env%sum(cdft_control%charges_fragment)
1726 : END IF
1727 24 : DO i = 1, nfrag_spins
1728 24 : CALL auxbas_pw_pool%give_back_pw(rho_frag(i))
1729 : END DO
1730 10 : DEALLOCATE (rho_frag)
1731 10 : cdft_control%fragments_integrated = .TRUE.
1732 :
1733 10 : CALL timestop(handle)
1734 :
1735 10 : END SUBROUTINE prepare_fragment_constraint
1736 :
1737 : END MODULE qs_cdft_methods
|
