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 Pointwise CDFT partition functions for nonuniform integration grids.
10 : ! **************************************************************************************************
11 : MODULE qs_cdft_grid
12 : USE atomic_kind_types, ONLY: atomic_kind_type,&
13 : get_atomic_kind
14 : USE cell_types, ONLY: cell_type,&
15 : pbc
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_units, ONLY: cp_unit_from_cp2k
18 : USE hirshfeld_methods, ONLY: create_shape_function
19 : USE hirshfeld_types, ONLY: hirshfeld_type
20 : USE input_constants, ONLY: becke_cutoff_element,&
21 : becke_cutoff_global,&
22 : outer_scf_becke_constraint,&
23 : outer_scf_hirshfeld_constraint
24 : USE kinds, ONLY: dp
25 : USE particle_types, ONLY: particle_type
26 : USE qs_cdft_types, ONLY: becke_constraint_type,&
27 : cdft_control_type,&
28 : hirshfeld_constraint_type
29 : USE qs_cdft_utils, ONLY: hirshfeld_constraint_init
30 : USE qs_environment_types, ONLY: get_qs_env,&
31 : qs_environment_type
32 : USE qs_kind_types, ONLY: qs_kind_type
33 : #include "./base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 :
37 : PRIVATE
38 :
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_grid'
40 :
41 : TYPE cdft_point_context_type
42 : INTEGER :: method = -1, natom = 0, ngroup = 0
43 : INTEGER, ALLOCATABLE, DIMENSION(:) :: numexp, cavity_numexp
44 : LOGICAL :: calculate_derivatives = .FALSE., &
45 : cavity_confine = .FALSE.
46 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: constraint_atom
47 : REAL(KIND=dp) :: eps = 0.0_dp, eps_cavity = 0.0_dp
48 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cutoffs, distances, cell_functions
49 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atom_coords, aij, coefficients, &
50 : alpha, amplitude, cavity_alpha, &
51 : cavity_amplitude, displacement, &
52 : datom_numerator, datom_sum, dcell_point, &
53 : density_atom_derivative
54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pair_vectors, dcell_atom
55 : TYPE(cell_type), POINTER :: cell => NULL()
56 : END TYPE cdft_point_context_type
57 :
58 : PUBLIC :: cdft_point_context_create, &
59 : cdft_point_context_release, &
60 : cdft_point_context_type, &
61 : cdft_point_weights
62 :
63 : CONTAINS
64 :
65 : ! **************************************************************************************************
66 : !> \brief Initialize reusable data for pointwise CDFT partition evaluation.
67 : !> \param qs_env Quickstep environment
68 : !> \param context pointwise partition context
69 : !> \param calculate_derivatives allocate scratch space for coordinate derivatives
70 : ! **************************************************************************************************
71 424 : SUBROUTINE cdft_point_context_create(qs_env, context, calculate_derivatives)
72 : TYPE(qs_environment_type), POINTER :: qs_env
73 : TYPE(cdft_point_context_type), INTENT(OUT) :: context
74 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_derivatives
75 :
76 : INTEGER :: atom, iatom, igroup, ikind, jatom, natom
77 : REAL(KIND=dp) :: chi, ircov, jrcov, uij
78 : REAL(KIND=dp), DIMENSION(3) :: pair_vector
79 424 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii, radii_list
80 424 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
81 : TYPE(becke_constraint_type), POINTER :: becke_control
82 : TYPE(cdft_control_type), POINTER :: cdft_control
83 : TYPE(dft_control_type), POINTER :: dft_control
84 : TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control
85 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
86 424 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
87 424 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
88 :
89 424 : NULLIFY (becke_control, cdft_control, dft_control, hirshfeld_control, &
90 424 : hirshfeld_env, atomic_kind_set, particle_set, qs_kind_set, radii, radii_list)
91 : CALL get_qs_env(qs_env, cell=context%cell, dft_control=dft_control, &
92 : natom=natom, particle_set=particle_set, atomic_kind_set=atomic_kind_set, &
93 424 : qs_kind_set=qs_kind_set)
94 424 : CPASSERT(ASSOCIATED(context%cell))
95 424 : CPASSERT(ASSOCIATED(dft_control))
96 424 : CPASSERT(ASSOCIATED(particle_set))
97 424 : CPASSERT(ASSOCIATED(atomic_kind_set))
98 424 : CPASSERT(ASSOCIATED(qs_kind_set))
99 424 : cdft_control => dft_control%qs_control%cdft_control
100 424 : CPASSERT(ASSOCIATED(cdft_control))
101 :
102 424 : context%method = cdft_control%type
103 424 : context%natom = natom
104 424 : context%ngroup = SIZE(cdft_control%group)
105 424 : context%calculate_derivatives = .FALSE.
106 424 : IF (PRESENT(calculate_derivatives)) context%calculate_derivatives = calculate_derivatives
107 0 : ALLOCATE (context%atom_coords(3, natom), &
108 0 : context%coefficients(context%ngroup, natom), &
109 0 : context%constraint_atom(natom), &
110 0 : context%distances(natom), context%displacement(3, natom), &
111 5088 : context%cell_functions(natom))
112 2484 : context%coefficients = 0.0_dp
113 1382 : context%constraint_atom = .FALSE.
114 1382 : DO atom = 1, natom
115 4256 : context%atom_coords(:, atom) = particle_set(atom)%r
116 : END DO
117 920 : DO igroup = 1, context%ngroup
118 1804 : DO iatom = 1, SIZE(cdft_control%group(igroup)%atoms)
119 884 : atom = cdft_control%group(igroup)%atoms(iatom)
120 884 : context%coefficients(igroup, atom) = cdft_control%group(igroup)%coeff(iatom)
121 1380 : context%constraint_atom(atom) = .TRUE.
122 : END DO
123 : END DO
124 :
125 424 : SELECT CASE (context%method)
126 : CASE (outer_scf_becke_constraint)
127 334 : becke_control => cdft_control%becke_control
128 334 : CPASSERT(ASSOCIATED(becke_control))
129 0 : ALLOCATE (context%cutoffs(natom), context%aij(natom, natom), &
130 3006 : context%pair_vectors(3, natom, natom))
131 334 : IF (context%calculate_derivatives) THEN
132 0 : ALLOCATE (context%dcell_point(3, natom), context%dcell_atom(3, natom, natom), &
133 28 : context%datom_sum(3, natom), context%datom_numerator(3, natom))
134 : END IF
135 334 : IF (ASSOCIATED(becke_control%cutoffs)) THEN
136 890 : context%cutoffs(:) = becke_control%cutoffs
137 : ELSE
138 148 : SELECT CASE (becke_control%cutoff_type)
139 : CASE (becke_cutoff_global)
140 222 : context%cutoffs = becke_control%rglobal
141 : CASE (becke_cutoff_element)
142 0 : CPASSERT(ASSOCIATED(becke_control%cutoffs_tmp))
143 0 : CPASSERT(SIZE(becke_control%cutoffs_tmp) == SIZE(atomic_kind_set))
144 0 : DO atom = 1, natom
145 0 : CALL get_atomic_kind(particle_set(atom)%atomic_kind, kind_number=ikind)
146 0 : context%cutoffs(atom) = becke_control%cutoffs_tmp(ikind)
147 : END DO
148 : CASE DEFAULT
149 74 : CPABORT("Unknown Becke cutoff type.")
150 : END SELECT
151 : END IF
152 2998 : context%aij = 0.0_dp
153 334 : IF (becke_control%adjust) THEN
154 46 : IF (ASSOCIATED(becke_control%aij)) THEN
155 598 : context%aij(:, :) = becke_control%aij
156 : ELSE
157 0 : IF (ASSOCIATED(becke_control%radii)) THEN
158 0 : radii => becke_control%radii
159 : ELSE
160 0 : radii => becke_control%radii_tmp
161 : END IF
162 0 : CPASSERT(ASSOCIATED(radii))
163 0 : CPASSERT(SIZE(radii) == SIZE(atomic_kind_set))
164 0 : DO iatom = 1, natom - 1
165 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
166 0 : ircov = radii(ikind)
167 0 : DO jatom = iatom + 1, natom
168 0 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind)
169 0 : jrcov = radii(ikind)
170 0 : IF (ircov /= jrcov) THEN
171 0 : chi = ircov/jrcov
172 0 : uij = (chi - 1.0_dp)/(chi + 1.0_dp)
173 : context%aij(iatom, jatom) = MAX(-0.5_dp, MIN(0.5_dp, &
174 0 : uij/(uij**2 - 1.0_dp)))
175 0 : context%aij(jatom, iatom) = -context%aij(iatom, jatom)
176 : END IF
177 : END DO
178 : END DO
179 : END IF
180 : END IF
181 8656 : context%pair_vectors = 0.0_dp
182 778 : DO iatom = 1, natom - 1
183 1332 : DO jatom = iatom + 1, natom
184 : pair_vector = pbc(context%atom_coords(:, jatom), &
185 554 : context%atom_coords(:, iatom), context%cell)
186 2216 : context%pair_vectors(:, iatom, jatom) = pair_vector
187 2660 : context%pair_vectors(:, jatom, iatom) = -pair_vector
188 : END DO
189 : END DO
190 334 : context%cavity_confine = becke_control%cavity_confine
191 334 : context%eps_cavity = becke_control%eps_cavity
192 334 : IF (context%cavity_confine) THEN
193 46 : hirshfeld_env => becke_control%cavity_env
194 46 : CPASSERT(ASSOCIATED(hirshfeld_env))
195 46 : IF (.NOT. ASSOCIATED(hirshfeld_env%kind_shape_fn)) THEN
196 0 : IF (ASSOCIATED(becke_control%radii)) THEN
197 0 : radii => becke_control%radii
198 0 : ELSE IF (ASSOCIATED(becke_control%radii_tmp)) THEN
199 0 : radii => becke_control%radii_tmp
200 : END IF
201 0 : IF (ASSOCIATED(radii)) THEN
202 0 : ALLOCATE (radii_list(SIZE(radii)))
203 0 : DO ikind = 1, SIZE(radii)
204 0 : IF (hirshfeld_env%use_bohr) THEN
205 0 : radii_list(ikind) = radii(ikind)
206 : ELSE
207 0 : radii_list(ikind) = cp_unit_from_cp2k(radii(ikind), "angstrom")
208 : END IF
209 : END DO
210 : END IF
211 : CALL create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, &
212 0 : radius=becke_control%rcavity, radii_list=radii_list)
213 0 : IF (ASSOCIATED(radii_list)) DEALLOCATE (radii_list)
214 : END IF
215 : CALL store_shape_functions(hirshfeld_env, particle_set, context%cavity_numexp, &
216 : context%cavity_alpha, context%cavity_amplitude, &
217 46 : include_charge=.FALSE.)
218 : END IF
219 : CASE (outer_scf_hirshfeld_constraint)
220 90 : hirshfeld_control => cdft_control%hirshfeld_control
221 90 : CPASSERT(ASSOCIATED(hirshfeld_control))
222 90 : hirshfeld_env => hirshfeld_control%hirshfeld_env
223 90 : CPASSERT(ASSOCIATED(hirshfeld_env))
224 90 : IF (.NOT. ASSOCIATED(hirshfeld_env%kind_shape_fn) .OR. &
225 : .NOT. ASSOCIATED(hirshfeld_env%charges)) THEN
226 0 : CALL hirshfeld_constraint_init(qs_env)
227 : END IF
228 90 : context%eps = hirshfeld_control%eps_cutoff
229 : CALL store_shape_functions(hirshfeld_env, particle_set, context%numexp, &
230 90 : context%alpha, context%amplitude, include_charge=.TRUE.)
231 90 : IF (context%calculate_derivatives) THEN
232 12 : ALLOCATE (context%density_atom_derivative(3, natom))
233 : END IF
234 : CASE DEFAULT
235 424 : CPABORT("Unknown CDFT partition type.")
236 : END SELECT
237 :
238 : CONTAINS
239 :
240 : ! **************************************************************************************************
241 : !> \brief ...
242 : !> \param environment ...
243 : !> \param particles ...
244 : !> \param numexp ...
245 : !> \param alpha ...
246 : !> \param amplitude ...
247 : !> \param include_charge ...
248 : ! **************************************************************************************************
249 136 : SUBROUTINE store_shape_functions(environment, particles, numexp, alpha, amplitude, &
250 : include_charge)
251 : TYPE(hirshfeld_type), POINTER :: environment
252 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
253 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: numexp
254 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
255 : INTENT(OUT) :: alpha, amplitude
256 : LOGICAL, INTENT(IN) :: include_charge
257 :
258 : INTEGER :: atom, iexp, ikind, maxexp
259 : REAL(KIND=dp) :: charge
260 :
261 408 : ALLOCATE (numexp(natom))
262 136 : numexp = 0
263 136 : maxexp = 0
264 454 : DO atom = 1, natom
265 318 : CALL get_atomic_kind(particles(atom)%atomic_kind, kind_number=ikind)
266 318 : numexp(atom) = environment%kind_shape_fn(ikind)%numexp
267 454 : maxexp = MAX(maxexp, numexp(atom))
268 : END DO
269 816 : ALLOCATE (alpha(maxexp, natom), amplitude(maxexp, natom))
270 136 : alpha = 0.0_dp
271 136 : amplitude = 0.0_dp
272 454 : DO atom = 1, natom
273 318 : CALL get_atomic_kind(particles(atom)%atomic_kind, kind_number=ikind)
274 318 : charge = 1.0_dp
275 318 : IF (include_charge) charge = environment%charges(atom)
276 1090 : DO iexp = 1, numexp(atom)
277 318 : alpha(iexp, atom) = environment%kind_shape_fn(ikind)%zet(iexp)
278 636 : amplitude(iexp, atom) = charge*environment%kind_shape_fn(ikind)%coef(iexp)
279 : END DO
280 : END DO
281 136 : END SUBROUTINE store_shape_functions
282 :
283 : END SUBROUTINE cdft_point_context_create
284 :
285 : ! **************************************************************************************************
286 : !> \brief Release a pointwise CDFT partition context.
287 : !> \param context pointwise partition context
288 : ! **************************************************************************************************
289 424 : SUBROUTINE cdft_point_context_release(context)
290 : TYPE(cdft_point_context_type), INTENT(INOUT) :: context
291 :
292 424 : IF (ALLOCATED(context%numexp)) DEALLOCATE (context%numexp)
293 424 : IF (ALLOCATED(context%cavity_numexp)) DEALLOCATE (context%cavity_numexp)
294 424 : IF (ALLOCATED(context%constraint_atom)) DEALLOCATE (context%constraint_atom)
295 424 : IF (ALLOCATED(context%cutoffs)) DEALLOCATE (context%cutoffs)
296 424 : IF (ALLOCATED(context%distances)) DEALLOCATE (context%distances)
297 424 : IF (ALLOCATED(context%cell_functions)) DEALLOCATE (context%cell_functions)
298 424 : IF (ALLOCATED(context%datom_sum)) DEALLOCATE (context%datom_sum)
299 424 : IF (ALLOCATED(context%atom_coords)) DEALLOCATE (context%atom_coords)
300 424 : IF (ALLOCATED(context%aij)) DEALLOCATE (context%aij)
301 424 : IF (ALLOCATED(context%coefficients)) DEALLOCATE (context%coefficients)
302 424 : IF (ALLOCATED(context%alpha)) DEALLOCATE (context%alpha)
303 424 : IF (ALLOCATED(context%amplitude)) DEALLOCATE (context%amplitude)
304 424 : IF (ALLOCATED(context%cavity_alpha)) DEALLOCATE (context%cavity_alpha)
305 424 : IF (ALLOCATED(context%cavity_amplitude)) DEALLOCATE (context%cavity_amplitude)
306 424 : IF (ALLOCATED(context%displacement)) DEALLOCATE (context%displacement)
307 424 : IF (ALLOCATED(context%datom_numerator)) DEALLOCATE (context%datom_numerator)
308 424 : IF (ALLOCATED(context%dcell_point)) DEALLOCATE (context%dcell_point)
309 424 : IF (ALLOCATED(context%density_atom_derivative)) DEALLOCATE (context%density_atom_derivative)
310 424 : IF (ALLOCATED(context%pair_vectors)) DEALLOCATE (context%pair_vectors)
311 424 : IF (ALLOCATED(context%dcell_atom)) DEALLOCATE (context%dcell_atom)
312 424 : NULLIFY (context%cell)
313 424 : context%method = -1
314 424 : context%natom = 0
315 424 : context%ngroup = 0
316 424 : context%calculate_derivatives = .FALSE.
317 424 : END SUBROUTINE cdft_point_context_release
318 :
319 : ! **************************************************************************************************
320 : !> \brief Evaluate CDFT weights and coordinate derivatives at one point.
321 : !> \param context pointwise partition context
322 : !> \param point Cartesian point
323 : !> \param weights group weights
324 : !> \param point_derivative derivatives with respect to the Cartesian point
325 : !> \param atom_derivative derivatives with respect to atom positions at fixed point
326 : !> \param atomic_weights optional individual atomic partition weights
327 : ! **************************************************************************************************
328 1197500 : SUBROUTINE cdft_point_weights(context, point, weights, point_derivative, atom_derivative, &
329 1197500 : atomic_weights)
330 : TYPE(cdft_point_context_type), INTENT(INOUT) :: context
331 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: point
332 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: weights
333 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: point_derivative
334 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: atom_derivative
335 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: atomic_weights
336 :
337 1197500 : CPASSERT(SIZE(weights) == context%ngroup)
338 1197500 : CPASSERT(SIZE(point_derivative, 1) == 3)
339 1197500 : CPASSERT(SIZE(point_derivative, 2) == context%ngroup)
340 1197500 : CPASSERT(SIZE(atom_derivative, 1) == 3)
341 1197500 : CPASSERT(SIZE(atom_derivative, 2) == context%natom)
342 1197500 : CPASSERT(SIZE(atom_derivative, 3) == context%ngroup)
343 1197500 : IF (PRESENT(atomic_weights)) THEN
344 247500 : CPASSERT(SIZE(atomic_weights) == context%natom)
345 : END IF
346 2170000 : SELECT CASE (context%method)
347 : CASE (outer_scf_becke_constraint)
348 : CALL becke_point_weights(context, point, weights, point_derivative, atom_derivative, &
349 1772500 : atomic_weights, context%calculate_derivatives)
350 : CASE (outer_scf_hirshfeld_constraint)
351 : CALL hirshfeld_point_weights(context, point, weights, point_derivative, atom_derivative, &
352 375000 : atomic_weights, context%calculate_derivatives)
353 : CASE DEFAULT
354 1197500 : CPABORT("Unknown CDFT partition type.")
355 : END SELECT
356 1197500 : END SUBROUTINE cdft_point_weights
357 :
358 : ! **************************************************************************************************
359 : !> \brief Evaluate Becke weights at one point.
360 : !> \param context ...
361 : !> \param point ...
362 : !> \param weights ...
363 : !> \param point_derivative ...
364 : !> \param atom_derivative ...
365 : !> \param atomic_weights ...
366 : !> \param calculate_derivatives ...
367 : ! **************************************************************************************************
368 972500 : SUBROUTINE becke_point_weights(context, point, weights, point_derivative, atom_derivative, &
369 972500 : atomic_weights, calculate_derivatives)
370 : TYPE(cdft_point_context_type), INTENT(INOUT) :: context
371 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: point
372 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: weights
373 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: point_derivative
374 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: atom_derivative
375 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: atomic_weights
376 : LOGICAL, INTENT(IN) :: calculate_derivatives
377 :
378 : INTEGER :: atom, iatom, iexp, igroup, jatom
379 : REAL(KIND=dp) :: adjusted_mu, cavity_density, delta, &
380 : dmu_factor, f1, f2, f3, mu, numerator, &
381 : old_cell, pair_distance, s, sum_cell
382 : REAL(KIND=dp), DIMENSION(3) :: dmu_i, dmu_j, dmu_point, &
383 : dpoint_numerator, ds_i, ds_j, &
384 : ds_point, dsum_point, unit_i, unit_j
385 :
386 2125000 : weights = 0.0_dp
387 972500 : IF (calculate_derivatives) THEN
388 82500 : point_derivative = 0.0_dp
389 200000 : atom_derivative = 0.0_dp
390 : END IF
391 1490000 : IF (PRESENT(atomic_weights)) atomic_weights = 0.0_dp
392 :
393 972500 : IF (context%cavity_confine) THEN
394 172500 : cavity_density = 0.0_dp
395 690000 : DO atom = 1, context%natom
396 517500 : IF (.NOT. context%constraint_atom(atom)) CYCLE
397 2070000 : context%displacement(:, atom) = pbc(context%atom_coords(:, atom), point, context%cell)
398 1207500 : DO iexp = 1, context%cavity_numexp(atom)
399 : cavity_density = cavity_density + context%cavity_amplitude(iexp, atom)* &
400 : EXP(-context%cavity_alpha(iexp, atom)* &
401 2587500 : DOT_PRODUCT(context%displacement(:, atom), context%displacement(:, atom)))
402 : END DO
403 : END DO
404 227055 : IF (cavity_density < context%eps_cavity) RETURN
405 : END IF
406 :
407 3292464 : context%cell_functions = 1.0_dp
408 963116 : IF (calculate_derivatives) THEN
409 142500 : context%dcell_point = 0.0_dp
410 395000 : context%dcell_atom = 0.0_dp
411 : END IF
412 3292464 : DO atom = 1, context%natom
413 9317392 : context%displacement(:, atom) = pbc(context%atom_coords(:, atom), point, context%cell)
414 10280508 : context%distances(atom) = NORM2(context%displacement(:, atom))
415 : END DO
416 3292464 : DO iatom = 1, context%natom
417 2329348 : IF (context%distances(iatom) > context%cutoffs(iatom)) THEN
418 185920 : context%cell_functions(iatom) = 0.0_dp
419 185920 : CYCLE
420 : END IF
421 2143428 : IF (calculate_derivatives) THEN
422 32500 : unit_i = 0.0_dp
423 32500 : IF (context%distances(iatom) > 1.0E-14_dp) THEN
424 130000 : unit_i = -context%displacement(:, iatom)/context%distances(iatom)
425 : END IF
426 : END IF
427 8488944 : DO jatom = 1, context%natom
428 5382400 : IF (jatom == iatom) CYCLE
429 12955888 : pair_distance = NORM2(context%pair_vectors(:, iatom, jatom))
430 3238972 : IF (pair_distance <= 1.0E-14_dp) CYCLE
431 3238972 : delta = context%distances(iatom) - context%distances(jatom)
432 3238972 : mu = delta/pair_distance
433 3238972 : adjusted_mu = mu + context%aij(iatom, jatom)*(1.0_dp - mu**2)
434 3238972 : f1 = 1.5_dp*adjusted_mu - 0.5_dp*adjusted_mu**3
435 3238972 : f2 = 1.5_dp*f1 - 0.5_dp*f1**3
436 3238972 : f3 = 1.5_dp*f2 - 0.5_dp*f2**3
437 3238972 : s = 0.5_dp*(1.0_dp - f3)
438 3238972 : old_cell = context%cell_functions(iatom)
439 3238972 : IF (calculate_derivatives) THEN
440 55000 : dmu_factor = 1.0_dp - 2.0_dp*context%aij(iatom, jatom)*mu
441 55000 : unit_j = 0.0_dp
442 55000 : IF (context%distances(jatom) > 1.0E-14_dp) THEN
443 220000 : unit_j = -context%displacement(:, jatom)/context%distances(jatom)
444 : END IF
445 : dmu_i = unit_i/pair_distance - &
446 220000 : delta*context%pair_vectors(:, iatom, jatom)/pair_distance**3
447 : dmu_j = -unit_j/pair_distance + &
448 220000 : delta*context%pair_vectors(:, iatom, jatom)/pair_distance**3
449 220000 : dmu_point = (-unit_i + unit_j)/pair_distance
450 : dmu_factor = -0.5_dp*dmu_factor*1.5_dp*(1.0_dp - adjusted_mu**2)* &
451 55000 : 1.5_dp*(1.0_dp - f1**2)*1.5_dp*(1.0_dp - f2**2)
452 220000 : ds_i = dmu_factor*dmu_i
453 220000 : ds_j = dmu_factor*dmu_j
454 220000 : ds_point = dmu_factor*dmu_point
455 675000 : context%dcell_atom(:, :, iatom) = context%dcell_atom(:, :, iatom)*s
456 220000 : context%dcell_point(:, iatom) = context%dcell_point(:, iatom)*s
457 : context%dcell_atom(:, iatom, iatom) = &
458 220000 : context%dcell_atom(:, iatom, iatom) + old_cell*ds_i
459 : context%dcell_atom(:, jatom, iatom) = &
460 220000 : context%dcell_atom(:, jatom, iatom) + old_cell*ds_j
461 220000 : context%dcell_point(:, iatom) = context%dcell_point(:, iatom) + old_cell*ds_point
462 : END IF
463 7525828 : context%cell_functions(iatom) = old_cell*s
464 : END DO
465 : END DO
466 :
467 3292464 : sum_cell = SUM(context%cell_functions)
468 963116 : IF (sum_cell <= 1.0E-6_dp) RETURN
469 1318076 : IF (PRESENT(atomic_weights)) atomic_weights = context%cell_functions/sum_cell
470 908561 : IF (calculate_derivatives) THEN
471 147500 : dsum_point = SUM(context%dcell_point, DIM=2)
472 405000 : context%datom_sum(:, :) = SUM(context%dcell_atom, DIM=3)
473 : END IF
474 1997122 : DO igroup = 1, context%ngroup
475 3642188 : numerator = DOT_PRODUCT(context%coefficients(igroup, :), context%cell_functions)
476 1088561 : weights(igroup) = numerator/sum_cell
477 1997122 : IF (calculate_derivatives) THEN
478 442500 : dpoint_numerator = MATMUL(context%dcell_point, context%coefficients(igroup, :))
479 : point_derivative(:, igroup) = &
480 70000 : (dpoint_numerator*sum_cell - numerator*dsum_point)/sum_cell**2
481 60000 : DO atom = 1, context%natom
482 : context%datom_numerator(:, atom) = &
483 1567500 : MATMUL(context%dcell_atom(:, atom, :), context%coefficients(igroup, :))
484 : atom_derivative(:, atom, igroup) = &
485 : (context%datom_numerator(:, atom)*sum_cell - &
486 187500 : numerator*context%datom_sum(:, atom))/sum_cell**2
487 : END DO
488 : END IF
489 : END DO
490 : END SUBROUTINE becke_point_weights
491 :
492 : ! **************************************************************************************************
493 : !> \brief Evaluate Hirshfeld weights at one point.
494 : !> \param context ...
495 : !> \param point ...
496 : !> \param weights ...
497 : !> \param point_derivative ...
498 : !> \param atom_derivative ...
499 : !> \param atomic_weights ...
500 : !> \param calculate_derivatives ...
501 : ! **************************************************************************************************
502 225000 : SUBROUTINE hirshfeld_point_weights(context, point, weights, point_derivative, atom_derivative, &
503 225000 : atomic_weights, calculate_derivatives)
504 : TYPE(cdft_point_context_type), INTENT(INOUT) :: context
505 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: point
506 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: weights
507 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: point_derivative
508 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: atom_derivative
509 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: atomic_weights
510 : LOGICAL, INTENT(IN) :: calculate_derivatives
511 :
512 : INTEGER :: atom, iexp, igroup
513 : REAL(KIND=dp) :: exponential, numerator, sum_density
514 : REAL(KIND=dp), DIMENSION(3) :: dpoint_numerator, dsum_point
515 :
516 450000 : weights = 0.0_dp
517 225000 : IF (calculate_derivatives) THEN
518 50000 : point_derivative = 0.0_dp
519 100000 : atom_derivative = 0.0_dp
520 : END IF
521 375000 : IF (PRESENT(atomic_weights)) atomic_weights = 0.0_dp
522 675000 : context%cell_functions = 0.0_dp
523 305000 : IF (calculate_derivatives) context%density_atom_derivative = 0.0_dp
524 675000 : DO atom = 1, context%natom
525 1800000 : context%displacement(:, atom) = pbc(context%atom_coords(:, atom), point, context%cell)
526 1125000 : DO iexp = 1, context%numexp(atom)
527 : exponential = context%amplitude(iexp, atom)* &
528 : EXP(-context%alpha(iexp, atom)* &
529 1800000 : DOT_PRODUCT(context%displacement(:, atom), context%displacement(:, atom)))
530 450000 : context%cell_functions(atom) = context%cell_functions(atom) + exponential
531 900000 : IF (calculate_derivatives) THEN
532 : context%density_atom_derivative(:, atom) = &
533 : context%density_atom_derivative(:, atom) + &
534 80000 : 2.0_dp*context%alpha(iexp, atom)*context%displacement(:, atom)*exponential
535 : END IF
536 : END DO
537 : END DO
538 675000 : sum_density = SUM(context%cell_functions)
539 225000 : IF (sum_density <= context%eps) THEN
540 68602 : RETURN
541 : END IF
542 273218 : IF (PRESENT(atomic_weights)) atomic_weights = context%cell_functions/sum_density
543 215204 : IF (calculate_derivatives) dsum_point = -SUM(context%density_atom_derivative, DIM=2)
544 312796 : DO igroup = 1, context%ngroup
545 469194 : numerator = DOT_PRODUCT(context%coefficients(igroup, :), context%cell_functions)
546 156398 : weights(igroup) = numerator/sum_density
547 312796 : IF (calculate_derivatives) THEN
548 : dpoint_numerator(:) = &
549 137214 : MATMUL(context%density_atom_derivative, context%coefficients(igroup, :))
550 26136 : dpoint_numerator = -dpoint_numerator
551 : point_derivative(:, igroup) = &
552 26136 : (dpoint_numerator*sum_density - numerator*dsum_point)/sum_density**2
553 19602 : DO atom = 1, context%natom
554 : atom_derivative(:, atom, igroup) = &
555 : (context%coefficients(igroup, atom) - weights(igroup))* &
556 58806 : context%density_atom_derivative(:, atom)/sum_density
557 : END DO
558 : END IF
559 : END DO
560 : END SUBROUTINE hirshfeld_point_weights
561 :
562 0 : END MODULE qs_cdft_grid
|
