Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 30 : 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 30 : 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 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_nlcc
300 30 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_nlcc, hab, pab
301 30 : 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 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
306 : TYPE(pw_env_type), POINTER :: pw_env
307 30 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
308 30 : 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 30 : CALL timeset(routineN, handle)
313 :
314 30 : NULLIFY (pw_env, cores)
315 :
316 30 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
317 30 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
318 :
319 30 : 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 30 : force=force, virial=virial)
329 :
330 30 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
331 :
332 30 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
333 :
334 74 : DO ikind = 1, SIZE(atomic_kind_set)
335 :
336 44 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
337 44 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
338 :
339 44 : IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
340 : CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
341 44 : alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
342 :
343 44 : IF (.NOT. nlcc) CYCLE
344 :
345 202 : DO iexp_nlcc = 1, nexp_nlcc
346 :
347 42 : alpha = alpha_nlcc(iexp_nlcc)
348 42 : nc = nct_nlcc(iexp_nlcc)
349 :
350 42 : ni = ncoset(2*nc - 2)
351 :
352 42 : nthread = 1
353 42 : ithread = 0
354 :
355 168 : ALLOCATE (hab(ni, 1), pab(ni, 1))
356 270 : pab = 0._dp
357 :
358 42 : CALL reallocate(cores, 1, natom)
359 42 : npme = 0
360 172 : cores = 0
361 :
362 : ! prepare core function
363 100 : DO j = 1, nc
364 42 : SELECT CASE (j)
365 : CASE (1)
366 42 : 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 58 : CPABORT("")
410 : END SELECT
411 : END DO
412 42 : IF (dft_control%nspins == 2) pab = pab*0.5_dp
413 :
414 172 : DO iatom = 1, natom
415 130 : atom_a = atom_list(iatom)
416 130 : ra(:) = pbc(particle_set(atom_a)%r, cell)
417 172 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
418 : ! replicated realspace grid, split the atoms up between procs
419 130 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
420 65 : npme = npme + 1
421 65 : 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 107 : DO j = 1, npme
430 :
431 65 : iatom = cores(j)
432 65 : atom_a = atom_list(iatom)
433 65 : ra(:) = pbc(particle_set(atom_a)%r, cell)
434 274 : hab(:, 1) = 0.0_dp
435 65 : force_a(:) = 0.0_dp
436 65 : force_b(:) = 0.0_dp
437 65 : IF (use_virial) THEN
438 48 : my_virial_a = 0.0_dp
439 48 : my_virial_b = 0.0_dp
440 : END IF
441 65 : 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 65 : 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 65 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
456 :
457 : force(ikind)%gth_nlcc(:, iatom) = &
458 260 : force(ikind)%gth_nlcc(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
459 :
460 107 : 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 86 : DEALLOCATE (hab, pab)
467 :
468 : END DO
469 :
470 : END DO
471 :
472 30 : DEALLOCATE (cores)
473 :
474 30 : CALL timestop(handle)
475 :
476 30 : 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 7263 : 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 7263 : 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 7263 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
501 7263 : 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 7263 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
506 : TYPE(pw_env_type), POINTER :: pw_env
507 7263 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
508 7263 : 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 7263 : CALL timeset(routineN, handle)
513 7263 : NULLIFY (virial, force, atprop, dft_control)
514 :
515 7263 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
516 :
517 : !If gapw, check for gpw kinds
518 7263 : skip_fcore = .FALSE.
519 7263 : IF (dft_control%qs_control%gapw) THEN
520 496 : 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 6823 : NULLIFY (pw_env)
525 6823 : ALLOCATE (cores(1))
526 6823 : ALLOCATE (hab(1, 1))
527 6823 : ALLOCATE (pab(1, 1))
528 :
529 6823 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
530 6823 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
531 :
532 6823 : 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 6823 : atprop=atprop)
544 :
545 : ! atomic energy contributions
546 6823 : natom = SIZE(particle_set)
547 6823 : IF (ASSOCIATED(atprop)) THEN
548 6823 : CALL atprop_array_init(atprop%ateb, natom)
549 : END IF
550 :
551 6823 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
552 :
553 6823 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
554 :
555 18799 : DO ikind = 1, SIZE(atomic_kind_set)
556 :
557 11976 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
558 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
559 : alpha_core_charge=alpha_core_charge, &
560 11976 : ccore_charge=ccore_charge)
561 :
562 11976 : IF (dft_control%qs_control%gapw .AND. paw_atom) CYCLE
563 :
564 11914 : pab(1, 1) = -ccore_charge
565 11914 : IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
566 :
567 11878 : CALL reallocate(cores, 1, natom_of_kind)
568 11878 : npme = 0
569 36015 : cores = 0
570 :
571 36015 : DO iatom = 1, natom_of_kind
572 24137 : atom_a = atom_list(iatom)
573 24137 : ra(:) = pbc(particle_set(atom_a)%r, cell)
574 36015 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
575 : ! replicated realspace grid, split the atoms up between procs
576 23388 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
577 11694 : npme = npme + 1
578 11694 : cores(npme) = iatom
579 : END IF
580 : ELSE
581 749 : npme = npme + 1
582 749 : cores(npme) = iatom
583 : END IF
584 : END DO
585 :
586 43120 : DO j = 1, npme
587 :
588 12443 : iatom = cores(j)
589 12443 : atom_a = atom_list(iatom)
590 12443 : ra(:) = pbc(particle_set(atom_a)%r, cell)
591 12443 : hab(1, 1) = 0.0_dp
592 12443 : force_a(:) = 0.0_dp
593 12443 : force_b(:) = 0.0_dp
594 12443 : IF (use_virial) THEN
595 1558 : my_virial_a = 0.0_dp
596 1558 : my_virial_b = 0.0_dp
597 : END IF
598 :
599 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
600 : ra=ra, rb=ra, rp=ra, &
601 : zetp=alpha_core_charge, eps=eps_rho_rspace, &
602 : pab=pab, o1=0, o2=0, & ! without map_consistent
603 12443 : prefactor=1.0_dp, cutoff=1.0_dp)
604 :
605 : CALL integrate_pgf_product(0, alpha_core_charge, 0, &
606 : 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
607 : rs_v, hab, pab=pab, o1=0, o2=0, &
608 : radius=radius, &
609 : calculate_forces=.TRUE., force_a=force_a, &
610 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
611 12443 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
612 :
613 12443 : IF (ASSOCIATED(force)) THEN
614 49384 : force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
615 : END IF
616 12443 : IF (use_virial) THEN
617 20254 : virial%pv_ehartree = virial%pv_ehartree + my_virial_a
618 20254 : virial%pv_virial = virial%pv_virial + my_virial_a
619 : END IF
620 12443 : IF (ASSOCIATED(atprop)) THEN
621 12443 : atprop%ateb(atom_a) = atprop%ateb(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
622 : END IF
623 24419 : IF (PRESENT(atecc)) THEN
624 47 : atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
625 : END IF
626 :
627 : END DO
628 :
629 : END DO
630 :
631 6823 : DEALLOCATE (hab, pab, cores)
632 :
633 : END IF
634 :
635 7263 : CALL timestop(handle)
636 :
637 7263 : END SUBROUTINE integrate_v_core_rspace
638 :
639 : ! **************************************************************************************************
640 : !> \brief computes the overlap of a set of Gaussians with a potential on grid
641 : !> \param v_rspace ...
642 : !> \param qs_env ...
643 : !> \param alpha ...
644 : !> \param ccore ...
645 : !> \param atecc ...
646 : ! **************************************************************************************************
647 2 : SUBROUTINE integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
648 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
649 : TYPE(qs_environment_type), POINTER :: qs_env
650 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha, ccore
651 : REAL(KIND=dp), DIMENSION(:) :: atecc
652 :
653 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_gaussian_rspace'
654 :
655 : INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
656 : natom_of_kind, npme
657 2 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
658 : REAL(KIND=dp) :: alpha_core_charge, eps_rho_rspace, radius
659 : REAL(KIND=dp), DIMENSION(3) :: ra
660 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
661 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
662 : TYPE(cell_type), POINTER :: cell
663 : TYPE(dft_control_type), POINTER :: dft_control
664 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
665 : TYPE(pw_env_type), POINTER :: pw_env
666 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
667 : TYPE(realspace_grid_type), POINTER :: rs_v
668 :
669 2 : CALL timeset(routineN, handle)
670 :
671 2 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
672 :
673 : !If gapw, check for gpw kinds
674 2 : CPASSERT(.NOT. dft_control%qs_control%gapw)
675 :
676 2 : NULLIFY (pw_env)
677 2 : ALLOCATE (cores(1))
678 2 : ALLOCATE (hab(1, 1))
679 2 : ALLOCATE (pab(1, 1))
680 :
681 2 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
682 2 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
683 :
684 2 : CALL transfer_pw2rs(rs_v, v_rspace)
685 :
686 : CALL get_qs_env(qs_env=qs_env, &
687 : atomic_kind_set=atomic_kind_set, &
688 : qs_kind_set=qs_kind_set, &
689 : cell=cell, &
690 : dft_control=dft_control, &
691 : particle_set=particle_set, &
692 2 : pw_env=pw_env)
693 :
694 : ! atomic energy contributions
695 2 : natom = SIZE(particle_set)
696 2 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
697 :
698 6 : DO ikind = 1, SIZE(atomic_kind_set)
699 :
700 4 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
701 4 : pab(1, 1) = -ccore(ikind)
702 4 : alpha_core_charge = alpha(ikind)
703 4 : IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
704 :
705 4 : CALL reallocate(cores, 1, natom_of_kind)
706 4 : npme = 0
707 10 : cores = 0
708 :
709 10 : DO iatom = 1, natom_of_kind
710 6 : atom_a = atom_list(iatom)
711 6 : ra(:) = pbc(particle_set(atom_a)%r, cell)
712 10 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
713 : ! replicated realspace grid, split the atoms up between procs
714 6 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
715 3 : npme = npme + 1
716 3 : cores(npme) = iatom
717 : END IF
718 : ELSE
719 0 : npme = npme + 1
720 0 : cores(npme) = iatom
721 : END IF
722 : END DO
723 :
724 13 : DO j = 1, npme
725 :
726 3 : iatom = cores(j)
727 3 : atom_a = atom_list(iatom)
728 3 : ra(:) = pbc(particle_set(atom_a)%r, cell)
729 3 : hab(1, 1) = 0.0_dp
730 :
731 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
732 : ra=ra, rb=ra, rp=ra, &
733 : zetp=alpha_core_charge, eps=eps_rho_rspace, &
734 : pab=pab, o1=0, o2=0, & ! without map_consistent
735 3 : prefactor=1.0_dp, cutoff=1.0_dp)
736 :
737 : CALL integrate_pgf_product(0, alpha_core_charge, 0, &
738 : 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
739 : rs_v, hab, pab=pab, o1=0, o2=0, &
740 : radius=radius, calculate_forces=.FALSE., &
741 3 : use_subpatch=.TRUE., subpatch_pattern=0)
742 7 : atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
743 :
744 : END DO
745 :
746 : END DO
747 :
748 2 : DEALLOCATE (hab, pab, cores)
749 :
750 2 : CALL timestop(handle)
751 :
752 2 : END SUBROUTINE integrate_v_gaussian_rspace
753 : ! **************************************************************************************************
754 : !> \brief computes integrals of product of v_rspace times a one-center function
755 : !> required for LRIGPW
756 : !> \param v_rspace ...
757 : !> \param qs_env ...
758 : !> \param int_res ...
759 : !> \param calculate_forces ...
760 : !> \param basis_type ...
761 : !> \param atomlist ...
762 : !> \author Dorothea Golze
763 : ! **************************************************************************************************
764 862 : SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, &
765 862 : calculate_forces, basis_type, atomlist)
766 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
767 : TYPE(qs_environment_type), POINTER :: qs_env
768 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: int_res
769 : LOGICAL, INTENT(IN) :: calculate_forces
770 : CHARACTER(len=*), INTENT(IN) :: basis_type
771 : INTEGER, DIMENSION(:), OPTIONAL :: atomlist
772 :
773 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_one_center'
774 :
775 : INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, m1, maxco, &
776 : maxsgf_set, my_pos, na1, natom_of_kind, ncoa, nkind, nseta, offset, sgfa
777 862 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
778 862 : nsgf_seta
779 862 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
780 : LOGICAL :: use_virial
781 862 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it
782 : REAL(KIND=dp) :: eps_rho_rspace, radius
783 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
784 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
785 862 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
786 862 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab, rpgfa, sphi_a, work_f, work_i, &
787 862 : zeta
788 862 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
789 : TYPE(cell_type), POINTER :: cell
790 : TYPE(dft_control_type), POINTER :: dft_control
791 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
792 : TYPE(gto_basis_set_type), POINTER :: lri_basis_set
793 862 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
794 : TYPE(pw_env_type), POINTER :: pw_env
795 862 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
796 862 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
797 : TYPE(realspace_grid_type), POINTER :: rs_grid
798 : TYPE(virial_type), POINTER :: virial
799 :
800 862 : CALL timeset(routineN, handle)
801 :
802 862 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
803 862 : first_sgfa, gridlevel_info, hab, la_max, la_min, lri_basis_set, &
804 862 : npgfa, nsgf_seta, pab, particle_set, pw_env, rpgfa, &
805 862 : rs_grid, rs_v, virial, set_radius_a, sphi_a, work_f, &
806 862 : work_i, zeta)
807 :
808 862 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
809 :
810 862 : CALL pw_env_get(pw_env, rs_grids=rs_v)
811 4280 : DO i = 1, SIZE(rs_v)
812 4280 : CALL rs_grid_zero(rs_v(i))
813 : END DO
814 :
815 862 : gridlevel_info => pw_env%gridlevel_info
816 :
817 862 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
818 :
819 : CALL get_qs_env(qs_env=qs_env, &
820 : atomic_kind_set=atomic_kind_set, &
821 : qs_kind_set=qs_kind_set, &
822 : cell=cell, &
823 : dft_control=dft_control, &
824 : nkind=nkind, &
825 : particle_set=particle_set, &
826 : pw_env=pw_env, &
827 862 : virial=virial)
828 :
829 862 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
830 :
831 862 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
832 :
833 862 : offset = 0
834 862 : my_pos = v_rspace%pw_grid%para%group%mepos
835 862 : group_size = v_rspace%pw_grid%para%group%num_pe
836 :
837 2570 : DO ikind = 1, nkind
838 :
839 1708 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
840 1708 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
841 : CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
842 : first_sgf=first_sgfa, &
843 : lmax=la_max, &
844 : lmin=la_min, &
845 : maxco=maxco, &
846 : maxsgf_set=maxsgf_set, &
847 : npgf=npgfa, &
848 : nset=nseta, &
849 : nsgf_set=nsgf_seta, &
850 : pgf_radius=rpgfa, &
851 : set_radius=set_radius_a, &
852 : sphi=sphi_a, &
853 1708 : zet=zeta)
854 :
855 6832 : ALLOCATE (hab(maxco, 1), pab(maxco, 1))
856 45412 : hab = 0._dp
857 43704 : pab(:, 1) = 0._dp
858 :
859 4950 : DO iatom = 1, natom_of_kind
860 :
861 3242 : atom_a = atom_list(iatom)
862 3242 : IF (PRESENT(atomlist)) THEN
863 400 : IF (atomlist(atom_a) == 0) CYCLE
864 : END IF
865 3042 : ra(:) = pbc(particle_set(atom_a)%r, cell)
866 3042 : force_a(:) = 0._dp
867 3042 : force_b(:) = 0._dp
868 3042 : my_virial_a(:, :) = 0._dp
869 3042 : my_virial_b(:, :) = 0._dp
870 :
871 44304 : m1 = MAXVAL(npgfa(1:nseta))
872 9126 : ALLOCATE (map_it(m1))
873 :
874 44304 : DO iset = 1, nseta
875 : !
876 86728 : map_it = .FALSE.
877 86536 : DO ipgf = 1, npgfa(iset)
878 45274 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
879 45274 : rs_grid => rs_v(igrid_level)
880 86536 : map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
881 : END DO
882 41262 : offset = offset + 1
883 : !
884 66941 : IF (ANY(map_it(1:npgfa(iset)))) THEN
885 20631 : sgfa = first_sgfa(1, iset)
886 20631 : ncoa = npgfa(iset)*ncoset(la_max(iset))
887 435578 : hab(:, 1) = 0._dp
888 61893 : ALLOCATE (work_i(nsgf_seta(iset), 1))
889 202340 : work_i = 0.0_dp
890 :
891 : ! get fit coefficients for forces
892 20631 : IF (calculate_forces) THEN
893 895 : m1 = sgfa + nsgf_seta(iset) - 1
894 1790 : ALLOCATE (work_f(nsgf_seta(iset), 1))
895 6465 : 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 895 : SIZE(pab, 1))
899 895 : DEALLOCATE (work_f)
900 : END IF
901 :
902 43268 : DO ipgf = 1, npgfa(iset)
903 22637 : na1 = (ipgf - 1)*ncoset(la_max(iset))
904 22637 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
905 22637 : rs_grid => rs_v(igrid_level)
906 :
907 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
908 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
909 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
910 22637 : prefactor=1.0_dp, cutoff=1.0_dp)
911 :
912 43268 : IF (map_it(ipgf)) THEN
913 22637 : IF (.NOT. calculate_forces) THEN
914 : CALL integrate_pgf_product(la_max=la_max(iset), &
915 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
916 : lb_max=0, zetb=0.0_dp, lb_min=0, &
917 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
918 : rsgrid=rs_grid, &
919 : hab=hab, o1=na1, o2=0, radius=radius, &
920 21378 : calculate_forces=calculate_forces)
921 : ELSE
922 : CALL integrate_pgf_product(la_max=la_max(iset), &
923 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
924 : lb_max=0, zetb=0.0_dp, lb_min=0, &
925 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
926 : rsgrid=rs_grid, &
927 : hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
928 : calculate_forces=calculate_forces, &
929 : force_a=force_a, force_b=force_b, &
930 : use_virial=use_virial, &
931 1259 : my_virial_a=my_virial_a, my_virial_b=my_virial_b)
932 : END IF
933 : END IF
934 : END DO
935 : ! contract hab
936 : CALL dgemm("T", "N", nsgf_seta(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
937 20631 : SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work_i(1, 1), SIZE(work_i, 1))
938 :
939 : int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) = &
940 342787 : int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) + work_i(1:nsgf_seta(iset), 1)
941 20631 : DEALLOCATE (work_i)
942 : END IF
943 : END DO
944 : !
945 3042 : IF (calculate_forces) THEN
946 568 : int_res(ikind)%v_dfdr(iatom, :) = int_res(ikind)%v_dfdr(iatom, :) + force_a(:)
947 142 : IF (use_virial) THEN
948 156 : virial%pv_lrigpw = virial%pv_lrigpw + my_virial_a
949 156 : virial%pv_virial = virial%pv_virial + my_virial_a
950 : END IF
951 : END IF
952 :
953 4950 : DEALLOCATE (map_it)
954 :
955 : END DO
956 :
957 5986 : DEALLOCATE (hab, pab)
958 : END DO
959 :
960 862 : CALL timestop(handle)
961 :
962 1724 : END SUBROUTINE integrate_v_rspace_one_center
963 :
964 : ! **************************************************************************************************
965 : !> \brief computes integrals of product of v_rspace times the diagonal block basis functions
966 : !> required for LRIGPW with exact 1c terms
967 : !> \param v_rspace ...
968 : !> \param ksmat ...
969 : !> \param pmat ...
970 : !> \param qs_env ...
971 : !> \param calculate_forces ...
972 : !> \param basis_type ...
973 : !> \author JGH
974 : ! **************************************************************************************************
975 36 : SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
976 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
977 : TYPE(dbcsr_type), INTENT(INOUT) :: ksmat, pmat
978 : TYPE(qs_environment_type), POINTER :: qs_env
979 : LOGICAL, INTENT(IN) :: calculate_forces
980 : CHARACTER(len=*), INTENT(IN) :: basis_type
981 :
982 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_diagonal'
983 :
984 : INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
985 : m1, maxco, maxsgf_set, my_pos, na1, na2, natom_of_kind, nb1, nb2, ncoa, ncob, nkind, &
986 : nseta, nsgfa, offset, sgfa, sgfb
987 36 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
988 36 : nsgf_seta
989 36 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
990 : LOGICAL :: found, use_virial
991 36 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
992 : REAL(KIND=dp) :: eps_rho_rspace, radius, zetp
993 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
994 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
995 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
996 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, hab, hmat, p_block, pab, pblk, &
997 36 : rpgfa, sphi_a, work, zeta
998 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
999 : TYPE(cell_type), POINTER :: cell
1000 : TYPE(dft_control_type), POINTER :: dft_control
1001 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1002 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1003 : TYPE(mp_para_env_type), POINTER :: para_env
1004 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1005 : TYPE(pw_env_type), POINTER :: pw_env
1006 36 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1007 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1008 36 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1009 : TYPE(realspace_grid_type), POINTER :: rs_grid
1010 : TYPE(virial_type), POINTER :: virial
1011 :
1012 36 : CALL timeset(routineN, handle)
1013 :
1014 36 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1015 36 : CALL pw_env_get(pw_env, rs_grids=rs_v)
1016 36 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1017 :
1018 36 : gridlevel_info => pw_env%gridlevel_info
1019 :
1020 : CALL get_qs_env(qs_env=qs_env, &
1021 : atomic_kind_set=atomic_kind_set, &
1022 : qs_kind_set=qs_kind_set, &
1023 : cell=cell, &
1024 : dft_control=dft_control, &
1025 : nkind=nkind, &
1026 : particle_set=particle_set, &
1027 : force=force, &
1028 : virial=virial, &
1029 36 : para_env=para_env)
1030 :
1031 36 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1032 36 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1033 :
1034 36 : offset = 0
1035 36 : my_pos = v_rspace%pw_grid%para%group%mepos
1036 36 : group_size = v_rspace%pw_grid%para%group%num_pe
1037 :
1038 108 : DO ikind = 1, nkind
1039 :
1040 72 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
1041 72 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1042 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1043 : lmax=la_max, lmin=la_min, maxco=maxco, maxsgf_set=maxsgf_set, &
1044 : npgf=npgfa, nset=nseta, nsgf_set=nsgf_seta, nsgf=nsgfa, &
1045 : first_sgf=first_sgfa, pgf_radius=rpgfa, set_radius=set_radius_a, &
1046 72 : sphi=sphi_a, zet=zeta)
1047 :
1048 720 : ALLOCATE (hab(maxco, maxco), work(maxco, maxsgf_set), hmat(nsgfa, nsgfa))
1049 72 : IF (calculate_forces) ALLOCATE (pab(maxco, maxco), pblk(nsgfa, nsgfa))
1050 :
1051 288 : DO iatom = 1, natom_of_kind
1052 216 : atom_a = atom_list(iatom)
1053 216 : ra(:) = pbc(particle_set(atom_a)%r, cell)
1054 17640 : hmat = 0.0_dp
1055 216 : IF (calculate_forces) THEN
1056 0 : CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, BLOCK=p_block, found=found)
1057 0 : IF (found) THEN
1058 0 : pblk(1:nsgfa, 1:nsgfa) = p_block(1:nsgfa, 1:nsgfa)
1059 : ELSE
1060 0 : pblk = 0.0_dp
1061 : END IF
1062 0 : CALL para_env%sum(pblk)
1063 0 : force_a(:) = 0._dp
1064 0 : force_b(:) = 0._dp
1065 0 : IF (use_virial) THEN
1066 0 : my_virial_a = 0.0_dp
1067 0 : my_virial_b = 0.0_dp
1068 : END IF
1069 : END IF
1070 432 : m1 = MAXVAL(npgfa(1:nseta))
1071 864 : ALLOCATE (map_it2(m1, m1))
1072 432 : DO iset = 1, nseta
1073 216 : sgfa = first_sgfa(1, iset)
1074 216 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1075 648 : DO jset = 1, nseta
1076 216 : sgfb = first_sgfa(1, jset)
1077 216 : ncob = npgfa(jset)*ncoset(la_max(jset))
1078 : !
1079 12312 : map_it2 = .FALSE.
1080 1728 : DO ipgf = 1, npgfa(iset)
1081 12312 : DO jpgf = 1, npgfa(jset)
1082 10584 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1083 10584 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1084 10584 : rs_grid => rs_v(igrid_level)
1085 12096 : map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
1086 : END DO
1087 : END DO
1088 216 : offset = offset + 1
1089 : !
1090 6480 : IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
1091 237492 : hab(:, :) = 0._dp
1092 108 : IF (calculate_forces) THEN
1093 : CALL dgemm("N", "N", ncoa, nsgf_seta(jset), nsgf_seta(iset), &
1094 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1095 : pblk(sgfa, sgfb), SIZE(pblk, 1), &
1096 0 : 0.0_dp, work(1, 1), SIZE(work, 1))
1097 : CALL dgemm("N", "T", ncoa, ncob, nsgf_seta(jset), &
1098 : 1.0_dp, work(1, 1), SIZE(work, 1), &
1099 : sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1100 0 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
1101 : END IF
1102 :
1103 864 : DO ipgf = 1, npgfa(iset)
1104 756 : na1 = (ipgf - 1)*ncoset(la_max(iset))
1105 756 : na2 = ipgf*ncoset(la_max(iset))
1106 6156 : DO jpgf = 1, npgfa(jset)
1107 5292 : nb1 = (jpgf - 1)*ncoset(la_max(jset))
1108 5292 : nb2 = jpgf*ncoset(la_max(jset))
1109 5292 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1110 5292 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1111 5292 : rs_grid => rs_v(igrid_level)
1112 :
1113 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1114 : lb_min=la_min(jset), lb_max=la_max(jset), &
1115 : ra=ra, rb=ra, rp=ra, &
1116 : zetp=zetp, eps=eps_rho_rspace, &
1117 5292 : prefactor=1.0_dp, cutoff=1.0_dp)
1118 :
1119 6048 : IF (map_it2(ipgf, jpgf)) THEN
1120 5292 : IF (calculate_forces) THEN
1121 : CALL integrate_pgf_product( &
1122 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1123 : la_max(jset), zeta(jpgf, jset), la_min(jset), &
1124 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1125 : rsgrid=rs_v(igrid_level), &
1126 : hab=hab, pab=pab, o1=na1, o2=nb1, &
1127 : radius=radius, &
1128 : calculate_forces=.TRUE., &
1129 : force_a=force_a, force_b=force_b, &
1130 0 : use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1131 : ELSE
1132 : CALL integrate_pgf_product( &
1133 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1134 : la_max(jset), zeta(jpgf, jset), la_min(jset), &
1135 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1136 : rsgrid=rs_v(igrid_level), &
1137 : hab=hab, o1=na1, o2=nb1, &
1138 : radius=radius, &
1139 5292 : calculate_forces=.FALSE.)
1140 : END IF
1141 : END IF
1142 : END DO
1143 : END DO
1144 : ! contract hab
1145 : CALL dgemm("N", "N", ncoa, nsgf_seta(jset), ncob, &
1146 : 1.0_dp, hab(1, 1), SIZE(hab, 1), &
1147 : sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1148 108 : 0.0_dp, work(1, 1), SIZE(work, 1))
1149 : CALL dgemm("T", "N", nsgf_seta(iset), nsgf_seta(jset), ncoa, &
1150 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1151 : work(1, 1), SIZE(work, 1), &
1152 108 : 1.0_dp, hmat(sgfa, sgfb), SIZE(hmat, 1))
1153 : END IF
1154 : END DO
1155 : END DO
1156 216 : DEALLOCATE (map_it2)
1157 : ! update KS matrix
1158 35064 : CALL para_env%sum(hmat)
1159 216 : CALL dbcsr_get_block_p(matrix=ksmat, row=atom_a, col=atom_a, BLOCK=h_block, found=found)
1160 216 : IF (found) THEN
1161 17532 : h_block(1:nsgfa, 1:nsgfa) = h_block(1:nsgfa, 1:nsgfa) + hmat(1:nsgfa, 1:nsgfa)
1162 : END IF
1163 504 : IF (calculate_forces) THEN
1164 0 : force(ikind)%rho_elec(:, iatom) = force(ikind)%rho_elec(:, iatom) + 2.0_dp*force_a(:)
1165 0 : IF (use_virial) THEN
1166 : IF (use_virial .AND. calculate_forces) THEN
1167 0 : virial%pv_lrigpw = virial%pv_lrigpw + 2.0_dp*my_virial_a
1168 0 : virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a
1169 : END IF
1170 : END IF
1171 : END IF
1172 : END DO
1173 72 : DEALLOCATE (hab, work, hmat)
1174 252 : IF (calculate_forces) DEALLOCATE (pab, pblk)
1175 : END DO
1176 :
1177 36 : CALL timestop(handle)
1178 :
1179 72 : END SUBROUTINE integrate_v_rspace_diagonal
1180 :
1181 : ! **************************************************************************************************
1182 : !> \brief computes integrals of product of v_rspace times a basis function (vector function)
1183 : !> and possible forces
1184 : !> \param qs_env ...
1185 : !> \param v_rspace ...
1186 : !> \param f_coef ...
1187 : !> \param f_integral ...
1188 : !> \param calculate_forces ...
1189 : !> \param basis_type ...
1190 : !> \author JGH [8.2024]
1191 : ! **************************************************************************************************
1192 4 : SUBROUTINE integrate_function(qs_env, v_rspace, f_coef, f_integral, calculate_forces, basis_type)
1193 : TYPE(qs_environment_type), POINTER :: qs_env
1194 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
1195 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: f_coef
1196 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: f_integral
1197 : LOGICAL, INTENT(IN) :: calculate_forces
1198 : CHARACTER(len=*), INTENT(IN) :: basis_type
1199 :
1200 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_function'
1201 :
1202 : INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, &
1203 : maxsgf_set, my_pos, na1, natom, ncoa, nkind, nseta, offset, sgfa
1204 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1205 4 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
1206 4 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1207 : LOGICAL :: use_virial
1208 : REAL(KIND=dp) :: eps_rho_rspace, radius
1209 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
1210 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
1211 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab, sphi_a, work, zeta
1212 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1213 : TYPE(cell_type), POINTER :: cell
1214 : TYPE(dft_control_type), POINTER :: dft_control
1215 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1216 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1217 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1218 : TYPE(pw_env_type), POINTER :: pw_env
1219 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1220 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1221 4 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1222 : TYPE(realspace_grid_type), POINTER :: rs_grid
1223 : TYPE(virial_type), POINTER :: virial
1224 :
1225 4 : CALL timeset(routineN, handle)
1226 :
1227 4 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1228 4 : gridlevel_info => pw_env%gridlevel_info
1229 :
1230 4 : CALL pw_env_get(pw_env, rs_grids=rs_v)
1231 4 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1232 :
1233 : CALL get_qs_env(qs_env=qs_env, &
1234 : atomic_kind_set=atomic_kind_set, &
1235 : qs_kind_set=qs_kind_set, &
1236 : force=force, &
1237 : cell=cell, &
1238 : dft_control=dft_control, &
1239 : nkind=nkind, &
1240 : natom=natom, &
1241 : particle_set=particle_set, &
1242 : pw_env=pw_env, &
1243 4 : virial=virial)
1244 4 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1245 :
1246 4 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1247 4 : IF (use_virial) THEN
1248 0 : CPABORT("Virial NYA")
1249 : END IF
1250 :
1251 4 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1252 :
1253 : CALL get_qs_kind_set(qs_kind_set, &
1254 4 : maxco=maxco, maxsgf_set=maxsgf_set, basis_type=basis_type)
1255 20 : ALLOCATE (hab(maxco, 1), pab(maxco, 1), work(maxco, 1))
1256 :
1257 4 : offset = 0
1258 4 : my_pos = v_rspace%pw_grid%para%group%mepos
1259 4 : group_size = v_rspace%pw_grid%para%group%num_pe
1260 :
1261 20 : DO iatom = 1, natom
1262 16 : ikind = particle_set(iatom)%atomic_kind%kind_number
1263 16 : atom_a = atom_of_kind(iatom)
1264 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1265 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1266 : first_sgf=first_sgfa, &
1267 : lmax=la_max, &
1268 : lmin=la_min, &
1269 : npgf=npgfa, &
1270 : nset=nseta, &
1271 : nsgf_set=nsgfa, &
1272 : sphi=sphi_a, &
1273 16 : zet=zeta)
1274 16 : ra(:) = pbc(particle_set(iatom)%r, cell)
1275 :
1276 16 : force_a(:) = 0._dp
1277 16 : force_b(:) = 0._dp
1278 16 : my_virial_a(:, :) = 0._dp
1279 16 : my_virial_b(:, :) = 0._dp
1280 :
1281 32 : DO iset = 1, nseta
1282 :
1283 16 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1284 16 : sgfa = first_sgfa(1, iset)
1285 :
1286 144 : hab = 0._dp
1287 144 : pab = 0._dp
1288 :
1289 120 : DO i = 1, nsgfa(iset)
1290 120 : work(i, 1) = f_coef(offset + i)
1291 : END DO
1292 :
1293 : CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
1294 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1295 : work(1, 1), SIZE(work, 1), &
1296 16 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
1297 :
1298 120 : DO ipgf = 1, npgfa(iset)
1299 :
1300 104 : na1 = (ipgf - 1)*ncoset(la_max(iset))
1301 :
1302 104 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
1303 104 : rs_grid => rs_v(igrid_level)
1304 :
1305 120 : IF (map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)) THEN
1306 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1307 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
1308 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
1309 52 : prefactor=1.0_dp, cutoff=1.0_dp)
1310 :
1311 52 : IF (calculate_forces) THEN
1312 : CALL integrate_pgf_product(la_max=la_max(iset), &
1313 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1314 : lb_max=0, zetb=0.0_dp, lb_min=0, &
1315 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1316 : rsgrid=rs_grid, &
1317 : hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
1318 : calculate_forces=.TRUE., &
1319 : force_a=force_a, force_b=force_b, &
1320 : use_virial=use_virial, &
1321 52 : my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1322 : ELSE
1323 : CALL integrate_pgf_product(la_max=la_max(iset), &
1324 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1325 : lb_max=0, zetb=0.0_dp, lb_min=0, &
1326 : ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1327 : rsgrid=rs_grid, &
1328 : hab=hab, o1=na1, o2=0, radius=radius, &
1329 0 : calculate_forces=.FALSE.)
1330 : END IF
1331 :
1332 : END IF
1333 :
1334 : END DO
1335 : !
1336 : CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
1337 16 : SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work(1, 1), SIZE(work, 1))
1338 120 : DO i = 1, nsgfa(iset)
1339 120 : f_integral(offset + i) = work(i, 1)
1340 : END DO
1341 :
1342 32 : offset = offset + nsgfa(iset)
1343 :
1344 : END DO
1345 :
1346 36 : IF (calculate_forces) THEN
1347 64 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + force_a(:)
1348 16 : IF (use_virial) THEN
1349 0 : virial%pv_virial = virial%pv_virial + my_virial_a
1350 : END IF
1351 : END IF
1352 :
1353 : END DO
1354 :
1355 4 : DEALLOCATE (hab, pab, work)
1356 :
1357 4 : CALL timestop(handle)
1358 :
1359 12 : END SUBROUTINE integrate_function
1360 :
1361 : END MODULE qs_integrate_potential_single
|
