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