Line data Source code
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 : MODULE hartree_local_methods
9 : USE atomic_kind_types, ONLY: atomic_kind_type,&
10 : get_atomic_kind
11 : USE basis_set_types, ONLY: get_gto_basis_set,&
12 : gto_basis_set_type
13 : USE cell_types, ONLY: cell_type
14 : USE cp_control_types, ONLY: dft_control_type
15 : USE hartree_local_types, ONLY: allocate_ecoul_1center,&
16 : ecoul_1center_type,&
17 : hartree_local_type,&
18 : set_ecoul_1c
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: fourpi,&
21 : pi
22 : USE message_passing, ONLY: mp_para_env_type
23 : USE orbital_pointers, ONLY: indso,&
24 : nsoset
25 : USE pw_env_types, ONLY: pw_env_get,&
26 : pw_env_type
27 : USE pw_poisson_types, ONLY: pw_poisson_periodic,&
28 : pw_poisson_type
29 : USE qs_charges_types, ONLY: qs_charges_type
30 : USE qs_cneo_methods, ONLY: Vh_1c_nuc_integrals,&
31 : calculate_rhoz_cneo
32 : USE qs_cneo_types, ONLY: cneo_potential_type,&
33 : rhoz_cneo_type
34 : USE qs_environment_types, ONLY: get_qs_env,&
35 : qs_environment_type
36 : USE qs_grid_atom, ONLY: grid_atom_type
37 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
38 : harmonics_atom_type
39 : USE qs_kind_types, ONLY: get_qs_kind,&
40 : get_qs_kind_set,&
41 : qs_kind_type
42 : USE qs_local_rho_types, ONLY: get_local_rho,&
43 : local_rho_type,&
44 : rhoz_type
45 : USE qs_rho0_types, ONLY: get_rho0_mpole,&
46 : rho0_atom_type,&
47 : rho0_mpole_type
48 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
49 : rho_atom_coeff,&
50 : rho_atom_type
51 : USE util, ONLY: get_limit
52 : USE virial_types, ONLY: virial_type
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hartree_local_methods'
60 :
61 : ! Public Subroutine
62 :
63 : PUBLIC :: init_coulomb_local, calculate_Vh_1center, Vh_1c_gg_integrals
64 : PUBLIC :: gapw_background_potential
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief ...
70 : !> \param hartree_local ...
71 : !> \param natom ...
72 : ! **************************************************************************************************
73 3300 : SUBROUTINE init_coulomb_local(hartree_local, natom)
74 :
75 : TYPE(hartree_local_type), POINTER :: hartree_local
76 : INTEGER, INTENT(IN) :: natom
77 :
78 : CHARACTER(len=*), PARAMETER :: routineN = 'init_coulomb_local'
79 :
80 : INTEGER :: handle
81 3300 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
82 :
83 3300 : CALL timeset(routineN, handle)
84 :
85 3300 : NULLIFY (ecoul_1c)
86 : ! Allocate and Initialize 1-center Potentials and Integrals
87 3300 : CALL allocate_ecoul_1center(ecoul_1c, natom)
88 3300 : hartree_local%ecoul_1c => ecoul_1c
89 :
90 3300 : CALL timestop(handle)
91 :
92 3300 : END SUBROUTINE init_coulomb_local
93 :
94 : ! **************************************************************************************************
95 : !> \brief Calculates Hartree potential for hard and soft densities (including
96 : !> nuclear charge and compensation charges) using numerical integration
97 : !> \param vrad_h ...
98 : !> \param vrad_s ...
99 : !> \param rrad_h ...
100 : !> \param rrad_s ...
101 : !> \param rrad_0 ...
102 : !> \param rrad_z ...
103 : !> \param grid_atom ...
104 : !> \par History
105 : !> 05.2012 JGH refactoring
106 : !> \author ??
107 : ! **************************************************************************************************
108 15 : SUBROUTINE calculate_Vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom)
109 :
110 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vrad_h, vrad_s
111 : TYPE(rho_atom_coeff), DIMENSION(:), INTENT(IN) :: rrad_h, rrad_s
112 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: rrad_0
113 : REAL(dp), DIMENSION(:), INTENT(IN) :: rrad_z
114 : TYPE(grid_atom_type), POINTER :: grid_atom
115 :
116 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_Vh_1center'
117 :
118 : INTEGER :: handle, ir, iso, ispin, l_ang, &
119 : max_s_harm, nchannels, nr, nspins
120 : REAL(dp) :: I1_down, I1_up, I2_down, I2_up, prefactor
121 15 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: rho_1, rho_2
122 15 : REAL(dp), DIMENSION(:), POINTER :: wr
123 15 : REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l
124 :
125 15 : CALL timeset(routineN, handle)
126 :
127 15 : nr = grid_atom%nr
128 15 : max_s_harm = SIZE(vrad_h, 2)
129 15 : nspins = SIZE(rrad_h, 1)
130 15 : nchannels = SIZE(rrad_0, 2)
131 :
132 15 : r2l => grid_atom%rad2l
133 15 : oor2l => grid_atom%oorad2l
134 15 : wr => grid_atom%wr
135 :
136 60 : ALLOCATE (rho_1(nr), rho_2(nr))
137 :
138 198 : DO iso = 1, max_s_harm
139 183 : rho_1(:) = 0.0_dp
140 183 : rho_2(:) = 0.0_dp
141 948 : IF (iso == 1) rho_1(:) = rrad_z(:)
142 6933 : IF (iso <= nchannels) rho_2(:) = rrad_0(:, iso)
143 549 : DO ispin = 1, nspins
144 18666 : rho_1(:) = rho_1(:) + rrad_h(ispin)%r_coef(:, iso)
145 18849 : rho_2(:) = rho_2(:) + rrad_s(ispin)%r_coef(:, iso)
146 : END DO
147 :
148 183 : l_ang = indso(1, iso)
149 183 : prefactor = fourpi/(2._dp*l_ang + 1._dp)
150 :
151 9333 : rho_1(:) = rho_1(:)*wr(:)
152 9333 : rho_2(:) = rho_2(:)*wr(:)
153 :
154 183 : I1_up = 0.0_dp
155 183 : I1_down = 0.0_dp
156 183 : I2_up = 0.0_dp
157 183 : I2_down = 0.0_dp
158 :
159 183 : I1_up = r2l(nr, l_ang)*rho_1(nr)
160 183 : I2_up = r2l(nr, l_ang)*rho_2(nr)
161 :
162 9150 : DO ir = nr - 1, 1, -1
163 8967 : I1_down = I1_down + oor2l(ir, l_ang + 1)*rho_1(ir)
164 9150 : I2_down = I2_down + oor2l(ir, l_ang + 1)*rho_2(ir)
165 : END DO
166 :
167 : vrad_h(nr, iso) = vrad_h(nr, iso) + prefactor* &
168 183 : (oor2l(nr, l_ang + 1)*I1_up + r2l(nr, l_ang)*I1_down)
169 : vrad_s(nr, iso) = vrad_s(nr, iso) + prefactor* &
170 183 : (oor2l(nr, l_ang + 1)*I2_up + r2l(nr, l_ang)*I2_down)
171 :
172 9165 : DO ir = nr - 1, 1, -1
173 8967 : I1_up = I1_up + r2l(ir, l_ang)*rho_1(ir)
174 8967 : I1_down = I1_down - oor2l(ir, l_ang + 1)*rho_1(ir)
175 8967 : I2_up = I2_up + r2l(ir, l_ang)*rho_2(ir)
176 8967 : I2_down = I2_down - oor2l(ir, l_ang + 1)*rho_2(ir)
177 :
178 : vrad_h(ir, iso) = vrad_h(ir, iso) + prefactor* &
179 8967 : (oor2l(ir, l_ang + 1)*I1_up + r2l(ir, l_ang)*I1_down)
180 : vrad_s(ir, iso) = vrad_s(ir, iso) + prefactor* &
181 9150 : (oor2l(ir, l_ang + 1)*I2_up + r2l(ir, l_ang)*I2_down)
182 :
183 : END DO
184 :
185 : END DO
186 :
187 15 : DEALLOCATE (rho_1, rho_2)
188 :
189 15 : CALL timestop(handle)
190 :
191 15 : END SUBROUTINE calculate_Vh_1center
192 :
193 : ! **************************************************************************************************
194 : !> \brief Calculates one center GAPW Hartree energies and matrix elements
195 : !> Hartree potentials are input
196 : !> Takes possible background charge into account
197 : !> Special case for densities without core charge
198 : !> \param qs_env ...
199 : !> \param energy_hartree_1c ...
200 : !> \param ecoul_1c ...
201 : !> \param local_rho_set ...
202 : !> \param para_env ...
203 : !> \param tddft ...
204 : !> \param local_rho_set_2nd ...
205 : !> \param core_2nd ...
206 : !> \par History
207 : !> 05.2012 JGH refactoring
208 : !> \author ??
209 : ! **************************************************************************************************
210 28722 : SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, &
211 : core_2nd)
212 :
213 : TYPE(qs_environment_type), POINTER :: qs_env
214 : REAL(kind=dp), INTENT(out) :: energy_hartree_1c
215 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
216 : TYPE(local_rho_type), POINTER :: local_rho_set
217 : TYPE(mp_para_env_type), POINTER :: para_env
218 : LOGICAL, INTENT(IN) :: tddft
219 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set_2nd
220 : LOGICAL, INTENT(IN), OPTIONAL :: core_2nd
221 :
222 : CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_gg_integrals'
223 :
224 : INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, &
225 : llmax_nuc, lmax0, lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_iso_not0_nuc, &
226 : max_s_harm, max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, nat, nchan_0, &
227 : nkind, nr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
228 28722 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
229 28722 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
230 28722 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmax_nuc, lmin, &
231 28722 : lmin_nuc, npgf, npgf_nuc
232 : LOGICAL :: cneo, core_charge, l_2nd_local_rho, &
233 : my_core_2nd, my_periodic, paw_atom
234 : REAL(dp) :: back_ch, background_energy, &
235 : background_nuc, ecoul_1_z_cneo
236 28722 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: background_potential, gexp, sqrtwr
237 28722 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_00_nuc, aVh1b_hh, &
238 28722 : aVh1b_hh_nuc, aVh1b_ss, aVh1b_ss_nuc, &
239 28722 : g0_h_w
240 28722 : REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
241 28722 : REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, gsph_nuc, rrad_0, &
242 28722 : Vh1_h, Vh1_s, vrrad_0, zet, zet_nuc
243 28722 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg, Qlm_gg_2nd
244 28722 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
245 : TYPE(cell_type), POINTER :: cell
246 : TYPE(cneo_potential_type), POINTER :: cneo_potential
247 : TYPE(dft_control_type), POINTER :: dft_control
248 : TYPE(grid_atom_type), POINTER :: grid_atom
249 : TYPE(gto_basis_set_type), POINTER :: basis_1c, nuc_basis
250 : TYPE(harmonics_atom_type), POINTER :: harmonics
251 : TYPE(pw_env_type), POINTER :: pw_env
252 : TYPE(pw_poisson_type), POINTER :: poisson_env
253 : TYPE(qs_charges_type), POINTER :: qs_charges
254 28722 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
255 28722 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho0_atom_set_2nd
256 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole, rho0_mpole_2nd
257 28722 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
258 : TYPE(rho_atom_type), POINTER :: rho_atom
259 28722 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
260 : TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
261 28722 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
262 : TYPE(virial_type), POINTER :: virial
263 :
264 28722 : CALL timeset(routineN, handle)
265 :
266 28722 : NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
267 28722 : NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
268 28722 : NULLIFY (rho0_mpole, rhoz_set)
269 28722 : NULLIFY (atom_list, grid_atom, harmonics)
270 28722 : NULLIFY (basis_1c, lmin, lmax, npgf, zet)
271 : NULLIFY (gsph)
272 28722 : NULLIFY (rhoz_cneo, rhoz_cneo_set)
273 :
274 : CALL get_qs_env(qs_env=qs_env, &
275 : cell=cell, dft_control=dft_control, &
276 : atomic_kind_set=atomic_kind_set, &
277 : qs_kind_set=qs_kind_set, &
278 28722 : pw_env=pw_env, qs_charges=qs_charges, virial=virial)
279 :
280 28722 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
281 28722 : my_periodic = (poisson_env%method == pw_poisson_periodic)
282 :
283 : ! The compensating uniform charge is fixed by the electron number, not
284 : ! by the approximate hard/soft quadrature. A planar countercharge makes
285 : ! the full Poisson problem neutral. At fixed charge the correction is
286 : ! linear in the density, so it contributes no TDDFT Hartree kernel.
287 28722 : back_ch = 0.0_dp
288 28722 : IF (my_periodic .AND. .NOT. dft_control%do_pcc .AND. .NOT. tddft) back_ch = dft_control%charge
289 :
290 : ! rhoz_set is not accessed in TDDFT
291 : CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, &
292 28722 : rhoz_cneo_set) ! for integral space
293 :
294 : ! for forces we need a second local_rho_set
295 28722 : l_2nd_local_rho = .FALSE.
296 28722 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
297 348 : l_2nd_local_rho = .TRUE.
298 348 : NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
299 348 : CALL get_local_rho(local_rho_set_2nd, rho_atom_set_2nd, rho0_atom_set_2nd, rho0_mpole_2nd, rhoz_set=rhoz_set_2nd)
300 : END IF
301 :
302 28722 : nkind = SIZE(atomic_kind_set, 1)
303 28722 : nspins = dft_control%nspins
304 :
305 28722 : core_charge = .NOT. tddft ! for forces mixed version
306 28722 : my_core_2nd = .TRUE.
307 28722 : IF (PRESENT(core_2nd)) my_core_2nd = .NOT. core_2nd ! if my_core_2nd true, include core charge
308 :
309 : ! The aim of the following code was to return immediately if the subroutine
310 : ! was called for triplet excited states in spin-restricted case. This check
311 : ! is also performed before invocation of this subroutine. It should be save
312 : ! to remove the optional argument 'do_triplet' from the subroutine interface.
313 : !IF (tddft) THEN
314 : ! CPASSERT(PRESENT(do_triplet))
315 : ! IF (nspins == 1 .AND. do_triplet) RETURN
316 : !END IF
317 :
318 28722 : CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
319 28722 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax_0)
320 :
321 : ! Put to 0 the local hartree energy contribution from 1 center integrals
322 28722 : energy_hartree_1c = 0.0_dp
323 : ! Restore total quantum nuclear density to zero
324 28722 : qs_charges%total_rho1_hard_nuc = 0.0_dp
325 28722 : rho0_mpole%tot_rhoz_cneo_s = 0.0_dp
326 :
327 : ! Here starts the loop over all the atoms
328 83340 : DO ikind = 1, nkind
329 :
330 54618 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
331 54618 : NULLIFY (cneo_potential)
332 : CALL get_qs_kind(qs_kind_set(ikind), &
333 : grid_atom=grid_atom, &
334 : harmonics=harmonics, ngrid_rad=nr, &
335 : max_iso_not0=max_iso_not0, paw_atom=paw_atom, &
336 54618 : cneo_potential=cneo_potential)
337 : CALL get_qs_kind(qs_kind_set(ikind), &
338 54618 : basis_set=basis_1c, basis_type="GAPW_1C")
339 :
340 54618 : cneo = ASSOCIATED(cneo_potential)
341 54618 : IF (cneo .AND. tddft) THEN
342 0 : CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.")
343 : END IF
344 :
345 54618 : NULLIFY (nuc_basis)
346 54618 : max_iso_not0_nuc = 0
347 54618 : IF (cneo) THEN
348 48 : CPASSERT(paw_atom)
349 : CALL get_qs_kind(qs_kind_set(ikind), &
350 48 : basis_set=nuc_basis, basis_type="NUC")
351 48 : max_iso_not0_nuc = cneo_potential%harmonics%max_iso_not0
352 : END IF
353 :
354 137958 : IF (paw_atom) THEN
355 : !=========== PAW ===============
356 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
357 : maxso=maxso, npgf=npgf, maxl=maxl, &
358 51264 : nset=nset, zet=zet)
359 :
360 51264 : max_s_harm = harmonics%max_s_harm
361 51264 : llmax = harmonics%llmax
362 :
363 51264 : nsotot = maxso*nset
364 205056 : ALLOCATE (gsph(nr, nsotot))
365 153792 : ALLOCATE (gexp(nr))
366 307584 : ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0), background_potential(nr))
367 51264 : background_potential(:) = gapw_background_potential(back_ch, cell%deth, grid_atom%rad2(:))
368 :
369 205056 : ALLOCATE (aVh1b_hh(nsotot, nsotot))
370 153792 : ALLOCATE (aVh1b_ss(nsotot, nsotot))
371 153792 : ALLOCATE (aVh1b_00(nsotot, nsotot))
372 :
373 51264 : NULLIFY (Qlm_gg, g0_h)
374 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
375 : l0_ikind=lmax0, &
376 51264 : Qlm_gg=Qlm_gg, g0_h=g0_h) ! Qlm_gg of density
377 :
378 51264 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
379 708 : NULLIFY (Qlm_gg_2nd, g0_h_2nd)
380 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole_2nd, ikind=ikind, &
381 : l0_ikind=lmax0_2nd, &
382 708 : Qlm_gg=Qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
383 : END IF
384 51264 : nchan_0 = nsoset(lmax0)
385 :
386 51264 : IF (nchan_0 > MAX(max_iso_not0, max_iso_not0_nuc)) THEN
387 0 : CPABORT("channels for rho0 > # max of spherical harmonics")
388 : END IF
389 :
390 51264 : NULLIFY (Vh1_h, Vh1_s)
391 205056 : ALLOCATE (Vh1_h(nr, max_iso_not0))
392 205056 : ALLOCATE (Vh1_s(nr, MAX(max_iso_not0, max_iso_not0_nuc, nchan_0)))
393 :
394 51264 : NULLIFY (lmax_nuc, lmin_nuc, npgf_nuc, zet_nuc, gsph_nuc)
395 51264 : maxso_nuc = 0
396 51264 : maxl_nuc = -1
397 51264 : nset_nuc = 0
398 51264 : max_s_harm_nuc = 0
399 51264 : llmax_nuc = -1
400 51264 : nsotot_nuc = 0
401 51264 : IF (cneo) THEN
402 : CALL get_gto_basis_set(gto_basis_set=nuc_basis, lmax=lmax_nuc, &
403 : lmin=lmin_nuc, maxso=maxso_nuc, npgf=npgf_nuc, &
404 48 : maxl=maxl_nuc, nset=nset_nuc, zet=zet_nuc)
405 :
406 48 : max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
407 48 : llmax_nuc = cneo_potential%harmonics%llmax
408 48 : nsotot_nuc = maxso_nuc*nset_nuc
409 192 : ALLOCATE (gsph_nuc(nr, nsotot_nuc))
410 198336 : gsph_nuc = 0.0_dp
411 192 : ALLOCATE (aVh1b_hh_nuc(nsotot_nuc, nsotot_nuc))
412 144 : ALLOCATE (aVh1b_ss_nuc(nsotot_nuc, nsotot_nuc))
413 192 : ALLOCATE (aVh1b_00_nuc(nsotot_nuc, nsotot_nuc))
414 : END IF
415 :
416 0 : ALLOCATE (cg_list(2, nsoset(MAX(maxl, maxl_nuc))**2, &
417 : MAX(max_s_harm, max_s_harm_nuc)), &
418 307584 : cg_n_list(MAX(max_s_harm, max_s_harm_nuc)))
419 :
420 51264 : NULLIFY (rrad_z, my_CG)
421 51264 : my_CG => harmonics%my_CG
422 :
423 : ! set to zero temporary arrays
424 51264 : sqrtwr = 0.0_dp
425 51264 : g0_h_w = 0.0_dp
426 51264 : gexp = 0.0_dp
427 105675970 : gsph = 0.0_dp
428 :
429 2841988 : sqrtwr(1:nr) = SQRT(grid_atom%wr(1:nr))
430 193596 : DO l_ang = 0, lmax0
431 7851840 : g0_h_w(1:nr, l_ang) = g0_h(1:nr, l_ang)*grid_atom%wr(1:nr)
432 : END DO
433 :
434 : m1 = 0
435 182808 : DO iset1 = 1, nset
436 131544 : n1 = nsoset(lmax(iset1))
437 495642 : DO ipgf1 = 1, npgf(iset1)
438 20096714 : gexp(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
439 1518524 : DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1))
440 1022882 : iso = is1 + (ipgf1 - 1)*n1 + m1
441 1022882 : l_ang = indso(1, is1)
442 111233212 : gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
443 : END DO ! is1
444 : END DO ! ipgf1
445 182808 : m1 = m1 + maxso
446 : END DO ! iset1
447 :
448 51264 : IF (cneo) THEN
449 : ! initialize nuclear pmat, cpc and e_core to zero
450 140 : DO iat = 1, nat
451 92 : iatom = atom_list(iat)
452 50876 : rhoz_cneo_set(iatom)%pmat = 0.0_dp
453 50876 : rhoz_cneo_set(iatom)%cpc_h = 0.0_dp
454 50876 : rhoz_cneo_set(iatom)%cpc_s = 0.0_dp
455 92 : rhoz_cneo_set(iatom)%e_core = 0.0_dp
456 140 : rhoz_cneo_set(iatom)%ready = .FALSE.
457 : END DO
458 :
459 : ! calculate nuclear gsph
460 48 : m1 = 0
461 480 : DO iset1 = 1, nset_nuc
462 432 : n1 = nsoset(lmax_nuc(iset1))
463 864 : DO ipgf1 = 1, npgf_nuc(iset1)
464 22032 : gexp(1:nr) = EXP(-zet_nuc(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
465 1968 : DO is1 = nsoset(lmin_nuc(iset1) - 1) + 1, nsoset(lmax_nuc(iset1))
466 1104 : iso = is1 + (ipgf1 - 1)*n1 + m1
467 1104 : l_ang = indso(1, is1)
468 111936 : gsph_nuc(1:nr, iso) = cneo_potential%rad2l(1:nr, l_ang)*gexp(1:nr)
469 : END DO ! is1
470 : END DO ! ipgf1
471 480 : m1 = m1 + maxso_nuc
472 : END DO ! iset1
473 : END IF
474 :
475 : ! Distribute the atoms of this kind
476 51264 : num_pe = para_env%num_pe
477 51264 : mepos = para_env%mepos
478 51264 : bo = get_limit(nat, num_pe, mepos)
479 :
480 91510 : DO iat = bo(1), bo(2) !1,nat
481 40246 : iatom = atom_list(iat)
482 40246 : rho_atom => rho_atom_set(iatom)
483 :
484 40246 : NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
485 40246 : IF (core_charge .AND. .NOT. cneo) THEN
486 33608 : rrad_z => rhoz_set(ikind)%r_coef ! for density
487 : END IF
488 40246 : IF (my_core_2nd .AND. .NOT. cneo) THEN
489 33813 : IF (l_2nd_local_rho) THEN
490 260 : vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
491 : ELSE
492 33553 : vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
493 : END IF
494 : END IF
495 40246 : rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
496 40246 : vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
497 40246 : IF (l_2nd_local_rho) THEN
498 520 : rho_atom => rho_atom_set_2nd(iatom)
499 520 : vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
500 : END IF
501 : CALL Vh_1c_atom_potential(rho_atom, vrrad_0, &
502 : grid_atom, my_core_2nd .AND. .NOT. cneo, & ! core charge for potential (2nd)
503 : vrrad_z, Vh1_h, Vh1_s, &
504 40246 : nchan_0, nspins, max_iso_not0, background_potential)
505 :
506 40246 : IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom) ! rho_atom for density
507 :
508 40246 : ecoul_1_z_cneo = 0.0_dp
509 40246 : background_nuc = 0.0_dp
510 40246 : IF (cneo) THEN
511 46 : rhoz_cneo => rhoz_cneo_set(iatom)
512 : ! Add the soft tail of nuclear Hartree potential to total Vh1_s first.
513 : ! vrho already contains the -zeff factor.
514 1196 : DO iso = 1, max_iso_not0_nuc
515 116196 : Vh1_s(:, iso) = Vh1_s(:, iso) + rhoz_cneo%vrho_rad_s(:, iso)
516 : END DO
517 :
518 : ! Build nuclear 1c integrals according to Vh1_h from electronic density only
519 : ! and Vh1_s from electron density and last step (or initial guess) nuclear density
520 : ! Vh1_h_nuc = -Z*Vh1_h, Vh1_s_nuc = -Z*Vh1_s
521 : CALL Vh_1c_nuc_integrals(rhoz_cneo, cneo_potential%zeff, &
522 : aVh1b_hh_nuc, aVh1b_ss_nuc, aVh1b_00_nuc, Vh1_h, Vh1_s, &
523 : max_iso_not0, max_iso_not0_nuc, &
524 : max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
525 : nset_nuc, npgf_nuc, lmin_nuc, lmax_nuc, nsotot_nuc, maxso_nuc, &
526 : nchan_0, gsph_nuc, g0_h_w, cneo_potential%harmonics%my_CG, &
527 46 : cneo_potential%Qlm_gg)
528 :
529 : ! Solve the nuclear 1c problem
530 : CALL calculate_rhoz_cneo(rhoz_cneo, cneo_potential, cg_list, cg_n_list, nset_nuc, &
531 46 : npgf_nuc, lmin_nuc, lmax_nuc, maxl_nuc, maxso_nuc)
532 : ! slm_int(iso=1) = sqrt(4*Pi), slm_int(iso>1) = 0, without using Lebedev grid
533 : ! when printing, nuclear density is positive, thus needing the minus sign
534 : qs_charges%total_rho1_hard_nuc = qs_charges%total_rho1_hard_nuc - SQRT(fourpi) &
535 2346 : *SUM(rhoz_cneo%rho_rad_h(:, 1)*grid_atom%wr(:))
536 : rho0_mpole%tot_rhoz_cneo_s = rho0_mpole%tot_rhoz_cneo_s - SQRT(fourpi) &
537 2346 : *SUM(rhoz_cneo%rho_rad_s(:, 1)*grid_atom%wr(:))
538 :
539 : ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz_cneo_h and Vh1_s*rhoz_cneo_s.
540 : ! Self-interaction of the quantum nucleus is already removed in ecoul_1_z_cneo.
541 : ! rho already contains the -zeff factor.
542 460 : DO iso = 1, MIN(max_iso_not0, max_iso_not0_nuc)
543 : ecoul_1_z_cneo = ecoul_1_z_cneo + 0.5_dp* &
544 : (SUM(Vh1_h(:, iso)*rhoz_cneo%rho_rad_h(:, iso)*grid_atom%wr(:)) &
545 41860 : - SUM(Vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:)))
546 : END DO
547 782 : DO iso = max_iso_not0 + 1, max_iso_not0_nuc
548 : ecoul_1_z_cneo = ecoul_1_z_cneo - 0.5_dp* &
549 37582 : SUM(Vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:))
550 : END DO
551 :
552 : background_nuc = SUM(background_potential*(rhoz_cneo%rho_rad_h(:, 1) - &
553 2346 : rhoz_cneo%rho_rad_s(:, 1))*grid_atom%wr(:))
554 : ! The background interaction is linear, whereas the generic
555 : ! Coulomb contraction above includes only one half of it.
556 46 : ecoul_1_z_cneo = ecoul_1_z_cneo + 0.5_dp*background_nuc
557 :
558 : ! Add nuclear Hartree potential to total Vh1_h after solving the nuclear 1c problem
559 : ! to avoid nuclear self-interaction.
560 : ! vrho already contains the -zeff factor.
561 : ! Here the min of two max_iso_not0's is chosen, because even when the nuclear one
562 : ! is larger, it is meaningless to let Vh1_h have higher angular momentum components
563 : ! as Vh1_h now is only used for the electronic part.
564 460 : DO iso = 1, MIN(max_iso_not0, max_iso_not0_nuc)
565 41860 : Vh1_h(:, iso) = Vh1_h(:, iso) + rhoz_cneo%vrho_rad_h(:, iso)
566 : END DO
567 : END IF
568 :
569 : CALL Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
570 : grid_atom, iatom, core_charge .AND. .NOT. cneo, & ! core charge for density
571 : rrad_z, Vh1_h, Vh1_s, nchan_0, nspins, max_iso_not0, &
572 40246 : background_potential, background_energy)
573 :
574 40246 : IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)
575 :
576 40246 : IF (cneo) THEN
577 46 : CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1c(iatom)%ecoul_1_z + ecoul_1_z_cneo)
578 46 : energy_hartree_1c = energy_hartree_1c + ecoul_1_z_cneo
579 : END IF
580 :
581 : ! Besides the OCE coefficient derivatives, E_background has an
582 : ! explicit 1/volume dependence. At fixed local radial functions
583 : ! its contribution to -dE/d(strain) is E_background times I.
584 40246 : IF (virial%pv_calculate .AND. .NOT. tddft) THEN
585 2664 : DO iso = 1, 3
586 1998 : virial%pv_virial(iso, iso) = virial%pv_virial(iso, iso) + background_energy + background_nuc
587 2664 : virial%pv_ehartree(iso, iso) = virial%pv_ehartree(iso, iso) + background_energy + background_nuc
588 : END DO
589 : END IF
590 :
591 : CALL Vh_1c_atom_integrals(rho_atom, & ! results (int_local_h and int_local_s) written on rho_atom_2nd
592 : ! int_local_h and int_local_s are used in update_ks_atom
593 : ! on int_local_h mixed core / non-core
594 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
595 : max_s_harm, llmax, cg_list, cg_n_list, &
596 : nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
597 91510 : g0_h_w, my_CG, Qlm_gg) ! Qlm_gg for density from local_rho_set
598 :
599 : END DO ! iat
600 :
601 51264 : DEALLOCATE (aVh1b_hh)
602 51264 : DEALLOCATE (aVh1b_ss)
603 51264 : DEALLOCATE (aVh1b_00)
604 51264 : IF (ALLOCATED(aVh1b_hh_nuc)) DEALLOCATE (aVh1b_hh_nuc)
605 51264 : IF (ALLOCATED(aVh1b_ss_nuc)) DEALLOCATE (aVh1b_ss_nuc)
606 51264 : IF (ALLOCATED(aVh1b_00_nuc)) DEALLOCATE (aVh1b_00_nuc)
607 51264 : DEALLOCATE (Vh1_h, Vh1_s)
608 51264 : DEALLOCATE (cg_list, cg_n_list)
609 51264 : DEALLOCATE (gsph)
610 51264 : IF (ASSOCIATED(gsph_nuc)) DEALLOCATE (gsph_nuc)
611 51264 : DEALLOCATE (gexp)
612 51264 : DEALLOCATE (sqrtwr, g0_h_w, background_potential)
613 :
614 153792 : IF (cneo) THEN
615 : ! broadcast nuclear pmat, cpc and e_core
616 140 : DO iat = 1, nat
617 92 : iatom = atom_list(iat)
618 101660 : CALL para_env%sum(rhoz_cneo_set(iatom)%pmat)
619 101660 : CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_h)
620 101660 : CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_s)
621 92 : CALL para_env%sum(rhoz_cneo_set(iatom)%e_core)
622 140 : rhoz_cneo_set(iatom)%ready = .TRUE.
623 : END DO
624 : END IF
625 : ELSE
626 : !=========== NO PAW ===============
627 : ! This term is taken care of using the core density as in GPW
628 : CYCLE
629 : END IF ! paw
630 : END DO ! ikind
631 :
632 28722 : CALL para_env%sum(energy_hartree_1c)
633 28722 : CALL para_env%sum(qs_charges%total_rho1_hard_nuc)
634 28722 : CALL para_env%sum(rho0_mpole%tot_rhoz_cneo_s)
635 :
636 28722 : CALL timestop(handle)
637 :
638 57444 : END SUBROUTINE Vh_1c_gg_integrals
639 :
640 : ! **************************************************************************************************
641 :
642 : ! **************************************************************************************************
643 : !> \brief ...
644 : !> \param rho_atom ...
645 : !> \param vrrad_0 ...
646 : !> \param grid_atom ...
647 : !> \param core_charge ...
648 : !> \param vrrad_z ...
649 : !> \param Vh1_h ...
650 : !> \param Vh1_s ...
651 : !> \param nchan_0 ...
652 : !> \param nspins ...
653 : !> \param max_iso_not0 ...
654 : !> \param background_potential Unweighted l=0 potential of the uniform background.
655 : ! **************************************************************************************************
656 40246 : SUBROUTINE Vh_1c_atom_potential(rho_atom, vrrad_0, &
657 : grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, &
658 40246 : nchan_0, nspins, max_iso_not0, background_potential)
659 :
660 : TYPE(rho_atom_type), POINTER :: rho_atom
661 : REAL(dp), DIMENSION(:, :), POINTER :: vrrad_0
662 : TYPE(grid_atom_type), POINTER :: grid_atom
663 : LOGICAL, INTENT(IN) :: core_charge
664 : REAL(dp), DIMENSION(:), POINTER :: vrrad_z
665 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
666 : INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
667 : REAL(dp), INTENT(IN) :: background_potential(:)
668 :
669 : INTEGER :: iso, ispin, nr
670 40246 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s
671 :
672 40246 : nr = grid_atom%nr
673 :
674 40246 : NULLIFY (vr_h, vr_s)
675 40246 : CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
676 :
677 29202180 : Vh1_h = 0.0_dp
678 29239716 : Vh1_s = 0.0_dp
679 :
680 3772167 : IF (core_charge) Vh1_h(:, 1) = vrrad_z(:)
681 :
682 355236 : DO iso = 1, nchan_0
683 33817714 : Vh1_s(:, iso) = vrrad_0(:, iso)
684 : END DO
685 :
686 86096 : DO ispin = 1, nspins
687 713058 : DO iso = 1, max_iso_not0
688 34570096 : Vh1_h(:, iso) = Vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
689 34615946 : Vh1_s(:, iso) = Vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
690 : END DO
691 : END DO
692 :
693 2210950 : Vh1_h(:, 1) = Vh1_h(:, 1) + background_potential
694 2210950 : Vh1_s(:, 1) = Vh1_s(:, 1) + background_potential
695 :
696 40246 : END SUBROUTINE Vh_1c_atom_potential
697 :
698 : ! **************************************************************************************************
699 :
700 : ! **************************************************************************************************
701 : !> \brief ...
702 : !> \param energy_hartree_1c ...
703 : !> \param ecoul_1c ...
704 : !> \param rho_atom ...
705 : !> \param rrad_0 ...
706 : !> \param grid_atom ...
707 : !> \param iatom ...
708 : !> \param core_charge ...
709 : !> \param rrad_z ...
710 : !> \param Vh1_h ...
711 : !> \param Vh1_s ...
712 : !> \param nchan_0 ...
713 : !> \param nspins ...
714 : !> \param max_iso_not0 ...
715 : !> \param background_potential Unweighted l=0 potential of the uniform background.
716 : !> \param background_energy Full linear background interaction for the volume derivative.
717 : ! **************************************************************************************************
718 40246 : SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
719 : grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
720 40246 : nchan_0, nspins, max_iso_not0, background_potential, background_energy)
721 :
722 : REAL(dp), INTENT(INOUT) :: energy_hartree_1c
723 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
724 : TYPE(rho_atom_type), POINTER :: rho_atom
725 : REAL(dp), DIMENSION(:, :), POINTER :: rrad_0
726 : TYPE(grid_atom_type), POINTER :: grid_atom
727 : INTEGER, INTENT(IN) :: iatom
728 : LOGICAL, INTENT(IN) :: core_charge
729 : REAL(dp), DIMENSION(:), POINTER :: rrad_z
730 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
731 : INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
732 : REAL(dp), INTENT(IN) :: background_potential(:)
733 : REAL(dp), INTENT(OUT) :: background_energy
734 :
735 : INTEGER :: iso, ispin, nr
736 : REAL(dp) :: bg_0, bg_h, bg_s, bg_z, ecoul_1_0, &
737 : ecoul_1_h, ecoul_1_s, ecoul_1_z
738 40246 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
739 :
740 40246 : nr = grid_atom%nr
741 :
742 40246 : NULLIFY (r_h, r_s)
743 40246 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
744 :
745 : ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz
746 40246 : ecoul_1_z = 0.0_dp
747 40246 : bg_z = 0.0_dp
748 40246 : IF (core_charge) THEN
749 1872412 : ecoul_1_z = 0.5_dp*SUM(Vh1_h(:, 1)*rrad_z(:)*grid_atom%wr(:))
750 1872412 : bg_z = SUM(background_potential*rrad_z(:)*grid_atom%wr(:))
751 : END IF
752 :
753 : ! Calculate the contributions to Ecoul coming from Vh1_s*rho0
754 40246 : ecoul_1_0 = 0.0_dp
755 355236 : DO iso = 1, nchan_0
756 16928980 : ecoul_1_0 = ecoul_1_0 + 0.5_dp*SUM(Vh1_s(:, iso)*rrad_0(:, iso)*grid_atom%wr(:))
757 : END DO
758 2210950 : bg_0 = SUM(background_potential*rrad_0(:, 1)*grid_atom%wr(:))
759 :
760 : ! Calculate the contributions to Ecoul coming from Vh1_h*rho1_h and Vh1_s*rho1_s
761 40246 : ecoul_1_s = 0.0_dp
762 40246 : ecoul_1_h = 0.0_dp
763 40246 : bg_h = 0.0_dp
764 40246 : bg_s = 0.0_dp
765 86096 : DO ispin = 1, nspins
766 2556904 : bg_h = bg_h + SUM(background_potential*r_h(ispin)%r_coef(:, 1)*grid_atom%wr(:))
767 2556904 : bg_s = bg_s + SUM(background_potential*r_s(ispin)%r_coef(:, 1)*grid_atom%wr(:))
768 713058 : DO iso = 1, max_iso_not0
769 34570096 : ecoul_1_s = ecoul_1_s + 0.5_dp*SUM(Vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
770 34615946 : ecoul_1_h = ecoul_1_h + 0.5_dp*SUM(Vh1_h(:, iso)*r_h(ispin)%r_coef(:, iso)*grid_atom%wr(:))
771 : END DO
772 : END DO
773 :
774 : ! Unlike the quadratic Hartree self-energy, the interaction with a fixed
775 : ! uniform background has no factor 1/2. Its potential is already in Vh1.
776 40246 : ecoul_1_z = ecoul_1_z + 0.5_dp*bg_z
777 40246 : ecoul_1_0 = ecoul_1_0 + 0.5_dp*bg_0
778 40246 : ecoul_1_h = ecoul_1_h + 0.5_dp*bg_h
779 40246 : ecoul_1_s = ecoul_1_s + 0.5_dp*bg_s
780 40246 : background_energy = bg_z - bg_0 + bg_h - bg_s
781 :
782 40246 : CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
783 40246 : CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
784 :
785 40246 : energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
786 40246 : energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
787 :
788 40246 : END SUBROUTINE Vh_1c_atom_energy
789 :
790 : ! **************************************************************************************************
791 : !> \brief l=0 potential for the periodic mean-potential correction.
792 : !> For a neutral localized difference delta_rho between hard and soft
793 : !> total charges, integral V(delta_rho) = -2*pi/3 integral r^2 delta_rho.
794 : !> CP2K stores electron-minus-nuclear densities, opposite to net_charge.
795 : !> Radial quadrature weights belong in the contraction, not in this potential.
796 : !> \param net_charge Physical cell charge, positive for missing electrons.
797 : !> \param volume Cell volume in bohr^3.
798 : !> \param rad2 Squared atom-centered radial coordinates.
799 : !> \return Real spherical-harmonic l=0 coefficient of the potential.
800 : ! **************************************************************************************************
801 51278 : PURE FUNCTION gapw_background_potential(net_charge, volume, rad2) RESULT(potential)
802 : REAL(dp), INTENT(IN) :: net_charge, volume, rad2(:)
803 : REAL(dp) :: potential(SIZE(rad2))
804 :
805 2864416 : potential = -(2.0_dp*pi/3.0_dp)*SQRT(fourpi)*net_charge/volume*rad2
806 51278 : END FUNCTION gapw_background_potential
807 :
808 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
809 :
810 : ! **************************************************************************************************
811 : !> \brief ...
812 : !> \param rho_atom ...
813 : !> \param aVh1b_hh ...
814 : !> \param aVh1b_ss ...
815 : !> \param aVh1b_00 ...
816 : !> \param Vh1_h ...
817 : !> \param Vh1_s ...
818 : !> \param max_iso_not0 ...
819 : !> \param max_s_harm ...
820 : !> \param llmax ...
821 : !> \param cg_list ...
822 : !> \param cg_n_list ...
823 : !> \param nset ...
824 : !> \param npgf ...
825 : !> \param lmin ...
826 : !> \param lmax ...
827 : !> \param nsotot ...
828 : !> \param maxso ...
829 : !> \param nspins ...
830 : !> \param nchan_0 ...
831 : !> \param gsph ...
832 : !> \param g0_h_w ...
833 : !> \param my_CG ...
834 : !> \param Qlm_gg ...
835 : ! **************************************************************************************************
836 40246 : SUBROUTINE Vh_1c_atom_integrals(rho_atom, &
837 40246 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
838 40246 : max_s_harm, llmax, cg_list, cg_n_list, &
839 : nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
840 40246 : g0_h_w, my_CG, Qlm_gg)
841 :
842 : TYPE(rho_atom_type), POINTER :: rho_atom
843 : REAL(dp), DIMENSION(:, :) :: aVh1b_hh, aVh1b_ss, aVh1b_00
844 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
845 : INTEGER, INTENT(IN) :: max_iso_not0, max_s_harm, llmax
846 : INTEGER, DIMENSION(:, :, :) :: cg_list
847 : INTEGER, DIMENSION(:) :: cg_n_list
848 : INTEGER, INTENT(IN) :: nset
849 : INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax
850 : INTEGER, INTENT(IN) :: nsotot, maxso, nspins, nchan_0
851 : REAL(dp), DIMENSION(:, :), POINTER :: gsph
852 : REAL(dp), DIMENSION(:, 0:) :: g0_h_w
853 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg
854 :
855 : INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, &
856 : iset2, iso, iso1, iso2, ispin, l_ang, &
857 : m1, m2, max_iso_not0_local, n1, n2, nr
858 : REAL(dp) :: gVg_0, gVg_h, gVg_s
859 40246 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
860 :
861 40246 : NULLIFY (int_local_h, int_local_s)
862 : CALL get_rho_atom(rho_atom=rho_atom, &
863 : ga_Vlocal_gb_h=int_local_h, &
864 40246 : ga_Vlocal_gb_s=int_local_s)
865 :
866 : ! Calculate the integrals of the potential with 2 primitives
867 137657068 : aVh1b_hh = 0.0_dp
868 137657068 : aVh1b_ss = 0.0_dp
869 137657068 : aVh1b_00 = 0.0_dp
870 :
871 40246 : nr = SIZE(gsph, 1)
872 :
873 40246 : m1 = 0
874 140431 : DO iset1 = 1, nset
875 : m2 = 0
876 432532 : DO iset2 = 1, nset
877 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
878 332347 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
879 :
880 332347 : n1 = nsoset(lmax(iset1))
881 1171643 : DO ipgf1 = 1, npgf(iset1)
882 839296 : n2 = nsoset(lmax(iset2))
883 3624868 : DO ipgf2 = 1, npgf(iset2)
884 : ! with contributions to V1_s*rho0
885 25620386 : DO iso = 1, MIN(nchan_0, max_iso_not0)
886 23167161 : l_ang = indso(1, iso)
887 1239345105 : gVg_0 = SUM(Vh1_s(:, iso)*g0_h_w(:, l_ang))
888 49835456 : DO icg = 1, cg_n_list(iso)
889 24215070 : is1 = cg_list(1, icg, iso)
890 24215070 : is2 = cg_list(2, icg, iso)
891 :
892 24215070 : iso1 = is1 + n1*(ipgf1 - 1) + m1
893 24215070 : iso2 = is2 + n2*(ipgf2 - 1) + m2
894 24215070 : gVg_h = 0.0_dp
895 24215070 : gVg_s = 0.0_dp
896 :
897 1286791864 : DO ir = 1, nr
898 1262576794 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
899 1286791864 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
900 : END DO ! ir
901 :
902 24215070 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
903 24215070 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
904 47382231 : aVh1b_00(iso1, iso2) = aVh1b_00(iso1, iso2) + gVg_0*Qlm_gg(iso1, iso2, iso)
905 :
906 : END DO !icg
907 : END DO ! iso
908 : ! without contributions to V1_s*rho0
909 32099873 : DO iso = nchan_0 + 1, max_iso_not0
910 42461174 : DO icg = 1, cg_n_list(iso)
911 11200597 : is1 = cg_list(1, icg, iso)
912 11200597 : is2 = cg_list(2, icg, iso)
913 :
914 11200597 : iso1 = is1 + n1*(ipgf1 - 1) + m1
915 11200597 : iso2 = is2 + n2*(ipgf2 - 1) + m2
916 11200597 : gVg_h = 0.0_dp
917 11200597 : gVg_s = 0.0_dp
918 :
919 581173787 : DO ir = 1, nr
920 569973190 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
921 581173787 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
922 : END DO ! ir
923 :
924 11200597 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
925 40007949 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
926 :
927 : END DO !icg
928 : END DO ! iso
929 : END DO ! ipgf2
930 : END DO ! ipgf1
931 432532 : m2 = m2 + maxso
932 : END DO ! iset2
933 140431 : m1 = m1 + maxso
934 : END DO !iset1
935 86096 : DO ispin = 1, nspins
936 45850 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
937 45850 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
938 45850 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_h(ispin)%r_coef, 1)
939 86096 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_s(ispin)%r_coef, 1)
940 : END DO ! ispin
941 :
942 40246 : END SUBROUTINE Vh_1c_atom_integrals
943 :
944 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
945 :
946 : END MODULE hartree_local_methods
|
