Line data Source code
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 : ! **************************************************************************************************
15 : MODULE 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 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param qs_env ...
77 : !> \param rho0 ...
78 : !> \param rhoz_cneo_set ...
79 : !> \param tot_rs_int ...
80 : ! **************************************************************************************************
81 48 : 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 48 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
94 48 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
95 : LOGICAL :: paw_atom
96 48 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
97 : REAL(KIND=dp) :: eps_rho_rspace, radius, scale, zeff, zetp
98 : REAL(KIND=dp), DIMENSION(3) :: ra
99 48 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, work, zeta
100 48 : 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 48 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
107 48 : 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 48 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
111 48 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
112 : TYPE(pw_r3d_rs_type), POINTER :: rhoz_cneo_s_rs
113 48 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
114 48 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
115 : TYPE(realspace_grid_type), POINTER :: rs_grid
116 :
117 48 : CALL timeset(routineN, handle)
118 :
119 48 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
120 48 : first_sgfa, gridlevel_info, la_max, la_min, nuc_soft_basis, &
121 48 : npgfa, nsgfa, p_block, pab, particle_set, pw_env, pw_pools, &
122 48 : 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 48 : 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 48 : basis_type="NUC_SOFT")
135 48 : IF (maxco > 0) THEN
136 336 : ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set))
137 :
138 48 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
139 :
140 : ! *** set up the pw multi-grids *** !
141 48 : CPASSERT(ASSOCIATED(pw_env))
142 : CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
143 48 : gridlevel_info=gridlevel_info)
144 :
145 48 : CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
146 :
147 48 : CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
148 :
149 : ! *** set up the rs multi-grids *** !
150 114 : DO igrid_level = 1, gridlevel_info%ngrid_levels
151 114 : CALL rs_grid_zero(rs_rho(igrid_level))
152 : END DO
153 :
154 48 : offset = 0
155 48 : my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
156 48 : group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
157 :
158 140 : DO ikind = 1, SIZE(atomic_kind_set)
159 92 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
160 92 : IF (.NOT. paw_atom) CYCLE
161 :
162 86 : NULLIFY (cneo_potential)
163 86 : CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
164 86 : IF (.NOT. ASSOCIATED(cneo_potential)) CYCLE
165 :
166 48 : NULLIFY (atom_list)
167 48 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
168 48 : NULLIFY (nuc_soft_basis)
169 48 : 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 48 : sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
173 48 : CALL get_cneo_potential(cneo_potential, zeff=zeff)
174 480 : m1 = MAXVAL(npgfa(1:nseta))
175 192 : ALLOCATE (map_it2(m1, m1))
176 :
177 140 : DO iatom = 1, natom
178 92 : atom_a = atom_list(iatom)
179 140 : IF (rhoz_cneo_set(atom_a)%ready) THEN
180 78 : ra(:) = pbc(particle_set(atom_a)%r, cell)
181 78 : p_block => rhoz_cneo_set(atom_a)%pmat
182 780 : DO iset = 1, nseta
183 4290 : DO jset = 1, iset
184 : ! processor mapping
185 10530 : map_it2 = .FALSE.
186 3708 : DO ipgf = 1, npgfa(iset)
187 198 : IF (jset == iset) THEN
188 : npgf2 = ipgf
189 : ELSE
190 96 : npgf2 = npgfa(jset)
191 : END IF
192 3858 : DO jpgf = 1, npgf2
193 150 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
194 150 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
195 150 : rs_grid => rs_rho(igrid_level)
196 348 : 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 3510 : IF (npgfa(iset) > 0 .AND. npgfa(jset) > 0) THEN
201 150 : offset = offset + 1
202 : END IF
203 : !
204 5034 : IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
205 75 : ncoa = npgfa(iset)*ncoset(la_max(iset))
206 75 : sgfa = first_sgfa(1, iset)
207 75 : ncob = npgfa(jset)*ncoset(la_max(jset))
208 75 : 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 75 : 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 75 : 0.0_dp, pab(1, 1), maxco)
218 150 : DO ipgf = 1, npgfa(iset)
219 75 : IF (jset == iset) THEN
220 : npgf2 = ipgf
221 : ELSE
222 24 : npgf2 = npgfa(jset)
223 : END IF
224 225 : DO jpgf = 1, npgf2
225 150 : IF (map_it2(ipgf, jpgf)) THEN
226 75 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
227 75 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
228 75 : rs_grid => rs_rho(igrid_level)
229 :
230 75 : na1 = (ipgf - 1)*ncoset(la_max(iset))
231 75 : 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 75 : prefactor=zeff, cutoff=1.0_dp)
240 :
241 75 : IF (jset == iset .AND. jpgf == ipgf) THEN
242 51 : scale = -zeff ! nuclear charge density is positive
243 : ELSE
244 24 : 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 75 : 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 328 : DEALLOCATE (map_it2)
261 : END DO
262 48 : DEALLOCATE (pab, work)
263 :
264 48 : 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 48 : rhoz_cneo_s_rs=rhoz_cneo_s_rs)
268 :
269 48 : CALL pw_zero(rhoz_cneo_s_gs)
270 48 : CALL pw_zero(rhoz_cneo_s_rs)
271 :
272 114 : DO igrid_level = 1, gridlevel_info%ngrid_levels
273 66 : CALL pw_zero(mgrid_rspace(igrid_level))
274 : CALL transfer_rs2pw(rs=rs_rho(igrid_level), &
275 114 : pw=mgrid_rspace(igrid_level))
276 : END DO
277 :
278 114 : DO igrid_level = 1, gridlevel_info%ngrid_levels
279 66 : CALL pw_zero(mgrid_gspace(igrid_level))
280 : CALL pw_transfer(mgrid_rspace(igrid_level), &
281 66 : mgrid_gspace(igrid_level))
282 114 : CALL pw_axpy(mgrid_gspace(igrid_level), rhoz_cneo_s_gs)
283 : END DO
284 48 : CALL pw_transfer(rhoz_cneo_s_gs, rhoz_cneo_s_rs)
285 48 : tot_rs_int = pw_integrate_function(rhoz_cneo_s_rs, isign=-1)
286 :
287 : ! *** give back the multi-grids *** !
288 48 : CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
289 144 : CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
290 : ELSE
291 0 : tot_rs_int = 0.0_dp
292 : END IF
293 :
294 48 : CALL timestop(handle)
295 :
296 96 : END SUBROUTINE put_rhoz_cneo_s_on_grid
297 :
298 : ! **************************************************************************************************
299 : !> \brief ...
300 : !> \param pw_env ...
301 : !> \param rho0_mpole ...
302 : ! **************************************************************************************************
303 16 : 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 16 : CALL timeset(routineN, handle)
314 :
315 16 : CPASSERT(ASSOCIATED(pw_env))
316 :
317 16 : NULLIFY (auxbas_pw_pool)
318 16 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
319 16 : CPASSERT(ASSOCIATED(auxbas_pw_pool))
320 :
321 : ! reallocate rho0 on the global grid in real and reciprocal space
322 16 : CPASSERT(ASSOCIATED(rho0_mpole))
323 :
324 : ! soft nuclear charge density in real space
325 16 : IF (ASSOCIATED(rho0_mpole%rhoz_cneo_s_rs)) THEN
326 8 : CALL rho0_mpole%rhoz_cneo_s_rs%release()
327 : ELSE
328 8 : ALLOCATE (rho0_mpole%rhoz_cneo_s_rs)
329 : END IF
330 16 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rhoz_cneo_s_rs)
331 :
332 : ! soft nuclear charge density in reciprocal space
333 16 : IF (ASSOCIATED(rho0_mpole%rhoz_cneo_s_gs)) THEN
334 8 : CALL rho0_mpole%rhoz_cneo_s_gs%release()
335 : ELSE
336 8 : ALLOCATE (rho0_mpole%rhoz_cneo_s_gs)
337 : END IF
338 16 : CALL auxbas_pw_pool%create_pw(rho0_mpole%rhoz_cneo_s_gs)
339 :
340 16 : CALL timestop(handle)
341 :
342 16 : 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 48 : 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 48 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
369 48 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
370 : LOGICAL :: paw_atom, use_virial
371 48 : 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 48 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, p_block, pab, sphi_a, vmat, work, &
377 48 : zeta
378 48 : 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 48 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
385 : TYPE(pw_env_type), POINTER :: pw_env
386 48 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
387 48 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
388 48 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
389 : TYPE(realspace_grid_type), POINTER :: rs_grid
390 : TYPE(virial_type), POINTER :: virial
391 :
392 48 : CALL timeset(routineN, handle)
393 :
394 48 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
395 48 : first_sgfa, force, gridlevel_info, hab, la_max, la_min, &
396 48 : nuc_soft_basis, npgfa, nsgfa, p_block, pab, particle_set, &
397 48 : 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 48 : 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 48 : basis_type="NUC_SOFT")
413 48 : IF (maxco > 0) THEN
414 48 : fscale = 1.0_dp
415 48 : IF (PRESENT(kforce)) THEN
416 0 : fscale = kforce
417 : END IF
418 432 : ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set), hab(maxco, maxco))
419 792 : pab = 0.0_dp
420 :
421 48 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
422 48 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
423 :
424 48 : CPASSERT(ASSOCIATED(pw_env))
425 48 : 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 48 : CALL potential_pw2rs(rs_rho, v_rspace, pw_env)
429 :
430 48 : offset = 0
431 48 : my_pos = rs_rho(1)%desc%my_pos
432 48 : group_size = rs_rho(1)%desc%group_size
433 :
434 140 : DO ikind = 1, SIZE(atomic_kind_set, 1)
435 92 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
436 92 : IF (.NOT. paw_atom) CYCLE
437 :
438 86 : NULLIFY (cneo_potential)
439 86 : CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
440 86 : IF (.NOT. ASSOCIATED(cneo_potential)) CYCLE
441 :
442 48 : NULLIFY (atom_list)
443 48 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
444 48 : NULLIFY (nuc_soft_basis)
445 48 : 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 48 : sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
449 48 : CALL get_cneo_potential(cneo_potential, zeff=zeff)
450 480 : m1 = MAXVAL(npgfa(1:nseta))
451 192 : ALLOCATE (map_it2(m1, m1))
452 :
453 140 : DO iatom = 1, natom
454 92 : atom_a = atom_list(iatom)
455 92 : ra(:) = pbc(particle_set(atom_a)%r, cell)
456 92 : IF (rhoz_cneo_set(atom_a)%ready) THEN
457 78 : p_block => rhoz_cneo_set(atom_a)%pmat
458 : ELSE
459 : NULLIFY (p_block)
460 : END IF
461 92 : vmat => rhoz_cneo_set(atom_a)%vmat
462 50876 : vmat = 0.0_dp
463 968 : DO iset = 1, nseta
464 5060 : DO jset = 1, iset
465 : ! processor mapping
466 12420 : map_it2 = .FALSE.
467 4412 : DO ipgf = 1, npgfa(iset)
468 272 : IF (jset == iset) THEN
469 : npgf2 = ipgf
470 : ELSE
471 144 : npgf2 = npgfa(jset)
472 : END IF
473 4612 : DO jpgf = 1, npgf2
474 200 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
475 200 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
476 200 : rs_grid => rs_rho(igrid_level)
477 472 : 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 4140 : IF (npgfa(iset) > 0 .AND. npgfa(jset) > 0) THEN
482 200 : offset = offset + 1
483 : END IF
484 : !
485 5976 : IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
486 6780 : hab = 0.0_dp
487 100 : IF (calculate_forces) THEN
488 24 : force_a = 0.0_dp
489 24 : force_b = 0.0_dp
490 : END IF
491 100 : IF (use_virial) THEN
492 0 : my_virial_a = 0.0_dp
493 0 : my_virial_b = 0.0_dp
494 : END IF
495 100 : ncoa = npgfa(iset)*ncoset(la_max(iset))
496 100 : sgfa = first_sgfa(1, iset)
497 100 : ncob = npgfa(jset)*ncoset(la_max(jset))
498 100 : sgfb = first_sgfa(1, jset)
499 100 : 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 24 : 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 24 : 0.0_dp, pab(1, 1), maxco)
509 : END IF
510 200 : DO ipgf = 1, npgfa(iset)
511 100 : IF (jset == iset) THEN
512 : npgf2 = ipgf
513 : ELSE
514 36 : npgf2 = npgfa(jset)
515 : END IF
516 300 : DO jpgf = 1, npgf2
517 200 : IF (map_it2(ipgf, jpgf)) THEN
518 100 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
519 100 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
520 100 : rs_grid => rs_rho(igrid_level)
521 :
522 100 : na1 = (ipgf - 1)*ncoset(la_max(iset))
523 100 : 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 100 : 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 100 : use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
541 100 : f0 = 0.0_dp
542 100 : IF (calculate_forces .OR. use_virial) THEN
543 24 : IF (jset == iset .AND. jpgf == ipgf) THEN
544 12 : f0 = -fscale*zeff
545 : ELSE
546 12 : f0 = -2.0_dp*fscale*zeff
547 : END IF
548 : END IF
549 100 : IF (calculate_forces) THEN
550 : force(ikind)%rho_cneo_nuc(1:3, iatom) = &
551 96 : force(ikind)%rho_cneo_nuc(1:3, iatom) + f0*(force_a + force_b)
552 : END IF
553 100 : IF (use_virial) THEN
554 0 : virial%pv_virial = virial%pv_virial + f0*(my_virial_a + my_virial_b)
555 : END IF
556 : ! symmetrize
557 100 : 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 0 : 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 100 : work(1:ncoa, 1:nsgfa(jset)) = MATMUL(hab(1:ncoa, 1:ncob), &
568 7960 : 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 400 : MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), &
572 5390 : work(1:ncoa, 1:nsgfa(jset)))
573 : ! symmetrize
574 100 : IF (iset /= jset) THEN
575 : vmat(sgfb:sgfb + nsgfa(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
576 684 : 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 328 : DEALLOCATE (map_it2)
583 : END DO
584 48 : DEALLOCATE (pab, work, hab)
585 :
586 140 : DO ikind = 1, SIZE(atomic_kind_set, 1)
587 92 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
588 92 : IF (.NOT. paw_atom) CYCLE
589 :
590 86 : NULLIFY (cneo_potential)
591 86 : CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
592 86 : IF (.NOT. ASSOCIATED(cneo_potential)) CYCLE
593 :
594 48 : NULLIFY (atom_list)
595 48 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
596 :
597 280 : DO iatom = 1, natom
598 92 : atom_a = atom_list(iatom)
599 92 : vmat => rhoz_cneo_set(atom_a)%vmat
600 101752 : CALL para_env%sum(vmat)
601 : END DO
602 : END DO
603 : END IF
604 :
605 48 : CALL timestop(handle)
606 :
607 96 : END SUBROUTINE integrate_vhgg_rspace
608 :
609 200 : END MODULE qs_cneo_ggrid
|
