(git:15c1bfc)
Loading...
Searching...
No Matches
qs_cneo_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
8! **************************************************************************************************
9!> \brief CNEO soft nuclear densities on the global grid
10!> (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
11!> \par History
12!> 08.2025 created [zc62]
13!> \author Zehua Chen
14! **************************************************************************************************
15MODULE qs_cneo_ggrid
16 USE ao_util, ONLY: exp_radius_very_extended
17 USE atomic_kind_types, ONLY: atomic_kind_type,&
18 get_atomic_kind
19 USE basis_set_types, ONLY: get_gto_basis_set,&
20 gto_basis_set_type
21 USE cell_types, ONLY: cell_type,&
22 pbc
23 USE cp_control_types, ONLY: dft_control_type
24 USE gaussian_gridlevels, ONLY: gaussian_gridlevel,&
25 gridlevel_info_type
26 USE grid_api, ONLY: grid_func_ab,&
27 collocate_pgf_product
28 USE kinds, ONLY: dp
29 USE message_passing, ONLY: mp_para_env_type
30 USE orbital_pointers, ONLY: ncoset
31 USE particle_types, ONLY: particle_type
32 USE pw_env_types, ONLY: pw_env_get,&
33 pw_env_type
34 USE pw_methods, ONLY: pw_axpy,&
35 pw_integrate_function,&
36 pw_transfer,&
37 pw_zero
38 USE pw_pool_types, ONLY: pw_pool_p_type,&
39 pw_pool_type,&
40 pw_pools_create_pws,&
41 pw_pools_give_back_pws
42 USE pw_types, ONLY: pw_c1d_gs_type,&
43 pw_r3d_rs_type
44 USE qs_cneo_types, ONLY: cneo_potential_type,&
45 get_cneo_potential,&
46 rhoz_cneo_type
47 USE qs_environment_types, ONLY: get_qs_env,&
48 qs_environment_type
49 USE qs_force_types, ONLY: qs_force_type
50 USE qs_integrate_potential, ONLY: integrate_pgf_product
51 USE qs_kind_types, ONLY: get_qs_kind,&
52 get_qs_kind_set,&
53 qs_kind_type
54 USE qs_rho0_types, ONLY: get_rho0_mpole,&
55 rho0_mpole_type
56 USE realspace_grid_types, ONLY: map_gaussian_here,&
57 realspace_grid_type,&
58 rs_grid_zero,&
59 transfer_rs2pw
60 USE rs_pw_interface, ONLY: potential_pw2rs
61 USE virial_types, ONLY: virial_type
62#include "./base/base_uses.f90"
63
64 IMPLICIT NONE
65
66 PRIVATE
67
68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_ggrid'
69
70 PUBLIC :: put_rhoz_cneo_s_on_grid, rhoz_cneo_s_grid_create, integrate_vhgg_rspace
71
72CONTAINS
73
74! **************************************************************************************************
75!> \brief ...
76!> \param qs_env ...
77!> \param rho0 ...
78!> \param rhoz_cneo_set ...
79!> \param tot_rs_int ...
80! **************************************************************************************************
81 SUBROUTINE put_rhoz_cneo_s_on_grid(qs_env, rho0, rhoz_cneo_set, tot_rs_int)
82
83 TYPE(qs_environment_type), POINTER :: qs_env
84 TYPE(rho0_mpole_type), POINTER :: rho0
85 TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
86 REAL(kind=dp), INTENT(OUT) :: tot_rs_int
87
88 CHARACTER(LEN=*), PARAMETER :: routinen = 'put_rhoz_cneo_s_on_grid'
89
90 INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
91 m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, npgf2, nseta, offset, sgfa, &
92 sgfb
93 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
94 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
95 LOGICAL :: paw_atom
96 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
97 REAL(kind=dp) :: eps_rho_rspace, radius, scale, zeff, zetp
98 REAL(kind=dp), DIMENSION(3) :: ra
99 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, work, zeta
100 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
101 TYPE(cell_type), POINTER :: cell
102 TYPE(cneo_potential_type), POINTER :: cneo_potential
103 TYPE(dft_control_type), POINTER :: dft_control
104 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
105 TYPE(gto_basis_set_type), POINTER :: nuc_soft_basis
106 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
107 TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
108 TYPE(pw_c1d_gs_type), POINTER :: rhoz_cneo_s_gs
109 TYPE(pw_env_type), POINTER :: pw_env
110 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
111 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
112 TYPE(pw_r3d_rs_type), POINTER :: rhoz_cneo_s_rs
113 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
114 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
115 TYPE(realspace_grid_type), POINTER :: rs_grid
116
117 CALL timeset(routinen, handle)
118
119 NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
120 first_sgfa, gridlevel_info, la_max, la_min, nuc_soft_basis, &
121 npgfa, nsgfa, p_block, pab, particle_set, pw_env, pw_pools, &
122 rs_grid, rs_rho, sphi_a, work, zeta)
123
124 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
125 atomic_kind_set=atomic_kind_set, &
126 cell=cell, particle_set=particle_set, &
127 pw_env=pw_env, &
128 dft_control=dft_control)
129
130 ! find maximum numbers
131 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
132 maxco=maxco, &
133 maxsgf_set=maxsgf_set, &
134 basis_type="NUC_SOFT")
135 IF (maxco > 0) THEN
136 ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set))
137
138 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
139
140 ! *** set up the pw multi-grids *** !
141 cpassert(ASSOCIATED(pw_env))
142 CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
143 gridlevel_info=gridlevel_info)
144
145 CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
146
147 CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
148
149 ! *** set up the rs multi-grids *** !
150 DO igrid_level = 1, gridlevel_info%ngrid_levels
151 CALL rs_grid_zero(rs_rho(igrid_level))
152 END DO
153
154 offset = 0
155 my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
156 group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
157
158 DO ikind = 1, SIZE(atomic_kind_set)
159 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
160 IF (.NOT. paw_atom) cycle
161
162 NULLIFY (cneo_potential)
163 CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
164 IF (.NOT. ASSOCIATED(cneo_potential)) cycle
165
166 NULLIFY (atom_list)
167 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
168 NULLIFY (nuc_soft_basis)
169 CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_soft_basis, basis_type="NUC_SOFT")
170 CALL get_gto_basis_set(gto_basis_set=nuc_soft_basis, lmax=la_max, &
171 lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
172 sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
173 CALL get_cneo_potential(cneo_potential, zeff=zeff)
174 m1 = maxval(npgfa(1:nseta))
175 ALLOCATE (map_it2(m1, m1))
176
177 DO iatom = 1, natom
178 atom_a = atom_list(iatom)
179 IF (rhoz_cneo_set(atom_a)%ready) THEN
180 ra(:) = pbc(particle_set(atom_a)%r, cell)
181 p_block => rhoz_cneo_set(atom_a)%pmat
182 DO iset = 1, nseta
183 DO jset = 1, iset
184 ! processor mapping
185 map_it2 = .false.
186 DO ipgf = 1, npgfa(iset)
187 IF (jset == iset) THEN
188 npgf2 = ipgf
189 ELSE
190 npgf2 = npgfa(jset)
191 END IF
192 DO jpgf = 1, npgf2
193 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
194 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
195 rs_grid => rs_rho(igrid_level)
196 map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
197 END DO
198 END DO
199 ! skip empty sets (not uncommon for soft nuclear basis)
200 IF (npgfa(iset) > 0 .AND. npgfa(jset) > 0) THEN
201 offset = offset + 1
202 END IF
203 !
204 IF (any(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
205 ncoa = npgfa(iset)*ncoset(la_max(iset))
206 sgfa = first_sgfa(1, iset)
207 ncob = npgfa(jset)*ncoset(la_max(jset))
208 sgfb = first_sgfa(1, jset)
209 ! decontract density block
210 CALL dgemm("N", "N", ncoa, nsgfa(jset), nsgfa(iset), &
211 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
212 p_block(sgfa, sgfb), SIZE(p_block, 1), &
213 0.0_dp, work(1, 1), maxco)
214 CALL dgemm("N", "T", ncoa, ncob, nsgfa(jset), &
215 1.0_dp, work(1, 1), maxco, &
216 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
217 0.0_dp, pab(1, 1), maxco)
218 DO ipgf = 1, npgfa(iset)
219 IF (jset == iset) THEN
220 npgf2 = ipgf
221 ELSE
222 npgf2 = npgfa(jset)
223 END IF
224 DO jpgf = 1, npgf2
225 IF (map_it2(ipgf, jpgf)) THEN
226 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
227 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
228 rs_grid => rs_rho(igrid_level)
229
230 na1 = (ipgf - 1)*ncoset(la_max(iset))
231 nb1 = (jpgf - 1)*ncoset(la_max(jset))
232
233 radius = exp_radius_very_extended(la_min=la_min(iset), &
234 la_max=la_max(iset), &
235 lb_min=la_min(jset), &
236 lb_max=la_max(jset), &
237 ra=ra, rb=ra, rp=ra, &
238 zetp=zetp, eps=eps_rho_rspace, &
239 prefactor=zeff, cutoff=1.0_dp)
240
241 IF (jset == iset .AND. jpgf == ipgf) THEN
242 scale = -zeff ! nuclear charge density is positive
243 ELSE
244 scale = -2.0_dp*zeff ! symmetric density matrix
245 END IF
246 CALL collocate_pgf_product( &
247 la_max(iset), zeta(ipgf, iset), la_min(iset), &
248 la_max(jset), zeta(jpgf, jset), la_min(jset), &
249 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
250 scale, pab, na1, nb1, rs_grid, &
251 radius=radius, ga_gb_function=grid_func_ab)
252 END IF
253 END DO
254 END DO
255 END IF
256 END DO
257 END DO
258 END IF
259 END DO
260 DEALLOCATE (map_it2)
261 END DO
262 DEALLOCATE (pab, work)
263
264 NULLIFY (rhoz_cneo_s_gs, rhoz_cneo_s_rs)
265 CALL get_rho0_mpole(rho0_mpole=rho0, &
266 rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
267 rhoz_cneo_s_rs=rhoz_cneo_s_rs)
268
269 CALL pw_zero(rhoz_cneo_s_gs)
270 CALL pw_zero(rhoz_cneo_s_rs)
271
272 DO igrid_level = 1, gridlevel_info%ngrid_levels
273 CALL pw_zero(mgrid_rspace(igrid_level))
274 CALL transfer_rs2pw(rs=rs_rho(igrid_level), &
275 pw=mgrid_rspace(igrid_level))
276 END DO
277
278 DO igrid_level = 1, gridlevel_info%ngrid_levels
279 CALL pw_zero(mgrid_gspace(igrid_level))
280 CALL pw_transfer(mgrid_rspace(igrid_level), &
281 mgrid_gspace(igrid_level))
282 CALL pw_axpy(mgrid_gspace(igrid_level), rhoz_cneo_s_gs)
283 END DO
284 CALL pw_transfer(rhoz_cneo_s_gs, rhoz_cneo_s_rs)
285 tot_rs_int = pw_integrate_function(rhoz_cneo_s_rs, isign=-1)
286
287 ! *** give back the multi-grids *** !
288 CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
289 CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
290 ELSE
291 tot_rs_int = 0.0_dp
292 END IF
293
294 CALL timestop(handle)
295
296 END SUBROUTINE put_rhoz_cneo_s_on_grid
297
298! **************************************************************************************************
299!> \brief ...
300!> \param pw_env ...
301!> \param rho0_mpole ...
302! **************************************************************************************************
303 SUBROUTINE rhoz_cneo_s_grid_create(pw_env, rho0_mpole)
304
305 TYPE(pw_env_type), POINTER :: pw_env
306 TYPE(rho0_mpole_type), POINTER :: rho0_mpole
307
308 CHARACTER(len=*), PARAMETER :: routinen = 'rhoz_cneo_s_grid_create'
309
310 INTEGER :: handle
311 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
312
313 CALL timeset(routinen, handle)
314
315 cpassert(ASSOCIATED(pw_env))
316
317 NULLIFY (auxbas_pw_pool)
318 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
319 cpassert(ASSOCIATED(auxbas_pw_pool))
320
321 ! reallocate rho0 on the global grid in real and reciprocal space
322 cpassert(ASSOCIATED(rho0_mpole))
323
324 ! soft nuclear charge density in real space
325 IF (ASSOCIATED(rho0_mpole%rhoz_cneo_s_rs)) THEN
326 CALL rho0_mpole%rhoz_cneo_s_rs%release()
327 ELSE
328 ALLOCATE (rho0_mpole%rhoz_cneo_s_rs)
329 END IF
330 CALL auxbas_pw_pool%create_pw(rho0_mpole%rhoz_cneo_s_rs)
331
332 ! soft nuclear charge density in reciprocal space
333 IF (ASSOCIATED(rho0_mpole%rhoz_cneo_s_gs)) THEN
334 CALL rho0_mpole%rhoz_cneo_s_gs%release()
335 ELSE
336 ALLOCATE (rho0_mpole%rhoz_cneo_s_gs)
337 END IF
338 CALL auxbas_pw_pool%create_pw(rho0_mpole%rhoz_cneo_s_gs)
339
340 CALL timestop(handle)
341
342 END SUBROUTINE rhoz_cneo_s_grid_create
343
344! **************************************************************************************************
345!> \brief ...
346!> \param qs_env ...
347!> \param v_rspace ...
348!> \param para_env ...
349!> \param calculate_forces ...
350!> \param rhoz_cneo_set ...
351!> \param kforce ...
352! **************************************************************************************************
353 SUBROUTINE integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, rhoz_cneo_set, &
354 kforce)
355
356 TYPE(qs_environment_type), POINTER :: qs_env
357 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
358 TYPE(mp_para_env_type), POINTER :: para_env
359 LOGICAL, INTENT(IN) :: calculate_forces
360 TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
361 REAL(kind=dp), INTENT(IN), OPTIONAL :: kforce
362
363 CHARACTER(LEN=*), PARAMETER :: routinen = 'integrate_vhgg_rspace'
364
365 INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
366 m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, npgf2, nseta, offset, sgfa, &
367 sgfb
368 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
369 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
370 LOGICAL :: paw_atom, use_virial
371 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
372 REAL(kind=dp) :: eps_rho_rspace, f0, fscale, radius, &
373 zeff, zetp
374 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
375 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
376 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, p_block, pab, sphi_a, vmat, work, &
377 zeta
378 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
379 TYPE(cell_type), POINTER :: cell
380 TYPE(cneo_potential_type), POINTER :: cneo_potential
381 TYPE(dft_control_type), POINTER :: dft_control
382 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
383 TYPE(gto_basis_set_type), POINTER :: nuc_soft_basis
384 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
385 TYPE(pw_env_type), POINTER :: pw_env
386 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
387 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
388 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
389 TYPE(realspace_grid_type), POINTER :: rs_grid
390 TYPE(virial_type), POINTER :: virial
391
392 CALL timeset(routinen, handle)
393
394 NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
395 first_sgfa, force, gridlevel_info, hab, la_max, la_min, &
396 nuc_soft_basis, npgfa, nsgfa, p_block, pab, particle_set, &
397 pw_env, rs_grid, rs_rho, sphi_a, virial, vmat, work, zeta)
398
399 CALL get_qs_env(qs_env=qs_env, &
400 atomic_kind_set=atomic_kind_set, &
401 qs_kind_set=qs_kind_set, &
402 cell=cell, &
403 dft_control=dft_control, &
404 force=force, pw_env=pw_env, &
405 particle_set=particle_set, &
406 virial=virial)
407
408 ! find maximum numbers
409 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
410 maxco=maxco, &
411 maxsgf_set=maxsgf_set, &
412 basis_type="NUC_SOFT")
413 IF (maxco > 0) THEN
414 fscale = 1.0_dp
415 IF (PRESENT(kforce)) THEN
416 fscale = kforce
417 END IF
418 ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set), hab(maxco, maxco))
419 pab = 0.0_dp
420
421 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
422 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
423
424 cpassert(ASSOCIATED(pw_env))
425 CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, gridlevel_info=gridlevel_info)
426
427 ! transform the potential on the rs_multigrids
428 CALL potential_pw2rs(rs_rho, v_rspace, pw_env)
429
430 offset = 0
431 my_pos = rs_rho(1)%desc%my_pos
432 group_size = rs_rho(1)%desc%group_size
433
434 DO ikind = 1, SIZE(atomic_kind_set, 1)
435 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
436 IF (.NOT. paw_atom) cycle
437
438 NULLIFY (cneo_potential)
439 CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
440 IF (.NOT. ASSOCIATED(cneo_potential)) cycle
441
442 NULLIFY (atom_list)
443 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
444 NULLIFY (nuc_soft_basis)
445 CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_soft_basis, basis_type="NUC_SOFT")
446 CALL get_gto_basis_set(gto_basis_set=nuc_soft_basis, lmax=la_max, &
447 lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
448 sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
449 CALL get_cneo_potential(cneo_potential, zeff=zeff)
450 m1 = maxval(npgfa(1:nseta))
451 ALLOCATE (map_it2(m1, m1))
452
453 DO iatom = 1, natom
454 atom_a = atom_list(iatom)
455 ra(:) = pbc(particle_set(atom_a)%r, cell)
456 IF (rhoz_cneo_set(atom_a)%ready) THEN
457 p_block => rhoz_cneo_set(atom_a)%pmat
458 ELSE
459 NULLIFY (p_block)
460 END IF
461 vmat => rhoz_cneo_set(atom_a)%vmat
462 vmat = 0.0_dp
463 DO iset = 1, nseta
464 DO jset = 1, iset
465 ! processor mapping
466 map_it2 = .false.
467 DO ipgf = 1, npgfa(iset)
468 IF (jset == iset) THEN
469 npgf2 = ipgf
470 ELSE
471 npgf2 = npgfa(jset)
472 END IF
473 DO jpgf = 1, npgf2
474 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
475 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
476 rs_grid => rs_rho(igrid_level)
477 map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
478 END DO
479 END DO
480 ! skip empty sets (not uncommon for soft nuclear basis)
481 IF (npgfa(iset) > 0 .AND. npgfa(jset) > 0) THEN
482 offset = offset + 1
483 END IF
484 !
485 IF (any(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
486 hab = 0.0_dp
487 IF (calculate_forces) THEN
488 force_a = 0.0_dp
489 force_b = 0.0_dp
490 END IF
491 IF (use_virial) THEN
492 my_virial_a = 0.0_dp
493 my_virial_b = 0.0_dp
494 END IF
495 ncoa = npgfa(iset)*ncoset(la_max(iset))
496 sgfa = first_sgfa(1, iset)
497 ncob = npgfa(jset)*ncoset(la_max(jset))
498 sgfb = first_sgfa(1, jset)
499 IF (calculate_forces .AND. ASSOCIATED(p_block)) THEN
500 ! decontract density block
501 CALL dgemm("N", "N", ncoa, nsgfa(jset), nsgfa(iset), &
502 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
503 p_block(sgfa, sgfb), SIZE(p_block, 1), &
504 0.0_dp, work(1, 1), maxco)
505 CALL dgemm("N", "T", ncoa, ncob, nsgfa(jset), &
506 1.0_dp, work(1, 1), maxco, &
507 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
508 0.0_dp, pab(1, 1), maxco)
509 END IF
510 DO ipgf = 1, npgfa(iset)
511 IF (jset == iset) THEN
512 npgf2 = ipgf
513 ELSE
514 npgf2 = npgfa(jset)
515 END IF
516 DO jpgf = 1, npgf2
517 IF (map_it2(ipgf, jpgf)) THEN
518 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
519 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
520 rs_grid => rs_rho(igrid_level)
521
522 na1 = (ipgf - 1)*ncoset(la_max(iset))
523 nb1 = (jpgf - 1)*ncoset(la_max(jset))
524
525 radius = exp_radius_very_extended(la_min=la_min(iset), &
526 la_max=la_max(iset), &
527 lb_min=la_min(jset), &
528 lb_max=la_max(jset), &
529 ra=ra, rb=ra, rp=ra, &
530 zetp=zetp, eps=eps_rho_rspace, &
531 prefactor=zeff, cutoff=1.0_dp)
532
533 CALL integrate_pgf_product( &
534 la_max(iset), zeta(ipgf, iset), la_min(iset), &
535 la_max(jset), zeta(jpgf, jset), la_min(jset), &
536 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), rs_grid, &
537 hab, pab=pab, o1=na1, o2=nb1, &
538 radius=radius, &
539 calculate_forces=calculate_forces, force_a=force_a, force_b=force_b, &
540 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
541 f0 = 0.0_dp
542 IF (calculate_forces .OR. use_virial) THEN
543 IF (jset == iset .AND. jpgf == ipgf) THEN
544 f0 = -fscale*zeff
545 ELSE
546 f0 = -2.0_dp*fscale*zeff
547 END IF
548 END IF
549 IF (calculate_forces) THEN
550 force(ikind)%rho_cneo_nuc(1:3, iatom) = &
551 force(ikind)%rho_cneo_nuc(1:3, iatom) + f0*(force_a + force_b)
552 END IF
553 IF (use_virial) THEN
554 virial%pv_virial = virial%pv_virial + f0*(my_virial_a + my_virial_b)
555 END IF
556 ! symmetrize
557 IF (iset == jset .AND. ipgf /= jpgf) THEN
558 hab(nb1 + 1:nb1 + ncoset(la_max(jset)), &
559 na1 + 1:na1 + ncoset(la_max(iset))) = &
560 transpose(hab(na1 + 1:na1 + ncoset(la_max(iset)), &
561 nb1 + 1:nb1 + ncoset(la_max(jset))))
562 END IF
563 END IF
564 END DO
565 END DO
566 ! contract the soft basis V_Hartree integral
567 work(1:ncoa, 1:nsgfa(jset)) = matmul(hab(1:ncoa, 1:ncob), &
568 sphi_a(1:ncob, sgfb:sgfb + nsgfa(jset) - 1))
569 vmat(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfa(jset) - 1) = &
570 vmat(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfa(jset) - 1) - zeff* &
571 matmul(transpose(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), &
572 work(1:ncoa, 1:nsgfa(jset)))
573 ! symmetrize
574 IF (iset /= jset) THEN
575 vmat(sgfb:sgfb + nsgfa(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
576 transpose(vmat(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfa(jset) - 1))
577 END IF
578 END IF
579 END DO
580 END DO
581 END DO
582 DEALLOCATE (map_it2)
583 END DO
584 DEALLOCATE (pab, work, hab)
585
586 DO ikind = 1, SIZE(atomic_kind_set, 1)
587 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
588 IF (.NOT. paw_atom) cycle
589
590 NULLIFY (cneo_potential)
591 CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
592 IF (.NOT. ASSOCIATED(cneo_potential)) cycle
593
594 NULLIFY (atom_list)
595 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
596
597 DO iatom = 1, natom
598 atom_a = atom_list(iatom)
599 vmat => rhoz_cneo_set(atom_a)%vmat
600 CALL para_env%sum(vmat)
601 END DO
602 END DO
603 END IF
604
605 CALL timestop(handle)
606
607 END SUBROUTINE integrate_vhgg_rspace
608
609END MODULE qs_cneo_ggrid
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:214
cell_types::pbc
Definition cell_types.F:92
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_integrate_function
Definition pw_methods.F:115
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
pw_pool_types::pw_pools_create_pws
Definition pw_pool_types.F:139
pw_pool_types::pw_pools_give_back_pws
Definition pw_pool_types.F:150
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition ao_util.F:209
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
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
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:630
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
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::grid_func_ab
integer, parameter, public grid_func_ab
Definition grid_api.F:29
grid_api::integrate_pgf_product
subroutine, public integrate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, rsgrid, hab, pab, o1, o2, radius, calculate_forces, force_a, force_b, compute_tau, use_virial, my_virial_a, my_virial_b, hdab, hadb, a_hdab, use_subpatch, subpatch_pattern)
low level function to compute matrix elements of primitive gaussian functions
Definition grid_api.F:279
grid_api::collocate_pgf_product
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
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_cneo_ggrid
CNEO soft nuclear densities on the global grid (see J. Chem. Theory Comput. 2025, 21,...
Definition qs_cneo_ggrid.F:15
qs_cneo_ggrid::put_rhoz_cneo_s_on_grid
subroutine, public put_rhoz_cneo_s_on_grid(qs_env, rho0, rhoz_cneo_set, tot_rs_int)
...
Definition qs_cneo_ggrid.F:82
qs_cneo_ggrid::integrate_vhgg_rspace
subroutine, public integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, rhoz_cneo_set, kforce)
...
Definition qs_cneo_ggrid.F:355
qs_cneo_ggrid::rhoz_cneo_s_grid_create
subroutine, public rhoz_cneo_s_grid_create(pw_env, rho0_mpole)
...
Definition qs_cneo_ggrid.F:304
qs_cneo_types
Types used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
Definition qs_cneo_types.F:15
qs_cneo_types::get_cneo_potential
subroutine, public get_cneo_potential(potential, z, zeff, mass, elec_conf, nsgf, nne, npsgf, nsotot, my_gcc_h, my_gcc_s, ovlp, kin, utrans, distance, harmonics, qlm_gg, gg, vgg, n2oindex, o2nindex, rad2l, oorad2l)
...
Definition qs_cneo_types.F:264
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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.
Definition qs_environment_types.F:526
qs_force_types
Definition qs_force_types.F:13
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
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, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:466
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:908
qs_rho0_types
Definition qs_rho0_types.F:8
qs_rho0_types::get_rho0_mpole
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)
...
Definition qs_rho0_types.F:449
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::map_gaussian_here
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
Definition realspace_grid_types.F:2144
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition realspace_grid_types.F:652
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1905
rs_pw_interface
Transfers densities from PW to RS grids and potentials from PW to RS.
Definition rs_pw_interface.F:14
rs_pw_interface::potential_pw2rs
subroutine, public potential_pw2rs(rs_v, v_rspace, pw_env)
transfers a potential from a pw_grid to a vector of realspace multigrids
Definition rs_pw_interface.F:152
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:598
gaussian_gridlevels::gridlevel_info_type
Definition gaussian_gridlevels.F:34
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_cneo_types::cneo_potential_type
Definition qs_cneo_types.F:53
qs_cneo_types::rhoz_cneo_type
Definition qs_cneo_types.F:29
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
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:177
qs_rho0_types::rho0_mpole_type
Definition qs_rho0_types.F:52
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136
virial_types::virial_type
Definition virial_types.F:26