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 Build up the plane wave density by collocating the primitive Gaussian
10 : !> functions (pgf).
11 : !> \par History
12 : !> Joost VandeVondele (02.2002)
13 : !> 1) rewrote collocate_pgf for increased accuracy and speed
14 : !> 2) collocate_core hack for PGI compiler
15 : !> 3) added multiple grid feature
16 : !> 4) new way to go over the grid
17 : !> Joost VandeVondele (05.2002)
18 : !> 1) prelim. introduction of the real space grid type
19 : !> JGH [30.08.02] multigrid arrays independent from potential
20 : !> JGH [17.07.03] distributed real space code
21 : !> JGH [23.11.03] refactoring and new loop ordering
22 : !> JGH [04.12.03] OpneMP parallelization of main loops
23 : !> Joost VandeVondele (12.2003)
24 : !> 1) modified to compute tau
25 : !> Joost removed incremental build feature
26 : !> Joost introduced map consistent
27 : !> Rewrote grid integration/collocation routines, [Joost VandeVondele,03.2007]
28 : !> \author Matthias Krack (03.04.2001)
29 : ! **************************************************************************************************
30 : MODULE qs_integrate_potential_single
31 : USE ao_util, ONLY: exp_radius_very_extended
32 : USE atomic_kind_types, ONLY: atomic_kind_type,&
33 : get_atomic_kind,&
34 : get_atomic_kind_set
35 : USE atprop_types, ONLY: atprop_array_init,&
36 : atprop_type
37 : USE basis_set_types, ONLY: get_gto_basis_set,&
38 : gto_basis_set_type
39 : USE cell_types, ONLY: cell_type,&
40 : pbc
41 : USE cp_control_types, ONLY: dft_control_type
42 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
43 : dbcsr_type
44 : USE cp_log_handling, ONLY: cp_get_default_logger,&
45 : cp_logger_get_default_unit_nr,&
46 : cp_logger_type
47 : USE external_potential_types, ONLY: get_potential,&
48 : gth_potential_type
49 : USE gaussian_gridlevels, ONLY: gaussian_gridlevel,&
50 : gridlevel_info_type
51 : USE grid_api, ONLY: integrate_pgf_product
52 : USE kinds, ONLY: dp
53 : USE lri_environment_types, ONLY: lri_kind_type
54 : USE mathlib, ONLY: det_3x3
55 : USE memory_utilities, ONLY: reallocate
56 : USE message_passing, ONLY: mp_para_env_type
57 : USE orbital_pointers, ONLY: coset,&
58 : ncoset
59 : USE particle_types, ONLY: particle_type
60 : USE physcon, ONLY: pascal
61 : USE pw_env_types, ONLY: pw_env_get,&
62 : pw_env_type
63 : USE pw_types, ONLY: pw_r3d_rs_type
64 : USE qs_environment_types, ONLY: get_qs_env,&
65 : qs_environment_type
66 : USE qs_force_types, ONLY: qs_force_type
67 : USE qs_kind_types, ONLY: get_qs_kind,&
68 : get_qs_kind_set,&
69 : qs_kind_type
70 : USE qs_nlcc_methods, ONLY: get_nlcc_gaussians
71 : USE realspace_grid_types, ONLY: map_gaussian_here,&
72 : realspace_grid_type,&
73 : rs_grid_zero,&
74 : transfer_pw2rs
75 : USE rs_pw_interface, ONLY: potential_pw2rs
76 : USE virial_methods, ONLY: one_third_sum_diag
77 : USE virial_types, ONLY: virial_type
78 : #include "./base/base_uses.f90"
79 :
80 : IMPLICIT NONE
81 :
82 : PRIVATE
83 :
84 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
85 :
86 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_integrate_potential_single'
87 :
88 : ! *** Public subroutines ***
89 : ! *** Don't include this routines directly, use the interface to
90 : ! *** qs_integrate_potential
91 :
92 : PUBLIC :: integrate_v_rspace_one_center, &
93 : integrate_v_rspace_diagonal, &
94 : integrate_v_core_rspace, &
95 : integrate_v_gaussian_rspace, &
96 : integrate_function, &
97 : integrate_ppl_rspace, &
98 : integrate_rho_nlcc
99 :
100 : CONTAINS
101 :
102 : ! **************************************************************************************************
103 : !> \brief computes the forces/virial due to the local pseudopotential
104 : !> \param rho_rspace ...
105 : !> \param qs_env ...
106 : ! **************************************************************************************************
107 14 : SUBROUTINE integrate_ppl_rspace(rho_rspace, qs_env)
108 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_rspace
109 : TYPE(qs_environment_type), POINTER :: qs_env
110 :
111 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_ppl_rspace'
112 :
113 : INTEGER :: atom_a, handle, iatom, ikind, j, lppl, &
114 : n, natom_of_kind, ni, npme
115 14 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
116 : LOGICAL :: use_virial
117 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
118 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
119 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
120 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl
121 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
122 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
123 : TYPE(cell_type), POINTER :: cell
124 : TYPE(dft_control_type), POINTER :: dft_control
125 : TYPE(gth_potential_type), POINTER :: gth_potential
126 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
127 : TYPE(pw_env_type), POINTER :: pw_env
128 14 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
129 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
130 : TYPE(realspace_grid_type), POINTER :: rs_v
131 : TYPE(virial_type), POINTER :: virial
132 :
133 14 : CALL timeset(routineN, handle)
134 :
135 14 : NULLIFY (pw_env, cores)
136 :
137 14 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
138 14 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
139 :
140 14 : CALL transfer_pw2rs(rs_v, rho_rspace)
141 :
142 : CALL get_qs_env(qs_env=qs_env, &
143 : atomic_kind_set=atomic_kind_set, &
144 : qs_kind_set=qs_kind_set, &
145 : cell=cell, &
146 : dft_control=dft_control, &
147 : particle_set=particle_set, &
148 : pw_env=pw_env, &
149 14 : force=force, virial=virial)
150 :
151 14 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
152 :
153 14 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
154 :
155 34 : DO ikind = 1, SIZE(atomic_kind_set)
156 :
157 20 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
158 20 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
159 :
160 20 : IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
161 20 : CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
162 :
163 20 : IF (lppl <= 0) CYCLE
164 :
165 20 : ni = ncoset(2*lppl - 2)
166 80 : ALLOCATE (hab(ni, 1), pab(ni, 1))
167 240 : pab = 0._dp
168 :
169 20 : CALL reallocate(cores, 1, natom_of_kind)
170 20 : npme = 0
171 70 : cores = 0
172 :
173 : ! prepare core function
174 60 : DO j = 1, lppl
175 20 : SELECT CASE (j)
176 : CASE (1)
177 20 : pab(1, 1) = cexp_ppl(1)
178 : CASE (2)
179 20 : n = coset(2, 0, 0)
180 20 : pab(n, 1) = cexp_ppl(2)
181 20 : n = coset(0, 2, 0)
182 20 : pab(n, 1) = cexp_ppl(2)
183 20 : n = coset(0, 0, 2)
184 20 : pab(n, 1) = cexp_ppl(2)
185 : CASE (3)
186 0 : n = coset(4, 0, 0)
187 0 : pab(n, 1) = cexp_ppl(3)
188 0 : n = coset(0, 4, 0)
189 0 : pab(n, 1) = cexp_ppl(3)
190 0 : n = coset(0, 0, 4)
191 0 : pab(n, 1) = cexp_ppl(3)
192 0 : n = coset(2, 2, 0)
193 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
194 0 : n = coset(2, 0, 2)
195 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
196 0 : n = coset(0, 2, 2)
197 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
198 : CASE (4)
199 0 : n = coset(6, 0, 0)
200 0 : pab(n, 1) = cexp_ppl(4)
201 0 : n = coset(0, 6, 0)
202 0 : pab(n, 1) = cexp_ppl(4)
203 0 : n = coset(0, 0, 6)
204 0 : pab(n, 1) = cexp_ppl(4)
205 0 : n = coset(4, 2, 0)
206 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
207 0 : n = coset(4, 0, 2)
208 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
209 0 : n = coset(2, 4, 0)
210 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
211 0 : n = coset(2, 0, 4)
212 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
213 0 : n = coset(0, 4, 2)
214 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
215 0 : n = coset(0, 2, 4)
216 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
217 0 : n = coset(2, 2, 2)
218 0 : pab(n, 1) = 6._dp*cexp_ppl(4)
219 : CASE DEFAULT
220 : CALL cp_abort(__LOCATION__, &
221 : "Only 1, 2, 3, 4 are supported as the "// &
222 40 : "value of j in integrate_ppl_rspace")
223 : END SELECT
224 : END DO
225 :
226 70 : DO iatom = 1, natom_of_kind
227 50 : atom_a = atom_list(iatom)
228 50 : ra(:) = pbc(particle_set(atom_a)%r, cell)
229 70 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
230 : ! replicated realspace grid, split the atoms up between procs
231 50 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
232 25 : npme = npme + 1
233 25 : cores(npme) = iatom
234 : END IF
235 : ELSE
236 0 : npme = npme + 1
237 0 : cores(npme) = iatom
238 : END IF
239 : END DO
240 :
241 45 : DO j = 1, npme
242 :
243 25 : iatom = cores(j)
244 25 : atom_a = atom_list(iatom)
245 25 : ra(:) = pbc(particle_set(atom_a)%r, cell)
246 275 : hab(:, 1) = 0.0_dp
247 25 : force_a(:) = 0.0_dp
248 25 : force_b(:) = 0.0_dp
249 25 : IF (use_virial) THEN
250 0 : my_virial_a = 0.0_dp
251 0 : my_virial_b = 0.0_dp
252 : END IF
253 25 : ni = 2*lppl - 2
254 :
255 : radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
256 : ra=ra, rb=ra, rp=ra, &
257 : zetp=alpha, eps=eps_rho_rspace, &
258 : pab=pab, o1=0, o2=0, & ! without map_consistent
259 25 : prefactor=1.0_dp, cutoff=1.0_dp)
260 :
261 : CALL integrate_pgf_product(ni, alpha, 0, &
262 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
263 : rs_v, hab, pab=pab, o1=0, o2=0, &
264 : radius=radius, &
265 : calculate_forces=.TRUE., force_a=force_a, &
266 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
267 25 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
268 :
269 : force(ikind)%gth_ppl(:, iatom) = &
270 100 : force(ikind)%gth_ppl(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
271 :
272 45 : IF (use_virial) THEN
273 0 : virial%pv_ppl = virial%pv_ppl + my_virial_a*rho_rspace%pw_grid%dvol
274 0 : virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
275 0 : CPABORT("Virial not debuged for CORE_PPL")
276 : END IF
277 : END DO
278 :
279 74 : DEALLOCATE (hab, pab)
280 :
281 : END DO
282 :
283 14 : DEALLOCATE (cores)
284 :
285 14 : CALL timestop(handle)
286 :
287 14 : END SUBROUTINE integrate_ppl_rspace
288 :
289 : ! **************************************************************************************************
290 : !> \brief computes the forces/virial due to the nlcc pseudopotential
291 : !> \param rho_rspace ...
292 : !> \param qs_env ...
293 : !> \param debug_forces ...
294 : !> \param debug_stress ...
295 : !> \param iotag ...
296 : ! **************************************************************************************************
297 66 : SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env, debug_forces, debug_stress, iotag)
298 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_rspace
299 : TYPE(qs_environment_type), POINTER :: qs_env
300 : LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress
301 : CHARACTER(len=*), OPTIONAL :: iotag
302 :
303 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_rho_nlcc'
304 :
305 : CHARACTER(len=30) :: ftag, stag
306 : INTEGER :: atom_a, handle, iatom, iexp_nlcc, ikind, &
307 : iounit, j, n, natom, nc, nexp_nlcc, &
308 : ni, npme
309 66 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nct_nlcc
310 66 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
311 : LOGICAL :: debfor, debstr, use_virial
312 : REAL(KIND=dp) :: alpha, eps_rho_rspace, fconv, radius
313 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha_nlcc
314 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cval_nlcc
315 : REAL(KIND=dp), DIMENSION(3) :: fodeb, force_a, force_b, ra
316 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b, stdeb
317 66 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
318 66 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
319 : TYPE(cell_type), POINTER :: cell
320 : TYPE(cp_logger_type), POINTER :: logger
321 : TYPE(dft_control_type), POINTER :: dft_control
322 : TYPE(mp_para_env_type), POINTER :: para_env
323 66 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
324 : TYPE(pw_env_type), POINTER :: pw_env
325 66 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
326 66 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
327 : TYPE(realspace_grid_type), POINTER :: rs_v
328 : TYPE(virial_type), POINTER :: virial
329 :
330 66 : CALL timeset(routineN, handle)
331 :
332 66 : logger => cp_get_default_logger()
333 66 : IF (logger%para_env%is_source()) THEN
334 33 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
335 : ELSE
336 : iounit = -1
337 : END IF
338 :
339 66 : NULLIFY (pw_env, cores)
340 :
341 66 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
342 66 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
343 :
344 66 : CALL transfer_pw2rs(rs_v, rho_rspace)
345 :
346 : CALL get_qs_env(qs_env=qs_env, &
347 : atomic_kind_set=atomic_kind_set, &
348 : qs_kind_set=qs_kind_set, &
349 : cell=cell, &
350 : dft_control=dft_control, &
351 : particle_set=particle_set, &
352 : pw_env=pw_env, &
353 : para_env=para_env, &
354 66 : force=force, virial=virial)
355 :
356 66 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
357 : !debug
358 66 : debfor = .FALSE.
359 66 : IF (PRESENT(debug_forces)) debfor = debug_forces
360 66 : debstr = .FALSE.
361 66 : IF (PRESENT(debug_stress)) debstr = debug_stress
362 66 : debstr = debstr .AND. use_virial
363 66 : IF (debstr) THEN
364 0 : fconv = 1.0E-9_dp*pascal/cell%deth
365 : END IF
366 66 : IF (PRESENT(iotag)) THEN
367 18 : ftag = "DEBUG :: "//ADJUSTL(iotag)
368 18 : stag = "STRESS| INT "//ADJUSTL(iotag)
369 : ELSE
370 48 : ftag = "DEBUG :: dnlcc*V(r) "
371 48 : stag = "STRESS| INT dnlcc*V(r) "
372 : END IF
373 :
374 90 : IF (debfor) fodeb(1:3) = force(1)%gth_nlcc(1:3, 1)
375 66 : IF (debstr) stdeb = virial%pv_nlcc
376 :
377 66 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
378 :
379 188 : DO ikind = 1, SIZE(atomic_kind_set)
380 :
381 122 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
382 122 : CALL get_nlcc_gaussians(qs_kind_set(ikind), alpha_nlcc, nct_nlcc, cval_nlcc)
383 122 : nexp_nlcc = SIZE(alpha_nlcc)
384 :
385 676 : DO iexp_nlcc = 1, nexp_nlcc
386 :
387 366 : alpha = alpha_nlcc(iexp_nlcc)
388 366 : nc = nct_nlcc(iexp_nlcc)
389 :
390 366 : ni = ncoset(2*nc - 2)
391 :
392 1464 : ALLOCATE (hab(ni, 1), pab(ni, 1))
393 1242 : pab = 0._dp
394 :
395 366 : CALL reallocate(cores, 1, natom)
396 366 : npme = 0
397 820 : cores = 0
398 :
399 : ! prepare core function
400 748 : DO j = 1, nc
401 366 : SELECT CASE (j)
402 : CASE (1)
403 366 : pab(1, 1) = cval_nlcc(1, iexp_nlcc)
404 : CASE (2)
405 16 : n = coset(2, 0, 0)
406 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
407 16 : n = coset(0, 2, 0)
408 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
409 16 : n = coset(0, 0, 2)
410 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
411 : CASE (3)
412 0 : n = coset(4, 0, 0)
413 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
414 0 : n = coset(0, 4, 0)
415 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
416 0 : n = coset(0, 0, 4)
417 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
418 0 : n = coset(2, 2, 0)
419 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
420 0 : n = coset(2, 0, 2)
421 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
422 0 : n = coset(0, 2, 2)
423 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
424 : CASE (4)
425 0 : n = coset(6, 0, 0)
426 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
427 0 : n = coset(0, 6, 0)
428 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
429 0 : n = coset(0, 0, 6)
430 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
431 0 : n = coset(4, 2, 0)
432 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
433 0 : n = coset(4, 0, 2)
434 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
435 0 : n = coset(2, 4, 0)
436 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
437 0 : n = coset(2, 0, 4)
438 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
439 0 : n = coset(0, 4, 2)
440 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
441 0 : n = coset(0, 2, 4)
442 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
443 0 : n = coset(2, 2, 2)
444 0 : pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
445 : CASE DEFAULT
446 : CALL cp_abort(__LOCATION__, &
447 : "Only 1, 2, 3, 4 are supported as the "// &
448 382 : "value of j in integrate_rho_nlcc")
449 : END SELECT
450 : END DO
451 366 : IF (dft_control%nspins == 2) pab = pab*0.5_dp
452 :
453 820 : DO iatom = 1, natom
454 454 : atom_a = atom_list(iatom)
455 454 : ra(:) = pbc(particle_set(atom_a)%r, cell)
456 820 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
457 : ! replicated realspace grid, split the atoms up between procs
458 454 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
459 227 : npme = npme + 1
460 227 : cores(npme) = iatom
461 : END IF
462 : ELSE
463 0 : npme = npme + 1
464 0 : cores(npme) = iatom
465 : END IF
466 : END DO
467 :
468 593 : DO j = 1, npme
469 :
470 227 : iatom = cores(j)
471 227 : atom_a = atom_list(iatom)
472 227 : ra(:) = pbc(particle_set(atom_a)%r, cell)
473 598 : hab(:, 1) = 0.0_dp
474 227 : force_a(:) = 0.0_dp
475 227 : force_b(:) = 0.0_dp
476 227 : IF (use_virial) THEN
477 48 : my_virial_a = 0.0_dp
478 48 : my_virial_b = 0.0_dp
479 : END IF
480 227 : ni = 2*nc - 2
481 :
482 : radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
483 : ra=ra, rb=ra, rp=ra, &
484 : zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
485 : pab=pab, o1=0, o2=0, & ! without map_consistent
486 227 : prefactor=1.0_dp, cutoff=1.0_dp)
487 :
488 : CALL integrate_pgf_product(ni, 1/(2*alpha**2), 0, &
489 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
490 : rs_v, hab, pab=pab, o1=0, o2=0, &
491 : radius=radius, &
492 : calculate_forces=.TRUE., force_a=force_a, &
493 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
494 227 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
495 :
496 : force(ikind)%gth_nlcc(:, iatom) = &
497 908 : force(ikind)%gth_nlcc(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
498 :
499 593 : IF (use_virial) THEN
500 624 : virial%pv_nlcc = virial%pv_nlcc + my_virial_a*rho_rspace%pw_grid%dvol
501 624 : virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
502 : END IF
503 : END DO
504 :
505 488 : DEALLOCATE (hab, pab)
506 :
507 : END DO
508 :
509 : END DO
510 :
511 66 : IF (ASSOCIATED(cores)) DEALLOCATE (cores)
512 :
513 66 : IF (debfor) THEN
514 32 : fodeb(1:3) = force(1)%gth_nlcc(1:3, 1) - fodeb(1:3)
515 8 : CALL para_env%sum(fodeb)
516 8 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") TRIM(ftag), fodeb
517 : END IF
518 66 : IF (debstr) THEN
519 0 : stdeb = fconv*(virial%pv_nlcc - stdeb)
520 0 : CALL para_env%sum(stdeb)
521 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
522 0 : TRIM(stag), one_third_sum_diag(stdeb), det_3x3(stdeb)
523 : END IF
524 :
525 66 : CALL timestop(handle)
526 :
527 132 : END SUBROUTINE integrate_rho_nlcc
528 :
529 : ! **************************************************************************************************
530 : !> \brief computes the forces/virial due to the ionic cores with a potential on
531 : !> grid
532 : !> \param v_rspace ...
533 : !> \param qs_env ...
534 : !> \param atecc ...
535 : !> \param calpha optional Gaussian exponents per kind replacing the core charge exponents
536 : !> \param ccore optional Gaussian coefficients per kind replacing the core charge coefficients
537 : ! **************************************************************************************************
538 8189 : SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env, atecc, calpha, ccore)
539 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
540 : TYPE(qs_environment_type), POINTER :: qs_env
541 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atecc, calpha, ccore
542 :
543 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_core_rspace'
544 :
545 : INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
546 : natom_of_kind, npme
547 8189 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
548 : LOGICAL :: paw_atom, skip_fcore, use_virial
549 : REAL(KIND=dp) :: alpha_core_charge, ccore_charge, &
550 : eps_rho_rspace, radius
551 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
552 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
553 8189 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
554 8189 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
555 : TYPE(atprop_type), POINTER :: atprop
556 : TYPE(cell_type), POINTER :: cell
557 : TYPE(dft_control_type), POINTER :: dft_control
558 8189 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
559 : TYPE(pw_env_type), POINTER :: pw_env
560 8189 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
561 8189 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
562 : TYPE(realspace_grid_type), POINTER :: rs_v
563 : TYPE(virial_type), POINTER :: virial
564 :
565 8189 : CALL timeset(routineN, handle)
566 8189 : NULLIFY (virial, force, atprop, dft_control)
567 8189 : IF (PRESENT(calpha)) THEN
568 2 : CPASSERT(PRESENT(ccore))
569 : END IF
570 :
571 8189 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
572 :
573 : !If gapw, check for gpw kinds. This filtering applies to the core charges only;
574 : !user defined Gaussians (calpha/ccore) are always integrated
575 8189 : skip_fcore = .FALSE.
576 8189 : IF (dft_control%qs_control%gapw .AND. (.NOT. PRESENT(calpha))) THEN
577 880 : IF (.NOT. dft_control%qs_control%gapw_control%nopaw_as_gpw) skip_fcore = .TRUE.
578 : END IF
579 :
580 : IF (.NOT. skip_fcore) THEN
581 7397 : NULLIFY (pw_env)
582 7397 : ALLOCATE (cores(1))
583 7397 : ALLOCATE (hab(1, 1))
584 7397 : ALLOCATE (pab(1, 1))
585 :
586 7397 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
587 7397 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
588 :
589 7397 : CALL transfer_pw2rs(rs_v, v_rspace)
590 :
591 : CALL get_qs_env(qs_env=qs_env, &
592 : atomic_kind_set=atomic_kind_set, &
593 : qs_kind_set=qs_kind_set, &
594 : cell=cell, &
595 : dft_control=dft_control, &
596 : particle_set=particle_set, &
597 : pw_env=pw_env, &
598 : force=force, &
599 : virial=virial, &
600 7397 : atprop=atprop)
601 :
602 : ! atomic energy contributions
603 7397 : natom = SIZE(particle_set)
604 7397 : IF (ASSOCIATED(atprop)) THEN
605 7397 : CALL atprop_array_init(atprop%ateb, natom)
606 : END IF
607 :
608 7397 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
609 :
610 7397 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
611 :
612 20291 : DO ikind = 1, SIZE(atomic_kind_set)
613 :
614 12894 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
615 12894 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
616 :
617 12894 : IF (dft_control%qs_control%gapw .AND. paw_atom .AND. (.NOT. PRESENT(calpha))) CYCLE
618 :
619 12808 : IF (PRESENT(calpha)) THEN
620 : ! User defined Gaussians (e.g. the SCCS core fill) instead of the core charges
621 4 : alpha_core_charge = calpha(ikind)
622 4 : ccore_charge = ccore(ikind)
623 : ELSE
624 : CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha_core_charge, &
625 12804 : ccore_charge=ccore_charge)
626 : END IF
627 :
628 12808 : pab(1, 1) = -ccore_charge
629 12808 : IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
630 :
631 12770 : CALL reallocate(cores, 1, natom_of_kind)
632 12770 : npme = 0
633 39149 : cores = 0
634 :
635 39149 : DO iatom = 1, natom_of_kind
636 26379 : atom_a = atom_list(iatom)
637 26379 : ra(:) = pbc(particle_set(atom_a)%r, cell)
638 39149 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
639 : ! replicated realspace grid, split the atoms up between procs
640 25630 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
641 12815 : npme = npme + 1
642 12815 : cores(npme) = iatom
643 : END IF
644 : ELSE
645 749 : npme = npme + 1
646 749 : cores(npme) = iatom
647 : END IF
648 : END DO
649 :
650 46625 : DO j = 1, npme
651 :
652 13564 : iatom = cores(j)
653 13564 : atom_a = atom_list(iatom)
654 13564 : ra(:) = pbc(particle_set(atom_a)%r, cell)
655 13564 : hab(1, 1) = 0.0_dp
656 13564 : force_a(:) = 0.0_dp
657 13564 : force_b(:) = 0.0_dp
658 13564 : IF (use_virial) THEN
659 1749 : my_virial_a = 0.0_dp
660 1749 : my_virial_b = 0.0_dp
661 : END IF
662 :
663 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
664 : ra=ra, rb=ra, rp=ra, &
665 : zetp=alpha_core_charge, eps=eps_rho_rspace, &
666 : pab=pab, o1=0, o2=0, & ! without map_consistent
667 13564 : prefactor=1.0_dp, cutoff=1.0_dp)
668 :
669 : CALL integrate_pgf_product(0, alpha_core_charge, 0, &
670 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
671 : rs_v, hab, pab=pab, o1=0, o2=0, &
672 : radius=radius, &
673 : calculate_forces=.TRUE., force_a=force_a, &
674 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
675 13564 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
676 :
677 13564 : IF (ASSOCIATED(force)) THEN
678 53868 : force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
679 : END IF
680 13564 : IF (use_virial) THEN
681 22737 : virial%pv_ehartree = virial%pv_ehartree + my_virial_a
682 22737 : virial%pv_virial = virial%pv_virial + my_virial_a
683 : END IF
684 13564 : IF (ASSOCIATED(atprop) .AND. (.NOT. PRESENT(calpha))) THEN
685 13563 : atprop%ateb(atom_a) = atprop%ateb(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
686 : END IF
687 26458 : IF (PRESENT(atecc)) THEN
688 47 : atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
689 : END IF
690 :
691 : END DO
692 :
693 : END DO
694 :
695 7397 : DEALLOCATE (hab, pab, cores)
696 :
697 : END IF
698 :
699 8189 : CALL timestop(handle)
700 :
701 8189 : END SUBROUTINE integrate_v_core_rspace
702 :
703 : ! **************************************************************************************************
704 : !> \brief computes the overlap of a set of Gaussians with a potential on grid
705 : !> \param v_rspace ...
706 : !> \param qs_env ...
707 : !> \param alpha ...
708 : !> \param ccore ...
709 : !> \param atecc ...
710 : ! **************************************************************************************************
711 2 : SUBROUTINE integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
712 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
713 : TYPE(qs_environment_type), POINTER :: qs_env
714 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha, ccore
715 : REAL(KIND=dp), DIMENSION(:) :: atecc
716 :
717 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_gaussian_rspace'
718 :
719 : INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
720 : natom_of_kind, npme
721 2 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
722 : REAL(KIND=dp) :: alpha_core_charge, eps_rho_rspace, radius
723 : REAL(KIND=dp), DIMENSION(3) :: ra
724 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
725 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
726 : TYPE(cell_type), POINTER :: cell
727 : TYPE(dft_control_type), POINTER :: dft_control
728 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
729 : TYPE(pw_env_type), POINTER :: pw_env
730 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
731 : TYPE(realspace_grid_type), POINTER :: rs_v
732 :
733 2 : CALL timeset(routineN, handle)
734 :
735 2 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
736 :
737 : !If gapw, check for gpw kinds
738 2 : CPASSERT(.NOT. dft_control%qs_control%gapw)
739 :
740 2 : NULLIFY (pw_env)
741 2 : ALLOCATE (cores(1))
742 2 : ALLOCATE (hab(1, 1))
743 2 : ALLOCATE (pab(1, 1))
744 :
745 2 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
746 2 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
747 :
748 2 : CALL transfer_pw2rs(rs_v, v_rspace)
749 :
750 : CALL get_qs_env(qs_env=qs_env, &
751 : atomic_kind_set=atomic_kind_set, &
752 : qs_kind_set=qs_kind_set, &
753 : cell=cell, &
754 : dft_control=dft_control, &
755 : particle_set=particle_set, &
756 2 : pw_env=pw_env)
757 :
758 : ! atomic energy contributions
759 2 : natom = SIZE(particle_set)
760 2 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
761 :
762 6 : DO ikind = 1, SIZE(atomic_kind_set)
763 :
764 4 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
765 4 : pab(1, 1) = -ccore(ikind)
766 4 : alpha_core_charge = alpha(ikind)
767 4 : IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
768 :
769 4 : CALL reallocate(cores, 1, natom_of_kind)
770 4 : npme = 0
771 10 : cores = 0
772 :
773 10 : DO iatom = 1, natom_of_kind
774 6 : atom_a = atom_list(iatom)
775 6 : ra(:) = pbc(particle_set(atom_a)%r, cell)
776 10 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
777 : ! replicated realspace grid, split the atoms up between procs
778 6 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
779 3 : npme = npme + 1
780 3 : cores(npme) = iatom
781 : END IF
782 : ELSE
783 0 : npme = npme + 1
784 0 : cores(npme) = iatom
785 : END IF
786 : END DO
787 :
788 13 : DO j = 1, npme
789 :
790 3 : iatom = cores(j)
791 3 : atom_a = atom_list(iatom)
792 3 : ra(:) = pbc(particle_set(atom_a)%r, cell)
793 3 : hab(1, 1) = 0.0_dp
794 :
795 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
796 : ra=ra, rb=ra, rp=ra, &
797 : zetp=alpha_core_charge, eps=eps_rho_rspace, &
798 : pab=pab, o1=0, o2=0, & ! without map_consistent
799 3 : prefactor=1.0_dp, cutoff=1.0_dp)
800 :
801 : CALL integrate_pgf_product(0, alpha_core_charge, 0, &
802 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
803 : rs_v, hab, pab=pab, o1=0, o2=0, &
804 : radius=radius, calculate_forces=.FALSE., &
805 3 : use_subpatch=.TRUE., subpatch_pattern=0)
806 7 : atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
807 :
808 : END DO
809 :
810 : END DO
811 :
812 2 : DEALLOCATE (hab, pab, cores)
813 :
814 2 : CALL timestop(handle)
815 :
816 2 : END SUBROUTINE integrate_v_gaussian_rspace
817 : ! **************************************************************************************************
818 : !> \brief computes integrals of product of v_rspace times a one-center function
819 : !> required for LRIGPW
820 : !> \param v_rspace ...
821 : !> \param qs_env ...
822 : !> \param int_res ...
823 : !> \param calculate_forces ...
824 : !> \param basis_type ...
825 : !> \param atomlist ...
826 : !> \author Dorothea Golze
827 : ! **************************************************************************************************
828 1030 : SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, &
829 1030 : calculate_forces, basis_type, atomlist)
830 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
831 : TYPE(qs_environment_type), POINTER :: qs_env
832 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: int_res
833 : LOGICAL, INTENT(IN) :: calculate_forces
834 : CHARACTER(len=*), INTENT(IN) :: basis_type
835 : INTEGER, DIMENSION(:), OPTIONAL :: atomlist
836 :
837 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_one_center'
838 :
839 : INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, m1, &
840 : max_npgf, maxco, maxsgf_set, my_pos, na1, natom_of_kind, ncoa, nkind, nseta, offset, sgfa
841 1030 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
842 1030 : nsgf_seta
843 1030 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
844 : LOGICAL :: use_virial
845 1030 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it
846 : REAL(KIND=dp) :: eps_rho_rspace, radius
847 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
848 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
849 1030 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
850 1030 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab, rpgfa, sphi_a, work_f, work_i, &
851 1030 : zeta
852 1030 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
853 : TYPE(cell_type), POINTER :: cell
854 : TYPE(dft_control_type), POINTER :: dft_control
855 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
856 : TYPE(gto_basis_set_type), POINTER :: lri_basis_set
857 1030 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
858 : TYPE(pw_env_type), POINTER :: pw_env
859 1030 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
860 1030 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
861 : TYPE(realspace_grid_type), POINTER :: rs_grid
862 : TYPE(virial_type), POINTER :: virial
863 :
864 1030 : CALL timeset(routineN, handle)
865 :
866 1030 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
867 1030 : first_sgfa, gridlevel_info, hab, la_max, la_min, lri_basis_set, &
868 1030 : npgfa, nsgf_seta, pab, particle_set, pw_env, rpgfa, &
869 1030 : rs_grid, rs_v, virial, set_radius_a, sphi_a, work_f, &
870 1030 : work_i, zeta)
871 :
872 1030 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
873 :
874 1030 : CALL pw_env_get(pw_env, rs_grids=rs_v)
875 5114 : DO i = 1, SIZE(rs_v)
876 5114 : CALL rs_grid_zero(rs_v(i))
877 : END DO
878 :
879 1030 : gridlevel_info => pw_env%gridlevel_info
880 :
881 1030 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
882 :
883 : CALL get_qs_env(qs_env=qs_env, &
884 : atomic_kind_set=atomic_kind_set, &
885 : qs_kind_set=qs_kind_set, &
886 : cell=cell, &
887 : dft_control=dft_control, &
888 : nkind=nkind, &
889 : particle_set=particle_set, &
890 : pw_env=pw_env, &
891 1030 : virial=virial)
892 :
893 1030 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
894 :
895 1030 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
896 :
897 1030 : offset = 0
898 1030 : my_pos = v_rspace%pw_grid%para%group%mepos
899 1030 : group_size = v_rspace%pw_grid%para%group%num_pe
900 :
901 3040 : DO ikind = 1, nkind
902 :
903 2010 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
904 2010 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
905 : CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
906 : first_sgf=first_sgfa, &
907 : lmax=la_max, &
908 : lmin=la_min, &
909 : maxco=maxco, &
910 : maxsgf_set=maxsgf_set, &
911 : npgf=npgfa, &
912 : nset=nseta, &
913 : nsgf_set=nsgf_seta, &
914 : pgf_radius=rpgfa, &
915 : set_radius=set_radius_a, &
916 : sphi=sphi_a, &
917 2010 : zet=zeta)
918 :
919 8040 : ALLOCATE (hab(maxco, 1), pab(maxco, 1))
920 50828 : hab = 0._dp
921 48818 : pab(:, 1) = 0._dp
922 30372 : max_npgf = MAXVAL(npgfa(1:nseta))
923 6030 : ALLOCATE (map_it(max_npgf))
924 6030 : ALLOCATE (work_i(maxsgf_set, 1))
925 2106 : IF (calculate_forces) ALLOCATE (work_f(maxsgf_set, 1))
926 :
927 5998 : DO iatom = 1, natom_of_kind
928 :
929 3988 : atom_a = atom_list(iatom)
930 3988 : IF (PRESENT(atomlist)) THEN
931 800 : IF (atomlist(atom_a) == 0) CYCLE
932 : END IF
933 3588 : ra(:) = pbc(particle_set(atom_a)%r, cell)
934 3588 : force_a(:) = 0._dp
935 3588 : force_b(:) = 0._dp
936 3588 : my_virial_a(:, :) = 0._dp
937 3588 : my_virial_b(:, :) = 0._dp
938 :
939 54974 : DO iset = 1, nseta
940 : !
941 51386 : map_it = .FALSE.
942 107264 : DO ipgf = 1, npgfa(iset)
943 55878 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
944 55878 : rs_grid => rs_v(igrid_level)
945 107264 : map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
946 : END DO
947 51386 : offset = offset + 1
948 : !
949 82913 : IF (ANY(map_it(1:npgfa(iset)))) THEN
950 25693 : sgfa = first_sgfa(1, iset)
951 25693 : ncoa = npgfa(iset)*ncoset(la_max(iset))
952 504220 : hab(:, 1) = 0._dp
953 212423 : work_i(1:nsgf_seta(iset), 1) = 0.0_dp
954 :
955 : ! get fit coefficients for forces
956 25693 : IF (calculate_forces) THEN
957 1379 : m1 = sgfa + nsgf_seta(iset) - 1
958 9649 : work_f(1:nsgf_seta(iset), 1) = int_res(ikind)%acoef(iatom, sgfa:m1)
959 : CALL dgemm("N", "N", ncoa, 1, nsgf_seta(iset), 1.0_dp, sphi_a(1, sgfa), &
960 : SIZE(sphi_a, 1), work_f(1, 1), SIZE(work_f, 1), 0.0_dp, pab(1, 1), &
961 1379 : SIZE(pab, 1))
962 : END IF
963 :
964 53632 : DO ipgf = 1, npgfa(iset)
965 27939 : na1 = (ipgf - 1)*ncoset(la_max(iset))
966 27939 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
967 27939 : rs_grid => rs_v(igrid_level)
968 :
969 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
970 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
971 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
972 27939 : prefactor=1.0_dp, cutoff=1.0_dp)
973 :
974 53632 : IF (map_it(ipgf)) THEN
975 27939 : IF (.NOT. calculate_forces) THEN
976 : CALL integrate_pgf_product(la_max=la_max(iset), &
977 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
978 : lb_max=0, zetb=0.0_dp, lb_min=0, &
979 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
980 : rsgrid=rs_grid, &
981 : hab=hab, o1=na1, o2=0, radius=radius, &
982 26186 : calculate_forces=calculate_forces)
983 : ELSE
984 : CALL integrate_pgf_product(la_max=la_max(iset), &
985 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
986 : lb_max=0, zetb=0.0_dp, lb_min=0, &
987 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
988 : rsgrid=rs_grid, &
989 : hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
990 : calculate_forces=calculate_forces, &
991 : force_a=force_a, force_b=force_b, &
992 : use_virial=use_virial, &
993 1753 : my_virial_a=my_virial_a, my_virial_b=my_virial_b)
994 : END IF
995 : END IF
996 : END DO
997 : ! contract hab
998 : CALL dgemm("T", "N", nsgf_seta(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
999 25693 : SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work_i(1, 1), SIZE(work_i, 1))
1000 :
1001 : int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) = &
1002 399153 : int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) + work_i(1:nsgf_seta(iset), 1)
1003 : END IF
1004 : END DO
1005 : !
1006 5598 : IF (calculate_forces) THEN
1007 808 : int_res(ikind)%v_dfdr(iatom, :) = int_res(ikind)%v_dfdr(iatom, :) + force_a(:)
1008 202 : IF (use_virial) THEN
1009 676 : virial%pv_lrigpw = virial%pv_lrigpw + my_virial_a
1010 676 : virial%pv_virial = virial%pv_virial + my_virial_a
1011 : END IF
1012 : END IF
1013 :
1014 : END DO
1015 :
1016 2010 : IF (calculate_forces) DEALLOCATE (work_f)
1017 2010 : DEALLOCATE (work_i, map_it)
1018 7060 : DEALLOCATE (hab, pab)
1019 : END DO
1020 :
1021 1030 : CALL timestop(handle)
1022 :
1023 2060 : END SUBROUTINE integrate_v_rspace_one_center
1024 :
1025 : ! **************************************************************************************************
1026 : !> \brief computes integrals of product of v_rspace times the diagonal block basis functions
1027 : !> required for LRIGPW with exact 1c terms
1028 : !> \param v_rspace ...
1029 : !> \param ksmat ...
1030 : !> \param pmat ...
1031 : !> \param qs_env ...
1032 : !> \param calculate_forces ...
1033 : !> \param basis_type ...
1034 : !> \author JGH
1035 : ! **************************************************************************************************
1036 36 : SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
1037 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
1038 : TYPE(dbcsr_type), INTENT(INOUT) :: ksmat, pmat
1039 : TYPE(qs_environment_type), POINTER :: qs_env
1040 : LOGICAL, INTENT(IN) :: calculate_forces
1041 : CHARACTER(len=*), INTENT(IN) :: basis_type
1042 :
1043 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_diagonal'
1044 :
1045 : INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
1046 : m1, maxco, maxsgf_set, my_pos, na1, na2, natom_of_kind, nb1, nb2, ncoa, ncob, nkind, &
1047 : nseta, nsgfa, offset, sgfa, sgfb
1048 36 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
1049 36 : nsgf_seta
1050 36 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1051 : LOGICAL :: found, use_virial
1052 36 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
1053 : REAL(KIND=dp) :: eps_rho_rspace, radius, zetp
1054 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
1055 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
1056 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
1057 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, hab, hmat, p_block, pab, pblk, &
1058 36 : rpgfa, sphi_a, work, zeta
1059 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1060 : TYPE(cell_type), POINTER :: cell
1061 : TYPE(dft_control_type), POINTER :: dft_control
1062 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1063 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1064 : TYPE(mp_para_env_type), POINTER :: para_env
1065 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1066 : TYPE(pw_env_type), POINTER :: pw_env
1067 36 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1068 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1069 36 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1070 : TYPE(realspace_grid_type), POINTER :: rs_grid
1071 : TYPE(virial_type), POINTER :: virial
1072 :
1073 36 : CALL timeset(routineN, handle)
1074 :
1075 36 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1076 36 : CALL pw_env_get(pw_env, rs_grids=rs_v)
1077 36 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1078 :
1079 36 : gridlevel_info => pw_env%gridlevel_info
1080 :
1081 : CALL get_qs_env(qs_env=qs_env, &
1082 : atomic_kind_set=atomic_kind_set, &
1083 : qs_kind_set=qs_kind_set, &
1084 : cell=cell, &
1085 : dft_control=dft_control, &
1086 : nkind=nkind, &
1087 : particle_set=particle_set, &
1088 : force=force, &
1089 : virial=virial, &
1090 36 : para_env=para_env)
1091 :
1092 36 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1093 36 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1094 :
1095 36 : offset = 0
1096 36 : my_pos = v_rspace%pw_grid%para%group%mepos
1097 36 : group_size = v_rspace%pw_grid%para%group%num_pe
1098 :
1099 108 : DO ikind = 1, nkind
1100 :
1101 72 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
1102 72 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1103 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1104 : lmax=la_max, lmin=la_min, maxco=maxco, maxsgf_set=maxsgf_set, &
1105 : npgf=npgfa, nset=nseta, nsgf_set=nsgf_seta, nsgf=nsgfa, &
1106 : first_sgf=first_sgfa, pgf_radius=rpgfa, set_radius=set_radius_a, &
1107 72 : sphi=sphi_a, zet=zeta)
1108 :
1109 720 : ALLOCATE (hab(maxco, maxco), work(maxco, maxsgf_set), hmat(nsgfa, nsgfa))
1110 72 : IF (calculate_forces) ALLOCATE (pab(maxco, maxco), pblk(nsgfa, nsgfa))
1111 :
1112 288 : DO iatom = 1, natom_of_kind
1113 216 : atom_a = atom_list(iatom)
1114 216 : ra(:) = pbc(particle_set(atom_a)%r, cell)
1115 17640 : hmat = 0.0_dp
1116 216 : IF (calculate_forces) THEN
1117 0 : CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, BLOCK=p_block, found=found)
1118 0 : IF (found) THEN
1119 0 : pblk(1:nsgfa, 1:nsgfa) = p_block(1:nsgfa, 1:nsgfa)
1120 : ELSE
1121 0 : pblk = 0.0_dp
1122 : END IF
1123 0 : CALL para_env%sum(pblk)
1124 0 : force_a(:) = 0._dp
1125 0 : force_b(:) = 0._dp
1126 0 : IF (use_virial) THEN
1127 0 : my_virial_a = 0.0_dp
1128 0 : my_virial_b = 0.0_dp
1129 : END IF
1130 : END IF
1131 432 : m1 = MAXVAL(npgfa(1:nseta))
1132 864 : ALLOCATE (map_it2(m1, m1))
1133 432 : DO iset = 1, nseta
1134 216 : sgfa = first_sgfa(1, iset)
1135 216 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1136 648 : DO jset = 1, nseta
1137 216 : sgfb = first_sgfa(1, jset)
1138 216 : ncob = npgfa(jset)*ncoset(la_max(jset))
1139 : !
1140 216 : map_it2 = .FALSE.
1141 1728 : DO ipgf = 1, npgfa(iset)
1142 12312 : DO jpgf = 1, npgfa(jset)
1143 10584 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1144 10584 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1145 10584 : rs_grid => rs_v(igrid_level)
1146 12096 : map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
1147 : END DO
1148 : END DO
1149 216 : offset = offset + 1
1150 : !
1151 6480 : IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
1152 237492 : hab(:, :) = 0._dp
1153 108 : IF (calculate_forces) THEN
1154 : CALL dgemm("N", "N", ncoa, nsgf_seta(jset), nsgf_seta(iset), &
1155 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1156 : pblk(sgfa, sgfb), SIZE(pblk, 1), &
1157 0 : 0.0_dp, work(1, 1), SIZE(work, 1))
1158 : CALL dgemm("N", "T", ncoa, ncob, nsgf_seta(jset), &
1159 : 1.0_dp, work(1, 1), SIZE(work, 1), &
1160 : sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1161 0 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
1162 : END IF
1163 :
1164 864 : DO ipgf = 1, npgfa(iset)
1165 756 : na1 = (ipgf - 1)*ncoset(la_max(iset))
1166 756 : na2 = ipgf*ncoset(la_max(iset))
1167 6156 : DO jpgf = 1, npgfa(jset)
1168 5292 : nb1 = (jpgf - 1)*ncoset(la_max(jset))
1169 5292 : nb2 = jpgf*ncoset(la_max(jset))
1170 5292 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1171 5292 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1172 5292 : rs_grid => rs_v(igrid_level)
1173 :
1174 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1175 : lb_min=la_min(jset), lb_max=la_max(jset), &
1176 : ra=ra, rb=ra, rp=ra, &
1177 : zetp=zetp, eps=eps_rho_rspace, &
1178 5292 : prefactor=1.0_dp, cutoff=1.0_dp)
1179 :
1180 6048 : IF (map_it2(ipgf, jpgf)) THEN
1181 5292 : IF (calculate_forces) THEN
1182 : CALL integrate_pgf_product( &
1183 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1184 : la_max(jset), zeta(jpgf, jset), la_min(jset), &
1185 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1186 : rsgrid=rs_v(igrid_level), &
1187 : hab=hab, pab=pab, o1=na1, o2=nb1, &
1188 : radius=radius, &
1189 : calculate_forces=.TRUE., &
1190 : force_a=force_a, force_b=force_b, &
1191 0 : use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1192 : ELSE
1193 : CALL integrate_pgf_product( &
1194 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1195 : la_max(jset), zeta(jpgf, jset), la_min(jset), &
1196 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1197 : rsgrid=rs_v(igrid_level), &
1198 : hab=hab, o1=na1, o2=nb1, &
1199 : radius=radius, &
1200 5292 : calculate_forces=.FALSE.)
1201 : END IF
1202 : END IF
1203 : END DO
1204 : END DO
1205 : ! contract hab
1206 : CALL dgemm("N", "N", ncoa, nsgf_seta(jset), ncob, &
1207 : 1.0_dp, hab(1, 1), SIZE(hab, 1), &
1208 : sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1209 108 : 0.0_dp, work(1, 1), SIZE(work, 1))
1210 : CALL dgemm("T", "N", nsgf_seta(iset), nsgf_seta(jset), ncoa, &
1211 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1212 : work(1, 1), SIZE(work, 1), &
1213 108 : 1.0_dp, hmat(sgfa, sgfb), SIZE(hmat, 1))
1214 : END IF
1215 : END DO
1216 : END DO
1217 216 : DEALLOCATE (map_it2)
1218 : ! update KS matrix
1219 35064 : CALL para_env%sum(hmat)
1220 216 : CALL dbcsr_get_block_p(matrix=ksmat, row=atom_a, col=atom_a, BLOCK=h_block, found=found)
1221 216 : IF (found) THEN
1222 17532 : h_block(1:nsgfa, 1:nsgfa) = h_block(1:nsgfa, 1:nsgfa) + hmat(1:nsgfa, 1:nsgfa)
1223 : END IF
1224 504 : IF (calculate_forces) THEN
1225 0 : force(ikind)%rho_elec(:, iatom) = force(ikind)%rho_elec(:, iatom) + 2.0_dp*force_a(:)
1226 0 : IF (use_virial) THEN
1227 : IF (use_virial .AND. calculate_forces) THEN
1228 0 : virial%pv_lrigpw = virial%pv_lrigpw + 2.0_dp*my_virial_a
1229 0 : virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a
1230 : END IF
1231 : END IF
1232 : END IF
1233 : END DO
1234 72 : DEALLOCATE (hab, work, hmat)
1235 252 : IF (calculate_forces) DEALLOCATE (pab, pblk)
1236 : END DO
1237 :
1238 36 : CALL timestop(handle)
1239 :
1240 72 : END SUBROUTINE integrate_v_rspace_diagonal
1241 :
1242 : ! **************************************************************************************************
1243 : !> \brief computes integrals of product of v_rspace times a basis function (vector function)
1244 : !> and possible forces
1245 : !> \param qs_env ...
1246 : !> \param v_rspace ...
1247 : !> \param f_coef ...
1248 : !> \param f_integral ...
1249 : !> \param calculate_forces ...
1250 : !> \param basis_type ...
1251 : !> \author JGH [8.2024]
1252 : ! **************************************************************************************************
1253 8 : SUBROUTINE integrate_function(qs_env, v_rspace, f_coef, f_integral, calculate_forces, basis_type)
1254 : TYPE(qs_environment_type), POINTER :: qs_env
1255 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
1256 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: f_coef
1257 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: f_integral
1258 : LOGICAL, INTENT(IN) :: calculate_forces
1259 : CHARACTER(len=*), INTENT(IN) :: basis_type
1260 :
1261 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_function'
1262 :
1263 : INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, &
1264 : maxsgf_set, my_pos, na1, natom, ncoa, nkind, nseta, offset, sgfa
1265 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1266 8 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
1267 8 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1268 : LOGICAL :: use_virial
1269 : REAL(KIND=dp) :: eps_rho_rspace, radius
1270 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
1271 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
1272 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab, sphi_a, work, zeta
1273 8 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1274 : TYPE(cell_type), POINTER :: cell
1275 : TYPE(dft_control_type), POINTER :: dft_control
1276 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1277 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1278 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1279 : TYPE(pw_env_type), POINTER :: pw_env
1280 8 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1281 8 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1282 8 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1283 : TYPE(realspace_grid_type), POINTER :: rs_grid
1284 : TYPE(virial_type), POINTER :: virial
1285 :
1286 8 : CALL timeset(routineN, handle)
1287 :
1288 8 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1289 8 : gridlevel_info => pw_env%gridlevel_info
1290 :
1291 8 : CALL pw_env_get(pw_env, rs_grids=rs_v)
1292 8 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1293 :
1294 : CALL get_qs_env(qs_env=qs_env, &
1295 : atomic_kind_set=atomic_kind_set, &
1296 : qs_kind_set=qs_kind_set, &
1297 : force=force, &
1298 : cell=cell, &
1299 : dft_control=dft_control, &
1300 : nkind=nkind, &
1301 : natom=natom, &
1302 : particle_set=particle_set, &
1303 : pw_env=pw_env, &
1304 8 : virial=virial)
1305 8 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1306 :
1307 8 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1308 8 : IF (use_virial) THEN
1309 0 : CPABORT("Virial NYA")
1310 : END IF
1311 :
1312 8 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1313 :
1314 : CALL get_qs_kind_set(qs_kind_set, &
1315 8 : maxco=maxco, maxsgf_set=maxsgf_set, basis_type=basis_type)
1316 40 : ALLOCATE (hab(maxco, 1), pab(maxco, 1), work(maxco, 1))
1317 :
1318 8 : offset = 0
1319 8 : my_pos = v_rspace%pw_grid%para%group%mepos
1320 8 : group_size = v_rspace%pw_grid%para%group%num_pe
1321 :
1322 40 : DO iatom = 1, natom
1323 32 : ikind = particle_set(iatom)%atomic_kind%kind_number
1324 32 : atom_a = atom_of_kind(iatom)
1325 32 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1326 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1327 : first_sgf=first_sgfa, &
1328 : lmax=la_max, &
1329 : lmin=la_min, &
1330 : npgf=npgfa, &
1331 : nset=nseta, &
1332 : nsgf_set=nsgfa, &
1333 : sphi=sphi_a, &
1334 32 : zet=zeta)
1335 32 : ra(:) = pbc(particle_set(iatom)%r, cell)
1336 :
1337 32 : force_a(:) = 0._dp
1338 32 : force_b(:) = 0._dp
1339 32 : my_virial_a(:, :) = 0._dp
1340 32 : my_virial_b(:, :) = 0._dp
1341 :
1342 64 : DO iset = 1, nseta
1343 :
1344 32 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1345 32 : sgfa = first_sgfa(1, iset)
1346 :
1347 288 : hab = 0._dp
1348 288 : pab = 0._dp
1349 :
1350 240 : DO i = 1, nsgfa(iset)
1351 240 : work(i, 1) = f_coef(offset + i)
1352 : END DO
1353 :
1354 : CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
1355 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1356 : work(1, 1), SIZE(work, 1), &
1357 32 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
1358 :
1359 240 : DO ipgf = 1, npgfa(iset)
1360 :
1361 208 : na1 = (ipgf - 1)*ncoset(la_max(iset))
1362 :
1363 208 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
1364 208 : rs_grid => rs_v(igrid_level)
1365 :
1366 240 : IF (map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)) THEN
1367 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1368 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
1369 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
1370 104 : prefactor=1.0_dp, cutoff=1.0_dp)
1371 :
1372 104 : IF (calculate_forces) THEN
1373 : CALL integrate_pgf_product(la_max=la_max(iset), &
1374 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1375 : lb_max=0, zetb=0.0_dp, lb_min=0, &
1376 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1377 : rsgrid=rs_grid, &
1378 : hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
1379 : calculate_forces=.TRUE., &
1380 : force_a=force_a, force_b=force_b, &
1381 : use_virial=use_virial, &
1382 104 : my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1383 : ELSE
1384 : CALL integrate_pgf_product(la_max=la_max(iset), &
1385 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1386 : lb_max=0, zetb=0.0_dp, lb_min=0, &
1387 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1388 : rsgrid=rs_grid, &
1389 : hab=hab, o1=na1, o2=0, radius=radius, &
1390 0 : calculate_forces=.FALSE.)
1391 : END IF
1392 :
1393 : END IF
1394 :
1395 : END DO
1396 : !
1397 : CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
1398 32 : SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work(1, 1), SIZE(work, 1))
1399 240 : DO i = 1, nsgfa(iset)
1400 240 : f_integral(offset + i) = work(i, 1)
1401 : END DO
1402 :
1403 64 : offset = offset + nsgfa(iset)
1404 :
1405 : END DO
1406 :
1407 72 : IF (calculate_forces) THEN
1408 128 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + force_a(:)
1409 32 : IF (use_virial) THEN
1410 0 : virial%pv_virial = virial%pv_virial + my_virial_a
1411 : END IF
1412 : END IF
1413 :
1414 : END DO
1415 :
1416 8 : DEALLOCATE (hab, pab, work)
1417 :
1418 8 : CALL timestop(handle)
1419 :
1420 24 : END SUBROUTINE integrate_function
1421 :
1422 : END MODULE qs_integrate_potential_single
|
