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