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