(git:33f85d8)
Loading...
Searching...
No Matches
qs_integrate_potential_single.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
30MODULE 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
93CONTAINS
94
95! **************************************************************************************************
96!> \brief computes the forces/virial due to the local pseudopotential
97!> \param rho_rspace ...
98!> \param qs_env ...
99! **************************************************************************************************
100 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: cexp_ppl
114 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab
115 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
120 TYPE(pw_env_type), POINTER :: pw_env
121 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
122 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 CALL timeset(routinen, handle)
127
128 NULLIFY (pw_env, cores)
129
130 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
131 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
132
133 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 force=force, virial=virial)
143
144 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
145
146 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
147
148 DO ikind = 1, SIZE(atomic_kind_set)
149
150 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
151 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
152
153 IF (.NOT. ASSOCIATED(gth_potential)) cycle
154 CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
155
156 IF (lppl <= 0) cycle
157
158 ni = ncoset(2*lppl - 2)
159 ALLOCATE (hab(ni, 1), pab(ni, 1))
160 pab = 0._dp
161
162 CALL reallocate(cores, 1, natom_of_kind)
163 npme = 0
164 cores = 0
165
166 ! prepare core function
167 DO j = 1, lppl
168 SELECT CASE (j)
169 CASE (1)
170 pab(1, 1) = cexp_ppl(1)
171 CASE (2)
172 n = coset(2, 0, 0)
173 pab(n, 1) = cexp_ppl(2)
174 n = coset(0, 2, 0)
175 pab(n, 1) = cexp_ppl(2)
176 n = coset(0, 0, 2)
177 pab(n, 1) = cexp_ppl(2)
178 CASE (3)
179 n = coset(4, 0, 0)
180 pab(n, 1) = cexp_ppl(3)
181 n = coset(0, 4, 0)
182 pab(n, 1) = cexp_ppl(3)
183 n = coset(0, 0, 4)
184 pab(n, 1) = cexp_ppl(3)
185 n = coset(2, 2, 0)
186 pab(n, 1) = 2._dp*cexp_ppl(3)
187 n = coset(2, 0, 2)
188 pab(n, 1) = 2._dp*cexp_ppl(3)
189 n = coset(0, 2, 2)
190 pab(n, 1) = 2._dp*cexp_ppl(3)
191 CASE (4)
192 n = coset(6, 0, 0)
193 pab(n, 1) = cexp_ppl(4)
194 n = coset(0, 6, 0)
195 pab(n, 1) = cexp_ppl(4)
196 n = coset(0, 0, 6)
197 pab(n, 1) = cexp_ppl(4)
198 n = coset(4, 2, 0)
199 pab(n, 1) = 3._dp*cexp_ppl(4)
200 n = coset(4, 0, 2)
201 pab(n, 1) = 3._dp*cexp_ppl(4)
202 n = coset(2, 4, 0)
203 pab(n, 1) = 3._dp*cexp_ppl(4)
204 n = coset(2, 0, 4)
205 pab(n, 1) = 3._dp*cexp_ppl(4)
206 n = coset(0, 4, 2)
207 pab(n, 1) = 3._dp*cexp_ppl(4)
208 n = coset(0, 2, 4)
209 pab(n, 1) = 3._dp*cexp_ppl(4)
210 n = coset(2, 2, 2)
211 pab(n, 1) = 6._dp*cexp_ppl(4)
212 CASE DEFAULT
213 cpabort("")
214 END SELECT
215 END DO
216
217 DO iatom = 1, natom_of_kind
218 atom_a = atom_list(iatom)
219 ra(:) = pbc(particle_set(atom_a)%r, cell)
220 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
221 ! replicated realspace grid, split the atoms up between procs
222 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
223 npme = npme + 1
224 cores(npme) = iatom
225 END IF
226 ELSE
227 npme = npme + 1
228 cores(npme) = iatom
229 END IF
230 END DO
231
232 DO j = 1, npme
233
234 iatom = cores(j)
235 atom_a = atom_list(iatom)
236 ra(:) = pbc(particle_set(atom_a)%r, cell)
237 hab(:, 1) = 0.0_dp
238 force_a(:) = 0.0_dp
239 force_b(:) = 0.0_dp
240 IF (use_virial) THEN
241 my_virial_a = 0.0_dp
242 my_virial_b = 0.0_dp
243 END IF
244 ni = 2*lppl - 2
245
246 radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
247 ra=ra, rb=ra, rp=ra, &
248 zetp=alpha, eps=eps_rho_rspace, &
249 pab=pab, o1=0, o2=0, & ! without map_consistent
250 prefactor=1.0_dp, cutoff=1.0_dp)
251
252 CALL integrate_pgf_product(ni, alpha, 0, &
253 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
254 rs_v, hab, pab=pab, o1=0, o2=0, &
255 radius=radius, &
256 calculate_forces=.true., force_a=force_a, &
257 force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
258 my_virial_b=my_virial_b, use_subpatch=.true., subpatch_pattern=0)
259
260 force(ikind)%gth_ppl(:, iatom) = &
261 force(ikind)%gth_ppl(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
262
263 IF (use_virial) THEN
264 virial%pv_ppl = virial%pv_ppl + my_virial_a*rho_rspace%pw_grid%dvol
265 virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
266 cpabort("Virial not debuged for CORE_PPL")
267 END IF
268 END DO
269
270 DEALLOCATE (hab, pab)
271
272 END DO
273
274 DEALLOCATE (cores)
275
276 CALL timestop(handle)
277
278 END SUBROUTINE integrate_ppl_rspace
279
280! **************************************************************************************************
281!> \brief computes the forces/virial due to the nlcc pseudopotential
282!> \param rho_rspace ...
283!> \param qs_env ...
284! **************************************************************************************************
285 SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env)
286 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_rspace
287 TYPE(qs_environment_type), POINTER :: qs_env
288
289 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_rho_nlcc'
290
291 INTEGER :: atom_a, handle, iatom, iexp_nlcc, ikind, &
292 ithread, j, n, natom, nc, nexp_nlcc, &
293 ni, npme, nthread
294 INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc
295 LOGICAL :: nlcc, use_virial
296 REAL(kind=dp) :: alpha, eps_rho_rspace, radius
297 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
298 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
299 REAL(kind=dp), DIMENSION(:), POINTER :: alpha_nlcc
300 REAL(kind=dp), DIMENSION(:, :), POINTER :: cval_nlcc, hab, pab
301 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
302 TYPE(cell_type), POINTER :: cell
303 TYPE(dft_control_type), POINTER :: dft_control
304 TYPE(gth_potential_type), POINTER :: gth_potential
305 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
306 TYPE(pw_env_type), POINTER :: pw_env
307 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
308 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
309 TYPE(realspace_grid_type), POINTER :: rs_v
310 TYPE(virial_type), POINTER :: virial
311
312 CALL timeset(routinen, handle)
313
314 NULLIFY (pw_env, cores)
315
316 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
317 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
318
319 CALL transfer_pw2rs(rs_v, rho_rspace)
320
321 CALL get_qs_env(qs_env=qs_env, &
322 atomic_kind_set=atomic_kind_set, &
323 qs_kind_set=qs_kind_set, &
324 cell=cell, &
325 dft_control=dft_control, &
326 particle_set=particle_set, &
327 pw_env=pw_env, &
328 force=force, virial=virial)
329
330 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
331
332 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
333
334 DO ikind = 1, SIZE(atomic_kind_set)
335
336 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
337 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
338
339 IF (.NOT. ASSOCIATED(gth_potential)) cycle
340 CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
341 alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
342
343 IF (.NOT. nlcc) cycle
344
345 DO iexp_nlcc = 1, nexp_nlcc
346
347 alpha = alpha_nlcc(iexp_nlcc)
348 nc = nct_nlcc(iexp_nlcc)
349
350 ni = ncoset(2*nc - 2)
351
352 nthread = 1
353 ithread = 0
354
355 ALLOCATE (hab(ni, 1), pab(ni, 1))
356 pab = 0._dp
357
358 CALL reallocate(cores, 1, natom)
359 npme = 0
360 cores = 0
361
362 ! prepare core function
363 DO j = 1, nc
364 SELECT CASE (j)
365 CASE (1)
366 pab(1, 1) = cval_nlcc(1, iexp_nlcc)
367 CASE (2)
368 n = coset(2, 0, 0)
369 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
370 n = coset(0, 2, 0)
371 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
372 n = coset(0, 0, 2)
373 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
374 CASE (3)
375 n = coset(4, 0, 0)
376 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
377 n = coset(0, 4, 0)
378 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
379 n = coset(0, 0, 4)
380 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
381 n = coset(2, 2, 0)
382 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
383 n = coset(2, 0, 2)
384 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
385 n = coset(0, 2, 2)
386 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
387 CASE (4)
388 n = coset(6, 0, 0)
389 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
390 n = coset(0, 6, 0)
391 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
392 n = coset(0, 0, 6)
393 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
394 n = coset(4, 2, 0)
395 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
396 n = coset(4, 0, 2)
397 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
398 n = coset(2, 4, 0)
399 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
400 n = coset(2, 0, 4)
401 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
402 n = coset(0, 4, 2)
403 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
404 n = coset(0, 2, 4)
405 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
406 n = coset(2, 2, 2)
407 pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
408 CASE DEFAULT
409 cpabort("")
410 END SELECT
411 END DO
412 IF (dft_control%nspins == 2) pab = pab*0.5_dp
413
414 DO iatom = 1, natom
415 atom_a = atom_list(iatom)
416 ra(:) = pbc(particle_set(atom_a)%r, cell)
417 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
418 ! replicated realspace grid, split the atoms up between procs
419 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
420 npme = npme + 1
421 cores(npme) = iatom
422 END IF
423 ELSE
424 npme = npme + 1
425 cores(npme) = iatom
426 END IF
427 END DO
428
429 DO j = 1, npme
430
431 iatom = cores(j)
432 atom_a = atom_list(iatom)
433 ra(:) = pbc(particle_set(atom_a)%r, cell)
434 hab(:, 1) = 0.0_dp
435 force_a(:) = 0.0_dp
436 force_b(:) = 0.0_dp
437 IF (use_virial) THEN
438 my_virial_a = 0.0_dp
439 my_virial_b = 0.0_dp
440 END IF
441 ni = 2*nc - 2
442
443 radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
444 ra=ra, rb=ra, rp=ra, &
445 zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
446 pab=pab, o1=0, o2=0, & ! without map_consistent
447 prefactor=1.0_dp, cutoff=1.0_dp)
448
449 CALL integrate_pgf_product(ni, 1/(2*alpha**2), 0, &
450 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
451 rs_v, hab, pab=pab, o1=0, o2=0, &
452 radius=radius, &
453 calculate_forces=.true., force_a=force_a, &
454 force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
455 my_virial_b=my_virial_b, use_subpatch=.true., subpatch_pattern=0)
456
457 force(ikind)%gth_nlcc(:, iatom) = &
458 force(ikind)%gth_nlcc(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
459
460 IF (use_virial) THEN
461 virial%pv_nlcc = virial%pv_nlcc + my_virial_a*rho_rspace%pw_grid%dvol
462 virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
463 END IF
464 END DO
465
466 DEALLOCATE (hab, pab)
467
468 END DO
469
470 END DO
471
472 DEALLOCATE (cores)
473
474 CALL timestop(handle)
475
476 END SUBROUTINE integrate_rho_nlcc
477
478! **************************************************************************************************
479!> \brief computes the forces/virial due to the ionic cores with a potential on
480!> grid
481!> \param v_rspace ...
482!> \param qs_env ...
483!> \param atecc ...
484! **************************************************************************************************
485 SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env, atecc)
486 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
487 TYPE(qs_environment_type), POINTER :: qs_env
488 REAL(kind=dp), DIMENSION(:), OPTIONAL :: atecc
489
490 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_core_rspace'
491
492 INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
493 natom_of_kind, npme
494 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
495 LOGICAL :: paw_atom, skip_fcore, use_virial
496 REAL(kind=dp) :: alpha_core_charge, ccore_charge, &
497 eps_rho_rspace, radius
498 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
499 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
500 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab
501 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
502 TYPE(atprop_type), POINTER :: atprop
503 TYPE(cell_type), POINTER :: cell
504 TYPE(dft_control_type), POINTER :: dft_control
505 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
506 TYPE(pw_env_type), POINTER :: pw_env
507 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
508 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
509 TYPE(realspace_grid_type), POINTER :: rs_v
510 TYPE(virial_type), POINTER :: virial
511
512 CALL timeset(routinen, handle)
513 NULLIFY (virial, force, atprop, dft_control)
514
515 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
516
517 !If gapw, check for gpw kinds
518 skip_fcore = .false.
519 IF (dft_control%qs_control%gapw) THEN
520 IF (.NOT. dft_control%qs_control%gapw_control%nopaw_as_gpw) skip_fcore = .true.
521 END IF
522
523 IF (.NOT. skip_fcore) THEN
524 NULLIFY (pw_env)
525 ALLOCATE (cores(1))
526 ALLOCATE (hab(1, 1))
527 ALLOCATE (pab(1, 1))
528
529 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
530 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
531
532 CALL transfer_pw2rs(rs_v, v_rspace)
533
534 CALL get_qs_env(qs_env=qs_env, &
535 atomic_kind_set=atomic_kind_set, &
536 qs_kind_set=qs_kind_set, &
537 cell=cell, &
538 dft_control=dft_control, &
539 particle_set=particle_set, &
540 pw_env=pw_env, &
541 force=force, &
542 virial=virial, &
543 atprop=atprop)
544
545 ! atomic energy contributions
546 natom = SIZE(particle_set)
547 IF (ASSOCIATED(atprop)) THEN
548 CALL atprop_array_init(atprop%ateb, natom)
549 END IF
550
551 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
552
553 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
554
555 DO ikind = 1, SIZE(atomic_kind_set)
556
557 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
558 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
559 alpha_core_charge=alpha_core_charge, &
560 ccore_charge=ccore_charge)
561
562 IF (dft_control%qs_control%gapw .AND. paw_atom) cycle
563
564 pab(1, 1) = -ccore_charge
565 IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
566
567 CALL reallocate(cores, 1, natom_of_kind)
568 npme = 0
569 cores = 0
570
571 DO iatom = 1, natom_of_kind
572 atom_a = atom_list(iatom)
573 ra(:) = pbc(particle_set(atom_a)%r, cell)
574 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
575 ! replicated realspace grid, split the atoms up between procs
576 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
577 npme = npme + 1
578 cores(npme) = iatom
579 END IF
580 ELSE
581 npme = npme + 1
582 cores(npme) = iatom
583 END IF
584 END DO
585
586 DO j = 1, npme
587
588 iatom = cores(j)
589 atom_a = atom_list(iatom)
590 ra(:) = pbc(particle_set(atom_a)%r, cell)
591 hab(1, 1) = 0.0_dp
592 force_a(:) = 0.0_dp
593 force_b(:) = 0.0_dp
594 IF (use_virial) THEN
595 my_virial_a = 0.0_dp
596 my_virial_b = 0.0_dp
597 END IF
598
599 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
600 ra=ra, rb=ra, rp=ra, &
601 zetp=alpha_core_charge, eps=eps_rho_rspace, &
602 pab=pab, o1=0, o2=0, & ! without map_consistent
603 prefactor=1.0_dp, cutoff=1.0_dp)
604
605 CALL integrate_pgf_product(0, alpha_core_charge, 0, &
606 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
607 rs_v, hab, pab=pab, o1=0, o2=0, &
608 radius=radius, &
609 calculate_forces=.true., force_a=force_a, &
610 force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
611 my_virial_b=my_virial_b, use_subpatch=.true., subpatch_pattern=0)
612
613 IF (ASSOCIATED(force)) THEN
614 force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
615 END IF
616 IF (use_virial) THEN
617 virial%pv_ehartree = virial%pv_ehartree + my_virial_a
618 virial%pv_virial = virial%pv_virial + my_virial_a
619 END IF
620 IF (ASSOCIATED(atprop)) THEN
621 atprop%ateb(atom_a) = atprop%ateb(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
622 END IF
623 IF (PRESENT(atecc)) THEN
624 atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
625 END IF
626
627 END DO
628
629 END DO
630
631 DEALLOCATE (hab, pab, cores)
632
633 END IF
634
635 CALL timestop(handle)
636
637 END SUBROUTINE integrate_v_core_rspace
638
639! **************************************************************************************************
640!> \brief computes the overlap of a set of Gaussians with a potential on grid
641!> \param v_rspace ...
642!> \param qs_env ...
643!> \param alpha ...
644!> \param ccore ...
645!> \param atecc ...
646! **************************************************************************************************
647 SUBROUTINE integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
648 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
649 TYPE(qs_environment_type), POINTER :: qs_env
650 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: alpha, ccore
651 REAL(kind=dp), DIMENSION(:) :: atecc
652
653 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_gaussian_rspace'
654
655 INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
656 natom_of_kind, npme
657 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
658 REAL(kind=dp) :: alpha_core_charge, eps_rho_rspace, radius
659 REAL(kind=dp), DIMENSION(3) :: ra
660 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab
661 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
662 TYPE(cell_type), POINTER :: cell
663 TYPE(dft_control_type), POINTER :: dft_control
664 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
665 TYPE(pw_env_type), POINTER :: pw_env
666 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
667 TYPE(realspace_grid_type), POINTER :: rs_v
668
669 CALL timeset(routinen, handle)
670
671 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
672
673 !If gapw, check for gpw kinds
674 cpassert(.NOT. dft_control%qs_control%gapw)
675
676 NULLIFY (pw_env)
677 ALLOCATE (cores(1))
678 ALLOCATE (hab(1, 1))
679 ALLOCATE (pab(1, 1))
680
681 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
682 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
683
684 CALL transfer_pw2rs(rs_v, v_rspace)
685
686 CALL get_qs_env(qs_env=qs_env, &
687 atomic_kind_set=atomic_kind_set, &
688 qs_kind_set=qs_kind_set, &
689 cell=cell, &
690 dft_control=dft_control, &
691 particle_set=particle_set, &
692 pw_env=pw_env)
693
694 ! atomic energy contributions
695 natom = SIZE(particle_set)
696 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
697
698 DO ikind = 1, SIZE(atomic_kind_set)
699
700 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
701 pab(1, 1) = -ccore(ikind)
702 alpha_core_charge = alpha(ikind)
703 IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
704
705 CALL reallocate(cores, 1, natom_of_kind)
706 npme = 0
707 cores = 0
708
709 DO iatom = 1, natom_of_kind
710 atom_a = atom_list(iatom)
711 ra(:) = pbc(particle_set(atom_a)%r, cell)
712 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
713 ! replicated realspace grid, split the atoms up between procs
714 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
715 npme = npme + 1
716 cores(npme) = iatom
717 END IF
718 ELSE
719 npme = npme + 1
720 cores(npme) = iatom
721 END IF
722 END DO
723
724 DO j = 1, npme
725
726 iatom = cores(j)
727 atom_a = atom_list(iatom)
728 ra(:) = pbc(particle_set(atom_a)%r, cell)
729 hab(1, 1) = 0.0_dp
730
731 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
732 ra=ra, rb=ra, rp=ra, &
733 zetp=alpha_core_charge, eps=eps_rho_rspace, &
734 pab=pab, o1=0, o2=0, & ! without map_consistent
735 prefactor=1.0_dp, cutoff=1.0_dp)
736
737 CALL integrate_pgf_product(0, alpha_core_charge, 0, &
738 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
739 rs_v, hab, pab=pab, o1=0, o2=0, &
740 radius=radius, calculate_forces=.false., &
741 use_subpatch=.true., subpatch_pattern=0)
742 atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
743
744 END DO
745
746 END DO
747
748 DEALLOCATE (hab, pab, cores)
749
750 CALL timestop(handle)
751
752 END SUBROUTINE integrate_v_gaussian_rspace
753! **************************************************************************************************
754!> \brief computes integrals of product of v_rspace times a one-center function
755!> required for LRIGPW
756!> \param v_rspace ...
757!> \param qs_env ...
758!> \param int_res ...
759!> \param calculate_forces ...
760!> \param basis_type ...
761!> \param atomlist ...
762!> \author Dorothea Golze
763! **************************************************************************************************
764 SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, &
765 calculate_forces, basis_type, atomlist)
766 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
767 TYPE(qs_environment_type), POINTER :: qs_env
768 TYPE(lri_kind_type), DIMENSION(:), POINTER :: int_res
769 LOGICAL, INTENT(IN) :: calculate_forces
770 CHARACTER(len=*), INTENT(IN) :: basis_type
771 INTEGER, DIMENSION(:), OPTIONAL :: atomlist
772
773 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_rspace_one_center'
774
775 INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, m1, maxco, &
776 maxsgf_set, my_pos, na1, natom_of_kind, ncoa, nkind, nseta, offset, sgfa
777 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
778 nsgf_seta
779 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
780 LOGICAL :: use_virial
781 LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it
782 REAL(kind=dp) :: eps_rho_rspace, radius
783 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
784 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
785 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a
786 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab, rpgfa, sphi_a, work_f, work_i, &
787 zeta
788 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
789 TYPE(cell_type), POINTER :: cell
790 TYPE(dft_control_type), POINTER :: dft_control
791 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
792 TYPE(gto_basis_set_type), POINTER :: lri_basis_set
793 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
794 TYPE(pw_env_type), POINTER :: pw_env
795 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
796 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
797 TYPE(realspace_grid_type), POINTER :: rs_grid
798 TYPE(virial_type), POINTER :: virial
799
800 CALL timeset(routinen, handle)
801
802 NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
803 first_sgfa, gridlevel_info, hab, la_max, la_min, lri_basis_set, &
804 npgfa, nsgf_seta, pab, particle_set, pw_env, rpgfa, &
805 rs_grid, rs_v, virial, set_radius_a, sphi_a, work_f, &
806 work_i, zeta)
807
808 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
809
810 CALL pw_env_get(pw_env, rs_grids=rs_v)
811 DO i = 1, SIZE(rs_v)
812 CALL rs_grid_zero(rs_v(i))
813 END DO
814
815 gridlevel_info => pw_env%gridlevel_info
816
817 CALL potential_pw2rs(rs_v, v_rspace, pw_env)
818
819 CALL get_qs_env(qs_env=qs_env, &
820 atomic_kind_set=atomic_kind_set, &
821 qs_kind_set=qs_kind_set, &
822 cell=cell, &
823 dft_control=dft_control, &
824 nkind=nkind, &
825 particle_set=particle_set, &
826 pw_env=pw_env, &
827 virial=virial)
828
829 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
830
831 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
832
833 offset = 0
834 my_pos = v_rspace%pw_grid%para%group%mepos
835 group_size = v_rspace%pw_grid%para%group%num_pe
836
837 DO ikind = 1, nkind
838
839 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
840 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
841 CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
842 first_sgf=first_sgfa, &
843 lmax=la_max, &
844 lmin=la_min, &
845 maxco=maxco, &
846 maxsgf_set=maxsgf_set, &
847 npgf=npgfa, &
848 nset=nseta, &
849 nsgf_set=nsgf_seta, &
850 pgf_radius=rpgfa, &
851 set_radius=set_radius_a, &
852 sphi=sphi_a, &
853 zet=zeta)
854
855 ALLOCATE (hab(maxco, 1), pab(maxco, 1))
856 hab = 0._dp
857 pab(:, 1) = 0._dp
858
859 DO iatom = 1, natom_of_kind
860
861 atom_a = atom_list(iatom)
862 IF (PRESENT(atomlist)) THEN
863 IF (atomlist(atom_a) == 0) cycle
864 END IF
865 ra(:) = pbc(particle_set(atom_a)%r, cell)
866 force_a(:) = 0._dp
867 force_b(:) = 0._dp
868 my_virial_a(:, :) = 0._dp
869 my_virial_b(:, :) = 0._dp
870
871 m1 = maxval(npgfa(1:nseta))
872 ALLOCATE (map_it(m1))
873
874 DO iset = 1, nseta
875 !
876 map_it = .false.
877 DO ipgf = 1, npgfa(iset)
878 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
879 rs_grid => rs_v(igrid_level)
880 map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
881 END DO
882 offset = offset + 1
883 !
884 IF (any(map_it(1:npgfa(iset)))) THEN
885 sgfa = first_sgfa(1, iset)
886 ncoa = npgfa(iset)*ncoset(la_max(iset))
887 hab(:, 1) = 0._dp
888 ALLOCATE (work_i(nsgf_seta(iset), 1))
889 work_i = 0.0_dp
890
891 ! get fit coefficients for forces
892 IF (calculate_forces) THEN
893 m1 = sgfa + nsgf_seta(iset) - 1
894 ALLOCATE (work_f(nsgf_seta(iset), 1))
895 work_f(1:nsgf_seta(iset), 1) = int_res(ikind)%acoef(iatom, sgfa:m1)
896 CALL dgemm("N", "N", ncoa, 1, nsgf_seta(iset), 1.0_dp, sphi_a(1, sgfa), &
897 SIZE(sphi_a, 1), work_f(1, 1), SIZE(work_f, 1), 0.0_dp, pab(1, 1), &
898 SIZE(pab, 1))
899 DEALLOCATE (work_f)
900 END IF
901
902 DO ipgf = 1, npgfa(iset)
903 na1 = (ipgf - 1)*ncoset(la_max(iset))
904 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
905 rs_grid => rs_v(igrid_level)
906
907 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
908 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
909 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
910 prefactor=1.0_dp, cutoff=1.0_dp)
911
912 IF (map_it(ipgf)) THEN
913 IF (.NOT. calculate_forces) THEN
914 CALL integrate_pgf_product(la_max=la_max(iset), &
915 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
916 lb_max=0, zetb=0.0_dp, lb_min=0, &
917 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
918 rsgrid=rs_grid, &
919 hab=hab, o1=na1, o2=0, radius=radius, &
920 calculate_forces=calculate_forces)
921 ELSE
922 CALL integrate_pgf_product(la_max=la_max(iset), &
923 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
924 lb_max=0, zetb=0.0_dp, lb_min=0, &
925 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
926 rsgrid=rs_grid, &
927 hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
928 calculate_forces=calculate_forces, &
929 force_a=force_a, force_b=force_b, &
930 use_virial=use_virial, &
931 my_virial_a=my_virial_a, my_virial_b=my_virial_b)
932 END IF
933 END IF
934 END DO
935 ! contract hab
936 CALL dgemm("T", "N", nsgf_seta(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
937 SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work_i(1, 1), SIZE(work_i, 1))
938
939 int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) = &
940 int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) + work_i(1:nsgf_seta(iset), 1)
941 DEALLOCATE (work_i)
942 END IF
943 END DO
944 !
945 IF (calculate_forces) THEN
946 int_res(ikind)%v_dfdr(iatom, :) = int_res(ikind)%v_dfdr(iatom, :) + force_a(:)
947 IF (use_virial) THEN
948 virial%pv_lrigpw = virial%pv_lrigpw + my_virial_a
949 virial%pv_virial = virial%pv_virial + my_virial_a
950 END IF
951 END IF
952
953 DEALLOCATE (map_it)
954
955 END DO
956
957 DEALLOCATE (hab, pab)
958 END DO
959
960 CALL timestop(handle)
961
962 END SUBROUTINE integrate_v_rspace_one_center
963
964! **************************************************************************************************
965!> \brief computes integrals of product of v_rspace times the diagonal block basis functions
966!> required for LRIGPW with exact 1c terms
967!> \param v_rspace ...
968!> \param ksmat ...
969!> \param pmat ...
970!> \param qs_env ...
971!> \param calculate_forces ...
972!> \param basis_type ...
973!> \author JGH
974! **************************************************************************************************
975 SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
976 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
977 TYPE(dbcsr_type), INTENT(INOUT) :: ksmat, pmat
978 TYPE(qs_environment_type), POINTER :: qs_env
979 LOGICAL, INTENT(IN) :: calculate_forces
980 CHARACTER(len=*), INTENT(IN) :: basis_type
981
982 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_rspace_diagonal'
983
984 INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
985 m1, maxco, maxsgf_set, my_pos, na1, na2, natom_of_kind, nb1, nb2, ncoa, ncob, nkind, &
986 nseta, nsgfa, offset, sgfa, sgfb
987 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
988 nsgf_seta
989 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
990 LOGICAL :: found, use_virial
991 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
992 REAL(kind=dp) :: eps_rho_rspace, radius, zetp
993 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
994 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
995 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a
996 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, hab, hmat, p_block, pab, pblk, &
997 rpgfa, sphi_a, work, zeta
998 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
999 TYPE(cell_type), POINTER :: cell
1000 TYPE(dft_control_type), POINTER :: dft_control
1001 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1002 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1003 TYPE(mp_para_env_type), POINTER :: para_env
1004 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1005 TYPE(pw_env_type), POINTER :: pw_env
1006 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1007 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1008 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1009 TYPE(realspace_grid_type), POINTER :: rs_grid
1010 TYPE(virial_type), POINTER :: virial
1011
1012 CALL timeset(routinen, handle)
1013
1014 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1015 CALL pw_env_get(pw_env, rs_grids=rs_v)
1016 CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1017
1018 gridlevel_info => pw_env%gridlevel_info
1019
1020 CALL get_qs_env(qs_env=qs_env, &
1021 atomic_kind_set=atomic_kind_set, &
1022 qs_kind_set=qs_kind_set, &
1023 cell=cell, &
1024 dft_control=dft_control, &
1025 nkind=nkind, &
1026 particle_set=particle_set, &
1027 force=force, &
1028 virial=virial, &
1029 para_env=para_env)
1030
1031 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1032 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1033
1034 offset = 0
1035 my_pos = v_rspace%pw_grid%para%group%mepos
1036 group_size = v_rspace%pw_grid%para%group%num_pe
1037
1038 DO ikind = 1, nkind
1039
1040 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
1041 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1042 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1043 lmax=la_max, lmin=la_min, maxco=maxco, maxsgf_set=maxsgf_set, &
1044 npgf=npgfa, nset=nseta, nsgf_set=nsgf_seta, nsgf=nsgfa, &
1045 first_sgf=first_sgfa, pgf_radius=rpgfa, set_radius=set_radius_a, &
1046 sphi=sphi_a, zet=zeta)
1047
1048 ALLOCATE (hab(maxco, maxco), work(maxco, maxsgf_set), hmat(nsgfa, nsgfa))
1049 IF (calculate_forces) ALLOCATE (pab(maxco, maxco), pblk(nsgfa, nsgfa))
1050
1051 DO iatom = 1, natom_of_kind
1052 atom_a = atom_list(iatom)
1053 ra(:) = pbc(particle_set(atom_a)%r, cell)
1054 hmat = 0.0_dp
1055 IF (calculate_forces) THEN
1056 CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, block=p_block, found=found)
1057 IF (found) THEN
1058 pblk(1:nsgfa, 1:nsgfa) = p_block(1:nsgfa, 1:nsgfa)
1059 ELSE
1060 pblk = 0.0_dp
1061 END IF
1062 CALL para_env%sum(pblk)
1063 force_a(:) = 0._dp
1064 force_b(:) = 0._dp
1065 IF (use_virial) THEN
1066 my_virial_a = 0.0_dp
1067 my_virial_b = 0.0_dp
1068 END IF
1069 END IF
1070 m1 = maxval(npgfa(1:nseta))
1071 ALLOCATE (map_it2(m1, m1))
1072 DO iset = 1, nseta
1073 sgfa = first_sgfa(1, iset)
1074 ncoa = npgfa(iset)*ncoset(la_max(iset))
1075 DO jset = 1, nseta
1076 sgfb = first_sgfa(1, jset)
1077 ncob = npgfa(jset)*ncoset(la_max(jset))
1078 !
1079 map_it2 = .false.
1080 DO ipgf = 1, npgfa(iset)
1081 DO jpgf = 1, npgfa(jset)
1082 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1083 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1084 rs_grid => rs_v(igrid_level)
1085 map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
1086 END DO
1087 END DO
1088 offset = offset + 1
1089 !
1090 IF (any(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
1091 hab(:, :) = 0._dp
1092 IF (calculate_forces) THEN
1093 CALL dgemm("N", "N", ncoa, nsgf_seta(jset), nsgf_seta(iset), &
1094 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1095 pblk(sgfa, sgfb), SIZE(pblk, 1), &
1096 0.0_dp, work(1, 1), SIZE(work, 1))
1097 CALL dgemm("N", "T", ncoa, ncob, nsgf_seta(jset), &
1098 1.0_dp, work(1, 1), SIZE(work, 1), &
1099 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1100 0.0_dp, pab(1, 1), SIZE(pab, 1))
1101 END IF
1102
1103 DO ipgf = 1, npgfa(iset)
1104 na1 = (ipgf - 1)*ncoset(la_max(iset))
1105 na2 = ipgf*ncoset(la_max(iset))
1106 DO jpgf = 1, npgfa(jset)
1107 nb1 = (jpgf - 1)*ncoset(la_max(jset))
1108 nb2 = jpgf*ncoset(la_max(jset))
1109 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1110 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1111 rs_grid => rs_v(igrid_level)
1112
1113 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1114 lb_min=la_min(jset), lb_max=la_max(jset), &
1115 ra=ra, rb=ra, rp=ra, &
1116 zetp=zetp, eps=eps_rho_rspace, &
1117 prefactor=1.0_dp, cutoff=1.0_dp)
1118
1119 IF (map_it2(ipgf, jpgf)) THEN
1120 IF (calculate_forces) THEN
1121 CALL integrate_pgf_product( &
1122 la_max(iset), zeta(ipgf, iset), la_min(iset), &
1123 la_max(jset), zeta(jpgf, jset), la_min(jset), &
1124 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1125 rsgrid=rs_v(igrid_level), &
1126 hab=hab, pab=pab, o1=na1, o2=nb1, &
1127 radius=radius, &
1128 calculate_forces=.true., &
1129 force_a=force_a, force_b=force_b, &
1130 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1131 ELSE
1132 CALL integrate_pgf_product( &
1133 la_max(iset), zeta(ipgf, iset), la_min(iset), &
1134 la_max(jset), zeta(jpgf, jset), la_min(jset), &
1135 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1136 rsgrid=rs_v(igrid_level), &
1137 hab=hab, o1=na1, o2=nb1, &
1138 radius=radius, &
1139 calculate_forces=.false.)
1140 END IF
1141 END IF
1142 END DO
1143 END DO
1144 ! contract hab
1145 CALL dgemm("N", "N", ncoa, nsgf_seta(jset), ncob, &
1146 1.0_dp, hab(1, 1), SIZE(hab, 1), &
1147 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1148 0.0_dp, work(1, 1), SIZE(work, 1))
1149 CALL dgemm("T", "N", nsgf_seta(iset), nsgf_seta(jset), ncoa, &
1150 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1151 work(1, 1), SIZE(work, 1), &
1152 1.0_dp, hmat(sgfa, sgfb), SIZE(hmat, 1))
1153 END IF
1154 END DO
1155 END DO
1156 DEALLOCATE (map_it2)
1157 ! update KS matrix
1158 CALL para_env%sum(hmat)
1159 CALL dbcsr_get_block_p(matrix=ksmat, row=atom_a, col=atom_a, block=h_block, found=found)
1160 IF (found) THEN
1161 h_block(1:nsgfa, 1:nsgfa) = h_block(1:nsgfa, 1:nsgfa) + hmat(1:nsgfa, 1:nsgfa)
1162 END IF
1163 IF (calculate_forces) THEN
1164 force(ikind)%rho_elec(:, iatom) = force(ikind)%rho_elec(:, iatom) + 2.0_dp*force_a(:)
1165 IF (use_virial) THEN
1166 IF (use_virial .AND. calculate_forces) THEN
1167 virial%pv_lrigpw = virial%pv_lrigpw + 2.0_dp*my_virial_a
1168 virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a
1169 END IF
1170 END IF
1171 END IF
1172 END DO
1173 DEALLOCATE (hab, work, hmat)
1174 IF (calculate_forces) DEALLOCATE (pab, pblk)
1175 END DO
1176
1177 CALL timestop(handle)
1178
1179 END SUBROUTINE integrate_v_rspace_diagonal
1180
1181! **************************************************************************************************
1182!> \brief computes integrals of product of v_rspace times a basis function (vector function)
1183!> and possible forces
1184!> \param qs_env ...
1185!> \param v_rspace ...
1186!> \param f_coef ...
1187!> \param f_integral ...
1188!> \param calculate_forces ...
1189!> \param basis_type ...
1190!> \author JGH [8.2024]
1191! **************************************************************************************************
1192 SUBROUTINE integrate_function(qs_env, v_rspace, f_coef, f_integral, calculate_forces, basis_type)
1193 TYPE(qs_environment_type), POINTER :: qs_env
1194 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
1195 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: f_coef
1196 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: f_integral
1197 LOGICAL, INTENT(IN) :: calculate_forces
1198 CHARACTER(len=*), INTENT(IN) :: basis_type
1199
1200 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_function'
1201
1202 INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, &
1203 maxsgf_set, my_pos, na1, natom, ncoa, nkind, nseta, offset, sgfa
1204 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1205 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
1206 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1207 LOGICAL :: use_virial
1208 REAL(kind=dp) :: eps_rho_rspace, radius
1209 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
1210 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
1211 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab, sphi_a, work, zeta
1212 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1213 TYPE(cell_type), POINTER :: cell
1214 TYPE(dft_control_type), POINTER :: dft_control
1215 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1216 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1217 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1218 TYPE(pw_env_type), POINTER :: pw_env
1219 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1220 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1221 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1222 TYPE(realspace_grid_type), POINTER :: rs_grid
1223 TYPE(virial_type), POINTER :: virial
1224
1225 CALL timeset(routinen, handle)
1226
1227 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1228 gridlevel_info => pw_env%gridlevel_info
1229
1230 CALL pw_env_get(pw_env, rs_grids=rs_v)
1231 CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1232
1233 CALL get_qs_env(qs_env=qs_env, &
1234 atomic_kind_set=atomic_kind_set, &
1235 qs_kind_set=qs_kind_set, &
1236 force=force, &
1237 cell=cell, &
1238 dft_control=dft_control, &
1239 nkind=nkind, &
1240 natom=natom, &
1241 particle_set=particle_set, &
1242 pw_env=pw_env, &
1243 virial=virial)
1244 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1245
1246 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1247 IF (use_virial) THEN
1248 cpabort("Virial NYA")
1249 END IF
1250
1251 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1252
1253 CALL get_qs_kind_set(qs_kind_set, &
1254 maxco=maxco, maxsgf_set=maxsgf_set, basis_type=basis_type)
1255 ALLOCATE (hab(maxco, 1), pab(maxco, 1), work(maxco, 1))
1256
1257 offset = 0
1258 my_pos = v_rspace%pw_grid%para%group%mepos
1259 group_size = v_rspace%pw_grid%para%group%num_pe
1260
1261 DO iatom = 1, natom
1262 ikind = particle_set(iatom)%atomic_kind%kind_number
1263 atom_a = atom_of_kind(iatom)
1264 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1265 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1266 first_sgf=first_sgfa, &
1267 lmax=la_max, &
1268 lmin=la_min, &
1269 npgf=npgfa, &
1270 nset=nseta, &
1271 nsgf_set=nsgfa, &
1272 sphi=sphi_a, &
1273 zet=zeta)
1274 ra(:) = pbc(particle_set(iatom)%r, cell)
1275
1276 force_a(:) = 0._dp
1277 force_b(:) = 0._dp
1278 my_virial_a(:, :) = 0._dp
1279 my_virial_b(:, :) = 0._dp
1280
1281 DO iset = 1, nseta
1282
1283 ncoa = npgfa(iset)*ncoset(la_max(iset))
1284 sgfa = first_sgfa(1, iset)
1285
1286 hab = 0._dp
1287 pab = 0._dp
1288
1289 DO i = 1, nsgfa(iset)
1290 work(i, 1) = f_coef(offset + i)
1291 END DO
1292
1293 CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
1294 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1295 work(1, 1), SIZE(work, 1), &
1296 0.0_dp, pab(1, 1), SIZE(pab, 1))
1297
1298 DO ipgf = 1, npgfa(iset)
1299
1300 na1 = (ipgf - 1)*ncoset(la_max(iset))
1301
1302 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
1303 rs_grid => rs_v(igrid_level)
1304
1305 IF (map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)) THEN
1306 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1307 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
1308 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
1309 prefactor=1.0_dp, cutoff=1.0_dp)
1310
1311 IF (calculate_forces) THEN
1312 CALL integrate_pgf_product(la_max=la_max(iset), &
1313 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1314 lb_max=0, zetb=0.0_dp, lb_min=0, &
1315 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1316 rsgrid=rs_grid, &
1317 hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
1318 calculate_forces=.true., &
1319 force_a=force_a, force_b=force_b, &
1320 use_virial=use_virial, &
1321 my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1322 ELSE
1323 CALL integrate_pgf_product(la_max=la_max(iset), &
1324 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1325 lb_max=0, zetb=0.0_dp, lb_min=0, &
1326 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1327 rsgrid=rs_grid, &
1328 hab=hab, o1=na1, o2=0, radius=radius, &
1329 calculate_forces=.false.)
1330 END IF
1331
1332 END IF
1333
1334 END DO
1335 !
1336 CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
1337 SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work(1, 1), SIZE(work, 1))
1338 DO i = 1, nsgfa(iset)
1339 f_integral(offset + i) = work(i, 1)
1340 END DO
1341
1342 offset = offset + nsgfa(iset)
1343
1344 END DO
1345
1346 IF (calculate_forces) THEN
1347 force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + force_a(:)
1348 IF (use_virial) THEN
1349 virial%pv_virial = virial%pv_virial + my_virial_a
1350 END IF
1351 END IF
1352
1353 END DO
1354
1355 DEALLOCATE (hab, pab, work)
1356
1357 CALL timestop(handle)
1358
1359 END SUBROUTINE integrate_function
1360
1361END MODULE qs_integrate_potential_single
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:206
cell_types::pbc
Definition cell_types.F:92
external_potential_types::get_potential
Definition external_potential_types.F:276
memory_utilities::reallocate
Definition memory_utilities.F:27
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition ao_util.F:209
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition atprop_types.F:87
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:624
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:679
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
gaussian_gridlevels
Definition gaussian_gridlevels.F:13
gaussian_gridlevels::gaussian_gridlevel
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Definition gaussian_gridlevels.F:148
grid_api
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
grid_api::integrate_pgf_product
subroutine, public integrate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, rsgrid, hab, pab, o1, o2, radius, calculate_forces, force_a, force_b, compute_tau, use_virial, my_virial_a, my_virial_b, hdab, hadb, a_hdab, use_subpatch, subpatch_pattern)
low level function to compute matrix elements of primitive gaussian functions
Definition grid_api.F:279
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition lri_environment_types.F:18
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_types
Definition pw_types.F:24
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_force_types
Definition qs_force_types.F:13
qs_integrate_potential_single
Build up the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_integrate_potential_single.F:30
qs_integrate_potential_single::integrate_function
subroutine, public integrate_function(qs_env, v_rspace, f_coef, f_integral, calculate_forces, basis_type)
computes integrals of product of v_rspace times a basis function (vector function) and possible force...
Definition qs_integrate_potential_single.F:1193
qs_integrate_potential_single::integrate_v_gaussian_rspace
subroutine, public integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
computes the overlap of a set of Gaussians with a potential on grid
Definition qs_integrate_potential_single.F:648
qs_integrate_potential_single::integrate_v_rspace_diagonal
subroutine, public integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
computes integrals of product of v_rspace times the diagonal block basis functions required for LRIGP...
Definition qs_integrate_potential_single.F:976
qs_integrate_potential_single::integrate_v_core_rspace
subroutine, public integrate_v_core_rspace(v_rspace, qs_env, atecc)
computes the forces/virial due to the ionic cores with a potential on grid
Definition qs_integrate_potential_single.F:486
qs_integrate_potential_single::integrate_v_rspace_one_center
subroutine, public integrate_v_rspace_one_center(v_rspace, qs_env, int_res, calculate_forces, basis_type, atomlist)
computes integrals of product of v_rspace times a one-center function required for LRIGPW
Definition qs_integrate_potential_single.F:766
qs_integrate_potential_single::integrate_ppl_rspace
subroutine, public integrate_ppl_rspace(rho_rspace, qs_env)
computes the forces/virial due to the local pseudopotential
Definition qs_integrate_potential_single.F:101
qs_integrate_potential_single::integrate_rho_nlcc
subroutine, public integrate_rho_nlcc(rho_rspace, qs_env)
computes the forces/virial due to the nlcc pseudopotential
Definition qs_integrate_potential_single.F:286
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::transfer_pw2rs
subroutine, public transfer_pw2rs(rs, pw)
...
Definition realspace_grid_types.F:694
realspace_grid_types::map_gaussian_here
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
Definition realspace_grid_types.F:2144
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1905
rs_pw_interface
Transfers densities from PW to RS grids and potentials from PW to RS.
Definition rs_pw_interface.F:14
rs_pw_interface::potential_pw2rs
subroutine, public potential_pw2rs(rs_v, v_rspace, pw_env)
transfers a potential from a pw_grid to a vector of realspace multigrids
Definition rs_pw_interface.F:152
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
gaussian_gridlevels::gridlevel_info_type
Definition gaussian_gridlevels.F:34
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:332
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136
virial_types::virial_type
Definition virial_types.F:26