(git:15c1bfc)
Loading...
Searching...
No Matches
qs_rho0_ggrid.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! **************************************************************************************************
13 USE cell_types, ONLY: cell_type,&
14 pbc
17 USE grid_api, ONLY: grid_func_ab,&
19 USE kinds, ONLY: dp
22 USE orbital_pointers, ONLY: indco,&
23 nco,&
24 ncoset,&
25 nso,&
26 nsoset
29 USE pw_env_types, ONLY: pw_env_get,&
31 USE pw_methods, ONLY: pw_axpy,&
32 pw_copy,&
34 pw_scale,&
37 USE pw_pool_types, ONLY: pw_pool_p_type,&
39 USE pw_types, ONLY: pw_c1d_gs_type,&
51 USE qs_kind_types, ONLY: get_qs_kind,&
55 USE qs_rho0_types, ONLY: get_rho0_mpole,&
68 USE util, ONLY: get_limit
69 USE virial_types, ONLY: virial_type
70#include "./base/base_uses.f90"
71
72 IMPLICIT NONE
73
74 PRIVATE
75
76 ! Global parameters (only in this module)
77
78 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_ggrid'
79
80 ! Public subroutines
81
83
84CONTAINS
85
86! **************************************************************************************************
87!> \brief ...
88!> \param qs_env ...
89!> \param rho0 ...
90!> \param tot_rs_int ...
91!> \param my_pools ...
92!> \param my_rs_grids ...
93!> \param my_rs_descs ...
94! **************************************************************************************************
95 SUBROUTINE put_rho0_on_grid(qs_env, rho0, tot_rs_int, my_pools, my_rs_grids, my_rs_descs)
96
97 TYPE(qs_environment_type), POINTER :: qs_env
98 TYPE(rho0_mpole_type), POINTER :: rho0
99 REAL(kind=dp), INTENT(OUT) :: tot_rs_int
100 TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
101 POINTER :: my_pools
102 TYPE(realspace_grid_type), DIMENSION(:), &
103 OPTIONAL, POINTER :: my_rs_grids
104 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
105 OPTIONAL, POINTER :: my_rs_descs
106
107 CHARACTER(LEN=*), PARAMETER :: routinen = 'put_rho0_on_grid'
108
109 INTEGER :: auxbas_grid, handle, iat, iatom, igrid, &
110 ikind, ithread, j, l0_ikind, lmax0, &
111 nat, nch_ik, nch_max, npme
112 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
113 LOGICAL :: paw_atom
114 REAL(kind=dp) :: eps_rho_rspace, rpgf0, zet0
115 REAL(kind=dp), DIMENSION(3) :: ra
116 REAL(kind=dp), DIMENSION(:), POINTER :: qlm_c
117 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
118 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
119 TYPE(cell_type), POINTER :: cell
120 TYPE(dft_control_type), POINTER :: dft_control
121 TYPE(mp_para_env_type), POINTER :: para_env
122 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
123 TYPE(pw_c1d_gs_type) :: coeff_gspace
124 TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs
125 TYPE(pw_env_type), POINTER :: pw_env
126 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
127 TYPE(pw_pool_type), POINTER :: pw_pool
128 TYPE(pw_r3d_rs_type) :: coeff_rspace, rho0_r_tmp
129 TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs
130 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
131 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
132 POINTER :: descs
133 TYPE(realspace_grid_desc_type), POINTER :: desc
134 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: grids
135 TYPE(realspace_grid_type), POINTER :: rs_grid
136
137 CALL timeset(routinen, handle)
138
139 NULLIFY (atomic_kind_set, qs_kind_set, cores, pab, qlm_c)
140
141 NULLIFY (dft_control, pw_env, particle_set, para_env, cell, rho0_s_gs, rho0_s_rs)
142 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
143 particle_set=particle_set, &
144 atomic_kind_set=atomic_kind_set, &
145 qs_kind_set=qs_kind_set, &
146 para_env=para_env, &
147 pw_env=pw_env, cell=cell)
148 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
149
150 NULLIFY (descs, pw_pools)
151 CALL pw_env_get(pw_env=pw_env, rs_descs=descs, rs_grids=grids, pw_pools=pw_pools)
152 auxbas_grid = pw_env%auxbas_grid
153
154 NULLIFY (rho0_s_gs, rho0_s_rs)
155 CALL get_rho0_mpole(rho0_mpole=rho0, lmax_0=lmax0, &
156 zet0_h=zet0, igrid_zet0_s=igrid, &
157 rho0_s_gs=rho0_s_gs, &
158 rho0_s_rs=rho0_s_rs)
159
160 ! *** set up the rs grid at level igrid
161 NULLIFY (rs_grid, desc, pw_pool)
162 ! IF present, overwrite qs grid for new pool
163 IF (PRESENT(my_pools)) THEN
164 desc => my_rs_descs(igrid)%rs_desc
165 rs_grid => my_rs_grids(igrid)
166 pw_pool => my_pools(igrid)%pool
167 ELSE
168 desc => descs(igrid)%rs_desc
169 rs_grid => grids(igrid)
170 pw_pool => pw_pools(igrid)%pool
171 END IF
172
173 cpassert(ASSOCIATED(desc))
174 cpassert(ASSOCIATED(pw_pool))
175
176 IF (igrid /= auxbas_grid) THEN
177 CALL pw_pool%create_pw(coeff_rspace)
178 CALL pw_pool%create_pw(coeff_gspace)
179 END IF
180 CALL rs_grid_zero(rs_grid)
181
182 nch_max = ncoset(lmax0)
183
184 ALLOCATE (pab(nch_max, 1))
185
186 DO ikind = 1, SIZE(atomic_kind_set)
187 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
188 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
189
190 IF (.NOT. paw_atom .AND. dft_control%qs_control%gapw_control%nopaw_as_gpw) cycle
191
192 CALL get_rho0_mpole(rho0_mpole=rho0, ikind=ikind, l0_ikind=l0_ikind, &
193 rpgf0_s=rpgf0)
194
195 nch_ik = ncoset(l0_ikind)
196 pab = 0.0_dp
197
198 CALL reallocate(cores, 1, nat)
199 npme = 0
200 cores = 0
201
202 DO iat = 1, nat
203 iatom = atom_list(iat)
204 ra(:) = pbc(particle_set(iatom)%r, cell)
205 IF (rs_grid%desc%parallel .AND. .NOT. rs_grid%desc%distributed) THEN
206 ! replicated realspace grid, split the atoms up between procs
207 IF (modulo(nat, rs_grid%desc%group_size) == rs_grid%desc%my_pos) THEN
208 npme = npme + 1
209 cores(npme) = iat
210 END IF
211 ELSE
212 npme = npme + 1
213 cores(npme) = iat
214 END IF
215
216 END DO
217
218 ithread = 0
219 DO j = 1, npme
220
221 iat = cores(j)
222 iatom = atom_list(iat)
223
224 CALL get_rho0_mpole(rho0_mpole=rho0, iat=iatom, qlm_car=qlm_c)
225
226 pab(1:nch_ik, 1) = qlm_c(1:nch_ik)
227
228 ra(:) = pbc(particle_set(iatom)%r, cell)
229
231 l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
232 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
233 rs_grid, ga_gb_function=grid_func_ab, radius=rpgf0, &
234 use_subpatch=.true., subpatch_pattern=0)
235
236 END DO ! j
237
238 END DO ! ikind
239
240 IF (ASSOCIATED(cores)) THEN
241 DEALLOCATE (cores)
242 END IF
243
244 DEALLOCATE (pab)
245
246 IF (igrid /= auxbas_grid) THEN
247 CALL transfer_rs2pw(rs_grid, coeff_rspace)
248 CALL pw_zero(rho0_s_gs)
249 CALL pw_transfer(coeff_rspace, coeff_gspace)
250 CALL pw_axpy(coeff_gspace, rho0_s_gs)
251
252 tot_rs_int = pw_integrate_function(coeff_rspace, isign=-1)
253
254 CALL pw_pool%give_back_pw(coeff_rspace)
255 CALL pw_pool%give_back_pw(coeff_gspace)
256 ELSE
257
258 CALL pw_pool%create_pw(rho0_r_tmp)
259
260 CALL transfer_rs2pw(rs_grid, rho0_r_tmp)
261
262 tot_rs_int = pw_integrate_function(rho0_r_tmp, isign=-1)
263
264 CALL pw_transfer(rho0_r_tmp, rho0_s_rs)
265 CALL pw_pool%give_back_pw(rho0_r_tmp)
266
267 CALL pw_zero(rho0_s_gs)
268 CALL pw_transfer(rho0_s_rs, rho0_s_gs)
269 END IF
270 CALL timestop(handle)
271
272 END SUBROUTINE put_rho0_on_grid
273
274! **************************************************************************************************
275!> \brief ...
276!> \param pw_env ...
277!> \param rho0_mpole ...
278! **************************************************************************************************
279 SUBROUTINE rho0_s_grid_create(pw_env, rho0_mpole)
280
281 TYPE(pw_env_type), POINTER :: pw_env
282 TYPE(rho0_mpole_type), POINTER :: rho0_mpole
283
284 CHARACTER(len=*), PARAMETER :: routinen = 'rho0_s_grid_create'
285
286 INTEGER :: handle
287 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
288
289 CALL timeset(routinen, handle)
290
291 cpassert(ASSOCIATED(pw_env))
292
293 NULLIFY (auxbas_pw_pool)
294 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
295 cpassert(ASSOCIATED(auxbas_pw_pool))
296
297 ! reallocate rho0 on the global grid in real and reciprocal space
298 cpassert(ASSOCIATED(rho0_mpole))
299
300 ! rho0 density in real space
301 IF (ASSOCIATED(rho0_mpole%rho0_s_rs)) THEN
302 CALL rho0_mpole%rho0_s_rs%release()
303 ELSE
304 ALLOCATE (rho0_mpole%rho0_s_rs)
305 END IF
306 CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_rs)
307
308 ! rho0 density in reciprocal space
309 IF (ASSOCIATED(rho0_mpole%rho0_s_gs)) THEN
310 CALL rho0_mpole%rho0_s_gs%release()
311 ELSE
312 ALLOCATE (rho0_mpole%rho0_s_gs)
313 END IF
314 CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_gs)
315
316 ! Find the grid level suitable for rho0_soft
317 rho0_mpole%igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*rho0_mpole%zet0_h)
318
319 CALL timestop(handle)
320
321 IF (rho0_mpole%do_cneo) THEN
322 CALL rhoz_cneo_s_grid_create(pw_env, rho0_mpole)
323 END IF
324
325 END SUBROUTINE rho0_s_grid_create
326
327! **************************************************************************************************
328!> \brief ...
329!> \param qs_env ...
330!> \param v_rspace ...
331!> \param para_env ...
332!> \param calculate_forces ...
333!> \param local_rho_set ...
334!> \param local_rho_set_2nd ...
335!> \param atener ...
336!> \param kforce ...
337!> \param my_pools ...
338!> \param my_rs_descs ...
339! **************************************************************************************************
340 SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
341 local_rho_set_2nd, atener, kforce, my_pools, my_rs_descs)
342
343 TYPE(qs_environment_type), POINTER :: qs_env
344 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
345 TYPE(mp_para_env_type), POINTER :: para_env
346 LOGICAL, INTENT(IN) :: calculate_forces
347 TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set, local_rho_set_2nd
348 REAL(kind=dp), DIMENSION(:), OPTIONAL :: atener
349 REAL(kind=dp), INTENT(IN), OPTIONAL :: kforce
350 TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
351 POINTER :: my_pools
352 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
353 OPTIONAL, POINTER :: my_rs_descs
354
355 CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_vhg0_rspace'
356
357 INTEGER :: auxbas_grid, bo(2), handle, i, iat, iatom, ic, icg, ico, ig1, ig2, igrid, ii, &
358 ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, &
359 llmax_nuc, lmax0, lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, &
360 max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, n2, nat, nch_ik, nch_max, &
361 ncurr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
362 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
363 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
364 INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmax_nuc, lmin, &
365 lmin_nuc, npgf, npgf_nuc
366 LOGICAL :: grid_distributed, paw_atom, use_virial
367 REAL(kind=dp) :: eps_rho_rspace, force_tmp(3), fscale, &
368 ra(3), rpgf0, zet0
369 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
370 REAL(kind=dp), DIMENSION(:), POINTER :: hab_sph, norm_l, qlm
371 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, intloc_nuc, pab
372 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: a_hdab_sph, hdab, qlm_gg
373 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: a_hdab
374 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
375 TYPE(cell_type), POINTER :: cell
376 TYPE(cneo_potential_type), POINTER :: cneo_potential
377 TYPE(dft_control_type), POINTER :: dft_control
378 TYPE(gto_basis_set_type), POINTER :: basis_1c_set, nuc_basis_set
379 TYPE(harmonics_atom_type), POINTER :: harmonics
380 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
381 TYPE(pw_c1d_gs_type) :: coeff_gaux, coeff_gspace
382 TYPE(pw_env_type), POINTER :: pw_env
383 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
384 TYPE(pw_pool_type), POINTER :: pw_aux, pw_pool
385 TYPE(pw_r3d_rs_type) :: coeff_raux, coeff_rspace
386 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
387 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
388 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
389 POINTER :: rs_descs
390 TYPE(realspace_grid_desc_type), POINTER :: rs_desc
391 TYPE(realspace_grid_type) :: rs_v
392 TYPE(rho0_mpole_type), POINTER :: rho0_mpole
393 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
394 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
395 TYPE(rho_atom_type), POINTER :: rho_atom
396 TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
397 TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
398 TYPE(virial_type), POINTER :: virial
399
400 CALL timeset(routinen, handle)
401
402 ! In case of the external density computed forces probably also
403 ! need to be stored outside qs_env. We can then remove the
404 ! attribute 'OPTIONAL' from the argument 'local_rho_set'.
405
406 ! CPASSERT(.NOT. (calculate_forces .AND. PRESENT(local_rho_set)))
407! IF (calculate_forces .AND. PRESENT(local_rho_set)) THEN
408! CPWARN("Forces and External Density!")
409! END IF
410
411 NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
412 NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set, rhoz_cneo_set)
413
414 CALL get_qs_env(qs_env=qs_env, &
415 atomic_kind_set=atomic_kind_set, &
416 qs_kind_set=qs_kind_set, &
417 cell=cell, &
418 dft_control=dft_control, &
419 force=force, pw_env=pw_env, &
420 rho0_mpole=rho0_mpole, &
421 rho_atom_set=rho_atom_set, &
422 rhoz_cneo_set=rhoz_cneo_set, &
423 particle_set=particle_set, &
424 virial=virial)
425
426 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
427
428 nspins = dft_control%nspins
429
430 ! The aim of the following code was to return immediately if the subroutine
431 ! was called for triplet excited states in spin-restricted case. This check
432 ! is also performed before invocation of this subroutine. It should be save
433 ! to remove the optional argument 'do_triplet' from the subroutine interface.
434 !my_tddft = PRESENT(local_rho_set)
435 !IF (my_tddft) THEN
436 ! IF (PRESENT(do_triplet)) THEN
437 ! IF (nspins == 1 .AND. do_triplet) RETURN
438 ! ELSE
439 ! IF (nspins == 1 .AND. dft_control%tddfpt_control%res_etype /= tddfpt_singlet) RETURN
440 ! END IF
441 !END IF
442
443 IF (PRESENT(local_rho_set)) &
444 CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set, &
445 rhoz_cneo_set=rhoz_cneo_set)
446 ! Q from rho0_mpole of local_rho_set
447 ! for TDDFT forces we need mixed potential / integral space
448 ! potential stored on local_rho_set_2nd
449 IF (PRESENT(local_rho_set_2nd)) THEN
450 CALL get_local_rho(local_rho_set_2nd, rho_atom_set=rho_atom_set)
451 END IF
452 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax0, &
453 zet0_h=zet0, igrid_zet0_s=igrid, &
454 norm_g0l_h=norm_l)
455
456 ! Setup of the potential on the multigrids
457 NULLIFY (rs_descs, pw_pools)
458 cpassert(ASSOCIATED(pw_env))
459 CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools)
460
461 ! Assign from pw_env
462 auxbas_grid = pw_env%auxbas_grid
463
464 ! IF present, overwrite qs grid for new pool
465 ! Get the potential on the right grid
466 IF (PRESENT(my_pools)) THEN
467 rs_desc => my_rs_descs(igrid)%rs_desc
468 pw_pool => my_pools(igrid)%pool
469 ELSE
470 rs_desc => rs_descs(igrid)%rs_desc
471 pw_pool => pw_pools(igrid)%pool
472 END IF
473
474 CALL pw_pool%create_pw(coeff_gspace)
475 CALL pw_pool%create_pw(coeff_rspace)
476
477 IF (igrid /= auxbas_grid) THEN
478 pw_aux => pw_pools(auxbas_grid)%pool
479 CALL pw_aux%create_pw(coeff_gaux)
480 CALL pw_transfer(v_rspace, coeff_gaux)
481 CALL pw_copy(coeff_gaux, coeff_gspace)
482 CALL pw_transfer(coeff_gspace, coeff_rspace)
483 CALL pw_aux%give_back_pw(coeff_gaux)
484 CALL pw_aux%create_pw(coeff_raux)
485 fscale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol
486 CALL pw_scale(coeff_rspace, fscale)
487 CALL pw_aux%give_back_pw(coeff_raux)
488 ELSE
489
490 IF (coeff_gspace%pw_grid%spherical) THEN
491 CALL pw_transfer(v_rspace, coeff_gspace)
492 CALL pw_transfer(coeff_gspace, coeff_rspace)
493 ELSE
494 CALL pw_copy(v_rspace, coeff_rspace)
495 END IF
496 END IF
497 CALL pw_pool%give_back_pw(coeff_gspace)
498
499 ! Setup the rs grid at level igrid
500 CALL rs_grid_create(rs_v, rs_desc)
501 CALL rs_grid_zero(rs_v)
502 CALL transfer_pw2rs(rs_v, coeff_rspace)
503
504 CALL pw_pool%give_back_pw(coeff_rspace)
505
506 ! Now the potential is on the right grid => integration
507
508 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
509
510 ! Allocate work storage
511
512 NULLIFY (hab, hab_sph, hdab, hdab_sph, pab, a_hdab, a_hdab_sph)
513 nch_max = ncoset(lmax0)
514 CALL reallocate(hab, 1, nch_max, 1, 1)
515 CALL reallocate(hab_sph, 1, nch_max)
516 CALL reallocate(hdab, 1, 3, 1, nch_max, 1, 1)
517 CALL reallocate(hdab_sph, 1, 3, 1, nch_max)
518 CALL reallocate(a_hdab, 1, 3, 1, 3, 1, nch_max, 1, 1)
519 CALL reallocate(a_hdab_sph, 1, 3, 1, 3, 1, nch_max)
520 CALL reallocate(pab, 1, nch_max, 1, 1)
521
522 ncurr = -1
523
524 grid_distributed = rs_v%desc%distributed
525
526 fscale = 1.0_dp
527 IF (PRESENT(kforce)) THEN
528 fscale = kforce
529 END IF
530
531 DO ikind = 1, SIZE(atomic_kind_set, 1)
532 NULLIFY (basis_1c_set, atom_list, harmonics, cneo_potential)
533 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
534 CALL get_qs_kind(qs_kind_set(ikind), &
535 basis_set=basis_1c_set, basis_type="GAPW_1C", &
536 paw_atom=paw_atom, &
537 harmonics=harmonics, &
538 cneo_potential=cneo_potential)
539
540 IF (.NOT. paw_atom) cycle
541
542 NULLIFY (qlm_gg, lmax, npgf)
543 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
544 l0_ikind=l0_ikind, qlm_gg=qlm_gg, & ! Qs different
545 rpgf0_s=rpgf0)
546
547 CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
548 lmax=lmax, lmin=lmin, &
549 maxso=maxso, maxl=maxl, &
550 nset=nset, npgf=npgf)
551
552 nsotot = maxso*nset
553 ALLOCATE (intloc(nsotot, nsotot))
554
555 ! Initialize the local KS integrals
556
557 nch_ik = ncoset(l0_ikind)
558 pab = 1.0_dp
559 max_s_harm = harmonics%max_s_harm
560 llmax = harmonics%llmax
561
562 NULLIFY (intloc_nuc)
563 maxl_nuc = -1
564 max_s_harm_nuc = 0
565 llmax_nuc = -1
566 IF (ASSOCIATED(cneo_potential)) THEN
567 NULLIFY (nuc_basis_set)
568 CALL get_qs_kind(qs_kind_set(ikind), &
569 basis_set=nuc_basis_set, &
570 basis_type="NUC")
571
572 CALL get_gto_basis_set(gto_basis_set=nuc_basis_set, &
573 lmax=lmax_nuc, lmin=lmin_nuc, &
574 maxso=maxso_nuc, maxl=maxl_nuc, &
575 nset=nset_nuc, npgf=npgf_nuc)
576 nsotot_nuc = maxso_nuc*nset_nuc
577 ALLOCATE (intloc_nuc(nsotot_nuc, nsotot_nuc))
578 max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
579 llmax_nuc = cneo_potential%harmonics%llmax
580 END IF
581
582 ALLOCATE (cg_list(2, nsoset(max(maxl, maxl_nuc))**2, &
583 max(max_s_harm, max_s_harm_nuc)), &
584 cg_n_list(max(max_s_harm, max_s_harm_nuc)))
585
586 num_pe = para_env%num_pe
587 mepos = para_env%mepos
588 DO j = 0, num_pe - 1
589 bo = get_limit(nat, num_pe, j)
590 IF (.NOT. grid_distributed .AND. j /= mepos) cycle
591
592 DO iat = bo(1), bo(2)
593 iatom = atom_list(iat)
594 ra(:) = pbc(particle_set(iatom)%r, cell)
595
596 NULLIFY (qlm)
597 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, iat=iatom, qlm_tot=qlm)
598
599 hab = 0.0_dp
600 hdab = 0.0_dp
601 intloc = 0._dp
602 IF (use_virial) THEN
603 my_virial_a = 0.0_dp
604 my_virial_b = 0.0_dp
605 a_hdab = 0.0_dp
606 END IF
607
609 l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
610 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), rs_v, &
611 hab, pab, o1=0, o2=0, &
612 radius=rpgf0, &
613 calculate_forces=calculate_forces, &
614 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, &
615 hdab=hdab, a_hdab=a_hdab, use_subpatch=.true., subpatch_pattern=0)
616
617 ! Convert from cartesian to spherical
618 DO lshell = 0, l0_ikind
619 DO is = 1, nso(lshell)
620 iso = is + nsoset(lshell - 1)
621 hab_sph(iso) = 0.0_dp
622 hdab_sph(1:3, iso) = 0.0_dp
623 a_hdab_sph(1:3, 1:3, iso) = 0.0_dp
624 DO ic = 1, nco(lshell)
625 ico = ic + ncoset(lshell - 1)
626 lx = indco(1, ico)
627 ly = indco(2, ico)
628 lz = indco(3, ico)
629 hab_sph(iso) = hab_sph(iso) + &
630 norm_l(lshell)* &
631 orbtramat(lshell)%slm(is, ic)* &
632 hab(ico, 1)
633 IF (calculate_forces) THEN
634 hdab_sph(1:3, iso) = hdab_sph(1:3, iso) + &
635 norm_l(lshell)* &
636 orbtramat(lshell)%slm(is, ic)* &
637 hdab(1:3, ico, 1)
638 END IF
639 IF (use_virial) THEN
640 DO ii = 1, 3
641 DO i = 1, 3
642 a_hdab_sph(i, ii, iso) = a_hdab_sph(i, ii, iso) + &
643 norm_l(lshell)* &
644 orbtramat(lshell)%slm(is, ic)* &
645 a_hdab(i, ii, ico, 1)
646 END DO
647 END DO
648 END IF
649
650 END DO ! ic
651 END DO ! is
652 END DO ! lshell
653
654 m1 = 0
655 DO iset1 = 1, nset
656
657 m2 = 0
658 DO iset2 = 1, nset
659 CALL get_none0_cg_list(harmonics%my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
660 max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
661 n1 = nsoset(lmax(iset1))
662 DO ipgf1 = 1, npgf(iset1)
663 n2 = nsoset(lmax(iset2))
664 DO ipgf2 = 1, npgf(iset2)
665
666 DO iso = 1, min(nsoset(l0_ikind), max_iso_not0_local)
667 DO icg = 1, cg_n_list(iso)
668 iso1 = cg_list(1, icg, iso)
669 iso2 = cg_list(2, icg, iso)
670
671 ig1 = iso1 + n1*(ipgf1 - 1) + m1
672 ig2 = iso2 + n2*(ipgf2 - 1) + m2
673
674 intloc(ig1, ig2) = intloc(ig1, ig2) + qlm_gg(ig1, ig2, iso)*hab_sph(iso) ! potential times Q
675
676 END DO ! icg
677 END DO ! iso
678
679 END DO ! ipgf2
680 END DO ! ipgf1
681 m2 = m2 + maxso
682 END DO ! iset2
683 m1 = m1 + maxso
684 END DO ! iset1
685
686 IF (ASSOCIATED(cneo_potential)) THEN
687 intloc_nuc = 0.0_dp
688 m1 = 0
689 DO iset1 = 1, nset_nuc
690 n1 = nsoset(lmax_nuc(iset1))
691 m2 = 0
692 DO iset2 = 1, nset_nuc
693 n2 = nsoset(lmax_nuc(iset2))
694 CALL get_none0_cg_list(cneo_potential%harmonics%my_CG, lmin_nuc(iset1), &
695 lmax_nuc(iset1), lmin_nuc(iset2), lmax_nuc(iset2), &
696 max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
697 max_iso_not0_local)
698 DO ipgf1 = 1, npgf_nuc(iset1)
699 DO ipgf2 = 1, npgf_nuc(iset2)
700
701 DO iso = 1, min(nsoset(l0_ikind), max_iso_not0_local)
702 DO icg = 1, cg_n_list(iso)
703 iso1 = cg_list(1, icg, iso)
704 iso2 = cg_list(2, icg, iso)
705
706 ig1 = iso1 + n1*(ipgf1 - 1) + m1
707 ig2 = iso2 + n2*(ipgf2 - 1) + m2
708
709 intloc_nuc(ig1, ig2) = intloc_nuc(ig1, ig2) - cneo_potential%zeff* &
710 cneo_potential%Qlm_gg(ig1, ig2, iso)*hab_sph(iso)
711
712 END DO ! icg
713 END DO ! iso
714
715 END DO ! ipgf2
716 END DO ! ipgf1
717 m2 = m2 + maxso_nuc
718 END DO ! iset2
719 m1 = m1 + maxso_nuc
720 END DO ! iset1
721 END IF
722
723 IF (grid_distributed) THEN
724 ! Sum result over all processors
725 CALL para_env%sum(intloc)
726 IF (ASSOCIATED(cneo_potential)) THEN
727 CALL para_env%sum(intloc_nuc)
728 END IF
729 END IF
730
731 IF (j == mepos) THEN
732 rho_atom => rho_atom_set(iatom)
733 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_local_h, ga_vlocal_gb_s=int_local_s)
734 DO ispin = 1, nspins
735 int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc
736 int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc
737 END DO
738 IF (ASSOCIATED(cneo_potential)) THEN
739 rhoz_cneo => rhoz_cneo_set(iatom)
740 rhoz_cneo%ga_Vlocal_gb_h = rhoz_cneo%ga_Vlocal_gb_h + intloc_nuc
741 rhoz_cneo%ga_Vlocal_gb_s = rhoz_cneo%ga_Vlocal_gb_s + intloc_nuc
742 END IF
743 END IF
744
745 IF (PRESENT(atener)) THEN
746 DO iso = 1, nsoset(l0_ikind)
747 atener(iatom) = atener(iatom) + 0.5_dp*qlm(iso)*hab_sph(iso)
748 END DO
749 END IF
750
751 IF (calculate_forces) THEN
752 force_tmp(1:3) = 0.0_dp
753 DO iso = 1, nsoset(l0_ikind)
754 force_tmp(1) = force_tmp(1) + qlm(iso)*hdab_sph(1, iso) ! Q here is from local_rho_set
755 force_tmp(2) = force_tmp(2) + qlm(iso)*hdab_sph(2, iso)
756 force_tmp(3) = force_tmp(3) + qlm(iso)*hdab_sph(3, iso)
757 END DO
758 force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + fscale*force_tmp(1:3)
759 END IF
760 IF (use_virial) THEN
761 my_virial_a = 0.0_dp
762 DO iso = 1, nsoset(l0_ikind)
763 DO ii = 1, 3
764 DO i = 1, 3
765 ! Q from local_rho_set
766 virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*qlm(iso)*a_hdab_sph(i, ii, iso)
767 virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + fscale*qlm(iso)*a_hdab_sph(i, ii, iso)
768 END DO
769 END DO
770 END DO
771 END IF
772
773 END DO
774 END DO
775
776 DEALLOCATE (intloc)
777 IF (ASSOCIATED(intloc_nuc)) DEALLOCATE (intloc_nuc)
778 DEALLOCATE (cg_list, cg_n_list)
779
780 END DO ! ikind
781
782 CALL rs_grid_release(rs_v)
783
784 DEALLOCATE (hab, hdab, hab_sph, hdab_sph, pab, a_hdab, a_hdab_sph)
785
786 CALL timestop(handle)
787
788 IF (rho0_mpole%do_cneo) THEN
789 IF (PRESENT(kforce)) THEN
790 CALL integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, &
791 rhoz_cneo_set, kforce)
792 ELSE
793 CALL integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, &
794 rhoz_cneo_set)
795 END IF
796 END IF
797
798 END SUBROUTINE integrate_vhg0_rspace
799
800END MODULE qs_rho0_ggrid
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Define the atomic kind types and their sub types.
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.
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)
...
Handles all functions related to the CELL.
Definition cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
integer, parameter, public grid_func_ab
Definition grid_api.F:29
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
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
Definition grid_api.F:119
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Utility routines for the memory handling.
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nco
integer, dimension(:), allocatable, public nsoset
integer, dimension(:), allocatable, public ncoset
integer, dimension(:, :), allocatable, public indco
integer, dimension(:), allocatable, public nso
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
type(orbtramat_type), dimension(:), pointer, public orbtramat
Define the data structure for the particle information.
container for various plainwaves related things
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
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
CNEO soft nuclear densities on the global grid (see J. Chem. Theory Comput. 2025, 21,...
subroutine, public integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, rhoz_cneo_set, kforce)
...
subroutine, public rhoz_cneo_s_grid_create(pw_env, rho0_mpole)
...
Types used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
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, 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, 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, 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.
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
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, 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.
subroutine, public get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set)
...
subroutine, public put_rho0_on_grid(qs_env, rho0, tot_rs_int, my_pools, my_rs_grids, my_rs_descs)
...
subroutine, public rho0_s_grid_create(pw_env, rho0_mpole)
...
subroutine, public integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, local_rho_set_2nd, atener, kforce, my_pools, my_rs_descs)
...
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, qlm_gg, qlm_car, qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs)
...
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
subroutine, public rs_grid_create(rs, desc)
...
subroutine, public transfer_pw2rs(rs, pw)
...
subroutine, public transfer_rs2pw(rs, pw)
...
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
All kind of helpful little routines.
Definition util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
stores all the informations relevant to an mpi environment
contained for different pw related things
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.