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