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 : 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 : #include "./base/base_uses.f90"
53 :
54 : IMPLICIT NONE
55 :
56 : PRIVATE
57 :
58 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hartree_local_methods'
59 :
60 : ! Public Subroutine
61 :
62 : PUBLIC :: init_coulomb_local, calculate_Vh_1center, Vh_1c_gg_integrals
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief ...
68 : !> \param hartree_local ...
69 : !> \param natom ...
70 : ! **************************************************************************************************
71 1892 : SUBROUTINE init_coulomb_local(hartree_local, natom)
72 :
73 : TYPE(hartree_local_type), POINTER :: hartree_local
74 : INTEGER, INTENT(IN) :: natom
75 :
76 : CHARACTER(len=*), PARAMETER :: routineN = 'init_coulomb_local'
77 :
78 : INTEGER :: handle
79 1892 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
80 :
81 1892 : CALL timeset(routineN, handle)
82 :
83 1892 : NULLIFY (ecoul_1c)
84 : ! Allocate and Initialize 1-center Potentials and Integrals
85 1892 : CALL allocate_ecoul_1center(ecoul_1c, natom)
86 1892 : hartree_local%ecoul_1c => ecoul_1c
87 :
88 1892 : CALL timestop(handle)
89 :
90 1892 : END SUBROUTINE init_coulomb_local
91 :
92 : ! **************************************************************************************************
93 : !> \brief Calculates Hartree potential for hard and soft densities (including
94 : !> nuclear charge and compensation charges) using numerical integration
95 : !> \param vrad_h ...
96 : !> \param vrad_s ...
97 : !> \param rrad_h ...
98 : !> \param rrad_s ...
99 : !> \param rrad_0 ...
100 : !> \param rrad_z ...
101 : !> \param grid_atom ...
102 : !> \par History
103 : !> 05.2012 JGH refactoring
104 : !> \author ??
105 : ! **************************************************************************************************
106 15 : SUBROUTINE calculate_Vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom)
107 :
108 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vrad_h, vrad_s
109 : TYPE(rho_atom_coeff), DIMENSION(:), INTENT(IN) :: rrad_h, rrad_s
110 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: rrad_0
111 : REAL(dp), DIMENSION(:), INTENT(IN) :: rrad_z
112 : TYPE(grid_atom_type), POINTER :: grid_atom
113 :
114 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_Vh_1center'
115 :
116 : INTEGER :: handle, ir, iso, ispin, l_ang, &
117 : max_s_harm, nchannels, nr, nspins
118 : REAL(dp) :: I1_down, I1_up, I2_down, I2_up, prefactor
119 15 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rho_1, rho_2
120 15 : REAL(dp), DIMENSION(:), POINTER :: wr
121 15 : REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l
122 :
123 15 : CALL timeset(routineN, handle)
124 :
125 15 : nr = grid_atom%nr
126 15 : max_s_harm = SIZE(vrad_h, 2)
127 15 : nspins = SIZE(rrad_h, 1)
128 15 : nchannels = SIZE(rrad_0, 2)
129 :
130 15 : r2l => grid_atom%rad2l
131 15 : oor2l => grid_atom%oorad2l
132 15 : wr => grid_atom%wr
133 :
134 90 : ALLOCATE (rho_1(nr, max_s_harm), rho_2(nr, max_s_harm))
135 9348 : rho_1 = 0.0_dp
136 9348 : rho_2 = 0.0_dp
137 :
138 : ! Case lm = 0
139 765 : rho_1(:, 1) = rrad_z(:)
140 765 : rho_2(:, 1) = rrad_0(:, 1)
141 :
142 135 : DO iso = 2, nchannels
143 6135 : rho_2(:, iso) = rrad_0(:, iso)
144 : END DO
145 :
146 198 : DO iso = 1, max_s_harm
147 549 : DO ispin = 1, nspins
148 18666 : rho_1(:, iso) = rho_1(:, iso) + rrad_h(ispin)%r_coef(:, iso)
149 18849 : rho_2(:, iso) = rho_2(:, iso) + rrad_s(ispin)%r_coef(:, iso)
150 : END DO
151 :
152 183 : l_ang = indso(1, iso)
153 183 : prefactor = fourpi/(2._dp*l_ang + 1._dp)
154 :
155 9333 : rho_1(:, iso) = rho_1(:, iso)*wr(:)
156 9333 : rho_2(:, iso) = rho_2(:, iso)*wr(:)
157 :
158 183 : I1_up = 0.0_dp
159 183 : I1_down = 0.0_dp
160 183 : I2_up = 0.0_dp
161 183 : I2_down = 0.0_dp
162 :
163 183 : I1_up = r2l(nr, l_ang)*rho_1(nr, iso)
164 183 : I2_up = r2l(nr, l_ang)*rho_2(nr, iso)
165 :
166 9150 : DO ir = nr - 1, 1, -1
167 8967 : I1_down = I1_down + oor2l(ir, l_ang + 1)*rho_1(ir, iso)
168 9150 : I2_down = I2_down + oor2l(ir, l_ang + 1)*rho_2(ir, iso)
169 : END DO
170 :
171 : vrad_h(nr, iso) = vrad_h(nr, iso) + prefactor* &
172 183 : (oor2l(nr, l_ang + 1)*I1_up + r2l(nr, l_ang)*I1_down)
173 : vrad_s(nr, iso) = vrad_s(nr, iso) + prefactor* &
174 183 : (oor2l(nr, l_ang + 1)*I2_up + r2l(nr, l_ang)*I2_down)
175 :
176 9165 : DO ir = nr - 1, 1, -1
177 8967 : I1_up = I1_up + r2l(ir, l_ang)*rho_1(ir, iso)
178 8967 : I1_down = I1_down - oor2l(ir, l_ang + 1)*rho_1(ir, iso)
179 8967 : I2_up = I2_up + r2l(ir, l_ang)*rho_2(ir, iso)
180 8967 : I2_down = I2_down - oor2l(ir, l_ang + 1)*rho_2(ir, iso)
181 :
182 : vrad_h(ir, iso) = vrad_h(ir, iso) + prefactor* &
183 8967 : (oor2l(ir, l_ang + 1)*I1_up + r2l(ir, l_ang)*I1_down)
184 : vrad_s(ir, iso) = vrad_s(ir, iso) + prefactor* &
185 9150 : (oor2l(ir, l_ang + 1)*I2_up + r2l(ir, l_ang)*I2_down)
186 :
187 : END DO
188 :
189 : END DO
190 :
191 15 : DEALLOCATE (rho_1, rho_2)
192 :
193 15 : CALL timestop(handle)
194 :
195 15 : END SUBROUTINE calculate_Vh_1center
196 :
197 : ! **************************************************************************************************
198 : !> \brief Calculates one center GAPW Hartree energies and matrix elements
199 : !> Hartree potentials are input
200 : !> Takes possible background charge into account
201 : !> Special case for densities without core charge
202 : !> \param qs_env ...
203 : !> \param energy_hartree_1c ...
204 : !> \param ecoul_1c ...
205 : !> \param local_rho_set ...
206 : !> \param para_env ...
207 : !> \param tddft ...
208 : !> \param local_rho_set_2nd ...
209 : !> \param core_2nd ...
210 : !> \par History
211 : !> 05.2012 JGH refactoring
212 : !> \author ??
213 : ! **************************************************************************************************
214 17104 : SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, &
215 : core_2nd)
216 :
217 : TYPE(qs_environment_type), POINTER :: qs_env
218 : REAL(kind=dp), INTENT(out) :: energy_hartree_1c
219 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
220 : TYPE(local_rho_type), POINTER :: local_rho_set
221 : TYPE(mp_para_env_type), POINTER :: para_env
222 : LOGICAL, INTENT(IN) :: tddft
223 : TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set_2nd
224 : LOGICAL, INTENT(IN), OPTIONAL :: core_2nd
225 :
226 : CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_gg_integrals'
227 :
228 : INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, &
229 : llmax_nuc, lmax0, lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_iso_not0_nuc, &
230 : max_s_harm, max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, nat, nchan_0, &
231 : nkind, nr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
232 17104 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
233 17104 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
234 17104 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmax_nuc, lmin, &
235 17104 : lmin_nuc, npgf, npgf_nuc
236 : LOGICAL :: cneo, core_charge, l_2nd_local_rho, &
237 : my_core_2nd, my_periodic, paw_atom
238 : REAL(dp) :: back_ch, ecoul_1_z_cneo, factor
239 17104 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr
240 17104 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_00_nuc, aVh1b_hh, &
241 17104 : aVh1b_hh_nuc, aVh1b_ss, aVh1b_ss_nuc, &
242 17104 : g0_h_w
243 17104 : REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
244 17104 : REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, gsph_nuc, rrad_0, &
245 17104 : Vh1_h, Vh1_s, vrrad_0, zet, zet_nuc
246 17104 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg, Qlm_gg_2nd
247 17104 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
248 : TYPE(cell_type), POINTER :: cell
249 : TYPE(cneo_potential_type), POINTER :: cneo_potential
250 : TYPE(dft_control_type), POINTER :: dft_control
251 : TYPE(grid_atom_type), POINTER :: grid_atom
252 : TYPE(gto_basis_set_type), POINTER :: basis_1c, nuc_basis
253 : TYPE(harmonics_atom_type), POINTER :: harmonics
254 : TYPE(pw_env_type), POINTER :: pw_env
255 : TYPE(pw_poisson_type), POINTER :: poisson_env
256 : TYPE(qs_charges_type), POINTER :: qs_charges
257 17104 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
258 17104 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho0_atom_set_2nd
259 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole, rho0_mpole_2nd
260 17104 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
261 : TYPE(rho_atom_type), POINTER :: rho_atom
262 17104 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
263 : TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
264 17104 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
265 :
266 17104 : CALL timeset(routineN, handle)
267 :
268 17104 : NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
269 17104 : NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
270 17104 : NULLIFY (rho0_mpole, rhoz_set)
271 17104 : NULLIFY (atom_list, grid_atom, harmonics)
272 17104 : NULLIFY (basis_1c, lmin, lmax, npgf, zet)
273 17104 : NULLIFY (gsph)
274 17104 : NULLIFY (rhoz_cneo, rhoz_cneo_set)
275 :
276 : CALL get_qs_env(qs_env=qs_env, &
277 : cell=cell, dft_control=dft_control, &
278 : atomic_kind_set=atomic_kind_set, &
279 : qs_kind_set=qs_kind_set, &
280 17104 : pw_env=pw_env, qs_charges=qs_charges)
281 :
282 17104 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
283 17104 : my_periodic = (poisson_env%method == pw_poisson_periodic)
284 :
285 17104 : back_ch = qs_charges%background*cell%deth
286 :
287 : ! rhoz_set is not accessed in TDDFT
288 : CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, &
289 17104 : rhoz_cneo_set) ! for integral space
290 :
291 : ! for forces we need a second local_rho_set
292 17104 : l_2nd_local_rho = .FALSE.
293 17104 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
294 164 : l_2nd_local_rho = .TRUE.
295 164 : NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
296 164 : CALL get_local_rho(local_rho_set_2nd, rho_atom_set_2nd, rho0_atom_set_2nd, rho0_mpole_2nd, rhoz_set=rhoz_set_2nd)
297 : END IF
298 :
299 17104 : nkind = SIZE(atomic_kind_set, 1)
300 17104 : nspins = dft_control%nspins
301 :
302 17104 : core_charge = .NOT. tddft ! for forces mixed version
303 17104 : my_core_2nd = .TRUE.
304 17104 : IF (PRESENT(core_2nd)) my_core_2nd = .NOT. core_2nd ! if my_core_2nd true, include core charge
305 :
306 : ! The aim of the following code was to return immediately if the subroutine
307 : ! was called for triplet excited states in spin-restricted case. This check
308 : ! is also performed before invocation of this subroutine. It should be save
309 : ! to remove the optional argument 'do_triplet' from the subroutine interface.
310 : !IF (tddft) THEN
311 : ! CPASSERT(PRESENT(do_triplet))
312 : ! IF (nspins == 1 .AND. do_triplet) RETURN
313 : !END IF
314 :
315 17104 : CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
316 17104 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax_0)
317 :
318 : ! Put to 0 the local hartree energy contribution from 1 center integrals
319 17104 : energy_hartree_1c = 0.0_dp
320 : ! Restore total quantum nuclear density to zero
321 17104 : qs_charges%total_rho1_hard_nuc = 0.0_dp
322 17104 : rho0_mpole%tot_rhoz_cneo_s = 0.0_dp
323 :
324 : ! Here starts the loop over all the atoms
325 50540 : DO ikind = 1, nkind
326 :
327 33436 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
328 33436 : NULLIFY (cneo_potential)
329 : CALL get_qs_kind(qs_kind_set(ikind), &
330 : grid_atom=grid_atom, &
331 : harmonics=harmonics, ngrid_rad=nr, &
332 : max_iso_not0=max_iso_not0, paw_atom=paw_atom, &
333 33436 : cneo_potential=cneo_potential)
334 : CALL get_qs_kind(qs_kind_set(ikind), &
335 33436 : basis_set=basis_1c, basis_type="GAPW_1C")
336 :
337 33436 : cneo = ASSOCIATED(cneo_potential)
338 33436 : IF (cneo .AND. tddft) &
339 0 : CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.")
340 :
341 33436 : NULLIFY (nuc_basis)
342 33436 : max_iso_not0_nuc = 0
343 33436 : IF (cneo) THEN
344 48 : CPASSERT(paw_atom)
345 : CALL get_qs_kind(qs_kind_set(ikind), &
346 48 : basis_set=nuc_basis, basis_type="NUC")
347 48 : max_iso_not0_nuc = cneo_potential%harmonics%max_iso_not0
348 : END IF
349 :
350 83976 : IF (paw_atom) THEN
351 : !=========== PAW ===============
352 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
353 : maxso=maxso, npgf=npgf, maxl=maxl, &
354 31846 : nset=nset, zet=zet)
355 :
356 31846 : max_s_harm = harmonics%max_s_harm
357 31846 : llmax = harmonics%llmax
358 :
359 31846 : nsotot = maxso*nset
360 127384 : ALLOCATE (gsph(nr, nsotot))
361 95538 : ALLOCATE (gexp(nr))
362 159230 : ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0))
363 :
364 127384 : ALLOCATE (aVh1b_hh(nsotot, nsotot))
365 95538 : ALLOCATE (aVh1b_ss(nsotot, nsotot))
366 95538 : ALLOCATE (aVh1b_00(nsotot, nsotot))
367 :
368 31846 : NULLIFY (Qlm_gg, g0_h)
369 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
370 : l0_ikind=lmax0, &
371 31846 : Qlm_gg=Qlm_gg, g0_h=g0_h) ! Qlm_gg of density
372 :
373 31846 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
374 332 : NULLIFY (Qlm_gg_2nd, g0_h_2nd)
375 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole_2nd, ikind=ikind, &
376 : l0_ikind=lmax0_2nd, &
377 332 : Qlm_gg=Qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
378 : END IF
379 31846 : nchan_0 = nsoset(lmax0)
380 :
381 31846 : IF (nchan_0 > MAX(max_iso_not0, max_iso_not0_nuc)) &
382 0 : CPABORT("channels for rho0 > # max of spherical harmonics")
383 :
384 31846 : NULLIFY (Vh1_h, Vh1_s)
385 127384 : ALLOCATE (Vh1_h(nr, max_iso_not0))
386 127384 : ALLOCATE (Vh1_s(nr, MAX(max_iso_not0, max_iso_not0_nuc, nchan_0)))
387 :
388 31846 : NULLIFY (lmax_nuc, lmin_nuc, npgf_nuc, zet_nuc, gsph_nuc)
389 31846 : maxso_nuc = 0
390 31846 : maxl_nuc = -1
391 31846 : nset_nuc = 0
392 31846 : max_s_harm_nuc = 0
393 31846 : llmax_nuc = -1
394 31846 : nsotot_nuc = 0
395 31846 : IF (cneo) THEN
396 : CALL get_gto_basis_set(gto_basis_set=nuc_basis, lmax=lmax_nuc, &
397 : lmin=lmin_nuc, maxso=maxso_nuc, npgf=npgf_nuc, &
398 48 : maxl=maxl_nuc, nset=nset_nuc, zet=zet_nuc)
399 :
400 48 : max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
401 48 : llmax_nuc = cneo_potential%harmonics%llmax
402 48 : nsotot_nuc = maxso_nuc*nset_nuc
403 192 : ALLOCATE (gsph_nuc(nr, nsotot_nuc))
404 198336 : gsph_nuc = 0.0_dp
405 192 : ALLOCATE (aVh1b_hh_nuc(nsotot_nuc, nsotot_nuc))
406 144 : ALLOCATE (aVh1b_ss_nuc(nsotot_nuc, nsotot_nuc))
407 192 : ALLOCATE (aVh1b_00_nuc(nsotot_nuc, nsotot_nuc))
408 : END IF
409 :
410 0 : ALLOCATE (cg_list(2, nsoset(MAX(maxl, maxl_nuc))**2, &
411 : MAX(max_s_harm, max_s_harm_nuc)), &
412 191076 : cg_n_list(MAX(max_s_harm, max_s_harm_nuc)))
413 :
414 31846 : NULLIFY (rrad_z, my_CG)
415 31846 : my_CG => harmonics%my_CG
416 :
417 : ! set to zero temporary arrays
418 1851466 : sqrtwr = 0.0_dp
419 5555348 : g0_h_w = 0.0_dp
420 1851466 : gexp = 0.0_dp
421 60431924 : gsph = 0.0_dp
422 :
423 1851466 : sqrtwr(1:nr) = SQRT(grid_atom%wr(1:nr))
424 117816 : DO l_ang = 0, lmax0
425 4952516 : g0_h_w(1:nr, l_ang) = g0_h(1:nr, l_ang)*grid_atom%wr(1:nr)
426 : END DO
427 :
428 : m1 = 0
429 115724 : DO iset1 = 1, nset
430 83878 : n1 = nsoset(lmax(iset1))
431 304942 : DO ipgf1 = 1, npgf(iset1)
432 12795964 : gexp(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
433 885628 : DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1))
434 580686 : iso = is1 + (ipgf1 - 1)*n1 + m1
435 580686 : l_ang = indso(1, is1)
436 66281030 : gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
437 : END DO ! is1
438 : END DO ! ipgf1
439 115724 : m1 = m1 + maxso
440 : END DO ! iset1
441 :
442 31846 : IF (cneo) THEN
443 : ! initialize nuclear pmat, cpc and e_core to zero
444 140 : DO iat = 1, nat
445 92 : iatom = atom_list(iat)
446 50876 : rhoz_cneo_set(iatom)%pmat = 0.0_dp
447 50876 : rhoz_cneo_set(iatom)%cpc_h = 0.0_dp
448 50876 : rhoz_cneo_set(iatom)%cpc_s = 0.0_dp
449 92 : rhoz_cneo_set(iatom)%e_core = 0.0_dp
450 140 : rhoz_cneo_set(iatom)%ready = .FALSE.
451 : END DO
452 :
453 : ! calculate nuclear gsph
454 48 : m1 = 0
455 480 : DO iset1 = 1, nset_nuc
456 432 : n1 = nsoset(lmax_nuc(iset1))
457 864 : DO ipgf1 = 1, npgf_nuc(iset1)
458 22032 : gexp(1:nr) = EXP(-zet_nuc(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
459 1968 : DO is1 = nsoset(lmin_nuc(iset1) - 1) + 1, nsoset(lmax_nuc(iset1))
460 1104 : iso = is1 + (ipgf1 - 1)*n1 + m1
461 1104 : l_ang = indso(1, is1)
462 111936 : gsph_nuc(1:nr, iso) = cneo_potential%rad2l(1:nr, l_ang)*gexp(1:nr)
463 : END DO ! is1
464 : END DO ! ipgf1
465 480 : m1 = m1 + maxso_nuc
466 : END DO ! iset1
467 : END IF
468 :
469 : ! Distribute the atoms of this kind
470 31846 : num_pe = para_env%num_pe
471 31846 : mepos = para_env%mepos
472 31846 : bo = get_limit(nat, num_pe, mepos)
473 :
474 57734 : DO iat = bo(1), bo(2) !1,nat
475 25888 : iatom = atom_list(iat)
476 25888 : rho_atom => rho_atom_set(iatom)
477 :
478 25888 : NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
479 25888 : IF (core_charge .AND. .NOT. cneo) THEN
480 22381 : rrad_z => rhoz_set(ikind)%r_coef ! for density
481 : END IF
482 25888 : IF (my_core_2nd .AND. .NOT. cneo) THEN
483 22484 : IF (l_2nd_local_rho) THEN
484 124 : vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
485 : ELSE
486 22360 : vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
487 : END IF
488 : END IF
489 25888 : rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
490 25888 : vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
491 25888 : IF (l_2nd_local_rho) THEN
492 248 : rho_atom => rho_atom_set_2nd(iatom)
493 248 : vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
494 : END IF
495 25888 : IF (my_periodic .AND. back_ch > 1.E-3_dp) THEN
496 896 : factor = -2.0_dp*pi/3.0_dp*SQRT(fourpi)*qs_charges%background
497 : ELSE
498 24992 : factor = 0._dp
499 : END IF
500 :
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 25888 : nchan_0, nspins, max_iso_not0, factor)
505 :
506 25888 : IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom) ! rho_atom for density
507 :
508 25888 : ecoul_1_z_cneo = 0.0_dp
509 25888 : IF (cneo) THEN
510 46 : rhoz_cneo => rhoz_cneo_set(iatom)
511 : ! Add the soft tail of nuclear Hartree potential to total Vh1_s first.
512 : ! vrho already contains the -zeff factor.
513 1196 : DO iso = 1, max_iso_not0_nuc
514 116196 : Vh1_s(:, iso) = Vh1_s(:, iso) + rhoz_cneo%vrho_rad_s(:, iso)
515 : END DO
516 :
517 : ! Build nuclear 1c integrals according to Vh1_h from electronic density only
518 : ! and Vh1_s from electron density and last step (or initial guess) nuclear density
519 : ! Vh1_h_nuc = -Z*Vh1_h, Vh1_s_nuc = -Z*Vh1_s
520 : CALL Vh_1c_nuc_integrals(rhoz_cneo, cneo_potential%zeff, &
521 : aVh1b_hh_nuc, aVh1b_ss_nuc, aVh1b_00_nuc, Vh1_h, Vh1_s, &
522 : max_iso_not0, max_iso_not0_nuc, &
523 : max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
524 : nset_nuc, npgf_nuc, lmin_nuc, lmax_nuc, nsotot_nuc, maxso_nuc, &
525 : nchan_0, gsph_nuc, g0_h_w, cneo_potential%harmonics%my_CG, &
526 46 : cneo_potential%Qlm_gg)
527 :
528 : ! Solve the nuclear 1c problem
529 : CALL calculate_rhoz_cneo(rhoz_cneo, cneo_potential, cg_list, cg_n_list, nset_nuc, &
530 46 : npgf_nuc, lmin_nuc, lmax_nuc, maxl_nuc, maxso_nuc)
531 : ! slm_int(iso=1) = sqrt(4*Pi), slm_int(iso>1) = 0, without using Lebedev grid
532 : ! when printing, nuclear density is positive, thus needing the minus sign
533 : qs_charges%total_rho1_hard_nuc = qs_charges%total_rho1_hard_nuc - SQRT(fourpi) &
534 2346 : *SUM(rhoz_cneo%rho_rad_h(:, 1)*grid_atom%wr(:))
535 : rho0_mpole%tot_rhoz_cneo_s = rho0_mpole%tot_rhoz_cneo_s - SQRT(fourpi) &
536 2346 : *SUM(rhoz_cneo%rho_rad_s(:, 1)*grid_atom%wr(:))
537 :
538 : ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz_cneo_h and Vh1_s*rhoz_cneo_s.
539 : ! Self-interaction of the quantum nucleus is already removed in ecoul_1_z_cneo.
540 : ! rho already contains the -zeff factor.
541 460 : DO iso = 1, MIN(max_iso_not0, max_iso_not0_nuc)
542 : ecoul_1_z_cneo = ecoul_1_z_cneo + 0.5_dp* &
543 : (SUM(Vh1_h(:, iso)*rhoz_cneo%rho_rad_h(:, iso)*grid_atom%wr(:)) &
544 41860 : - SUM(Vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:)))
545 : END DO
546 782 : DO iso = max_iso_not0 + 1, max_iso_not0_nuc
547 : ecoul_1_z_cneo = ecoul_1_z_cneo - 0.5_dp* &
548 37582 : SUM(Vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:))
549 : END DO
550 :
551 : ! Add nuclear Hartree potential to total Vh1_h after solving the nuclear 1c problem
552 : ! to avoid nuclear self-interaction.
553 : ! vrho already contains the -zeff factor.
554 : ! Here the min of two max_iso_not0's is chosen, because even when the nuclear one
555 : ! is larger, it is meaningless to let Vh1_h have higher angular momentum components
556 : ! as Vh1_h now is only used for the electronic part.
557 460 : DO iso = 1, MIN(max_iso_not0, max_iso_not0_nuc)
558 41860 : Vh1_h(:, iso) = Vh1_h(:, iso) + rhoz_cneo%vrho_rad_h(:, iso)
559 : END DO
560 : END IF
561 :
562 : CALL Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
563 : grid_atom, iatom, core_charge .AND. .NOT. cneo, & ! core charge for density
564 25888 : rrad_z, Vh1_h, Vh1_s, nchan_0, nspins, max_iso_not0)
565 :
566 25888 : IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)
567 :
568 25888 : IF (cneo) THEN
569 46 : CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1c(iatom)%ecoul_1_z + ecoul_1_z_cneo)
570 46 : energy_hartree_1c = energy_hartree_1c + ecoul_1_z_cneo
571 : END IF
572 :
573 : CALL Vh_1c_atom_integrals(rho_atom, & ! results (int_local_h and int_local_s) written on rho_atom_2nd
574 : ! int_local_h and int_local_s are used in update_ks_atom
575 : ! on int_local_h mixed core / non-core
576 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
577 : max_s_harm, llmax, cg_list, cg_n_list, &
578 : nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
579 57734 : g0_h_w, my_CG, Qlm_gg) ! Qlm_gg for density from local_rho_set
580 :
581 : END DO ! iat
582 :
583 31846 : DEALLOCATE (aVh1b_hh)
584 31846 : DEALLOCATE (aVh1b_ss)
585 31846 : DEALLOCATE (aVh1b_00)
586 31846 : IF (ALLOCATED(aVh1b_hh_nuc)) DEALLOCATE (aVh1b_hh_nuc)
587 31846 : IF (ALLOCATED(aVh1b_ss_nuc)) DEALLOCATE (aVh1b_ss_nuc)
588 31846 : IF (ALLOCATED(aVh1b_00_nuc)) DEALLOCATE (aVh1b_00_nuc)
589 31846 : DEALLOCATE (Vh1_h, Vh1_s)
590 31846 : DEALLOCATE (cg_list, cg_n_list)
591 31846 : DEALLOCATE (gsph)
592 31846 : IF (ASSOCIATED(gsph_nuc)) DEALLOCATE (gsph_nuc)
593 31846 : DEALLOCATE (gexp)
594 31846 : DEALLOCATE (sqrtwr, g0_h_w)
595 :
596 95538 : IF (cneo) THEN
597 : ! broadcast nuclear pmat, cpc and e_core
598 140 : DO iat = 1, nat
599 92 : iatom = atom_list(iat)
600 101660 : CALL para_env%sum(rhoz_cneo_set(iatom)%pmat)
601 101660 : CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_h)
602 101660 : CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_s)
603 92 : CALL para_env%sum(rhoz_cneo_set(iatom)%e_core)
604 140 : rhoz_cneo_set(iatom)%ready = .TRUE.
605 : END DO
606 : END IF
607 : ELSE
608 : !=========== NO PAW ===============
609 : ! This term is taken care of using the core density as in GPW
610 : CYCLE
611 : END IF ! paw
612 : END DO ! ikind
613 :
614 17104 : CALL para_env%sum(energy_hartree_1c)
615 17104 : CALL para_env%sum(qs_charges%total_rho1_hard_nuc)
616 17104 : CALL para_env%sum(rho0_mpole%tot_rhoz_cneo_s)
617 :
618 17104 : CALL timestop(handle)
619 :
620 34208 : END SUBROUTINE Vh_1c_gg_integrals
621 :
622 : ! **************************************************************************************************
623 :
624 : ! **************************************************************************************************
625 : !> \brief ...
626 : !> \param rho_atom ...
627 : !> \param vrrad_0 ...
628 : !> \param grid_atom ...
629 : !> \param core_charge ...
630 : !> \param vrrad_z ...
631 : !> \param Vh1_h ...
632 : !> \param Vh1_s ...
633 : !> \param nchan_0 ...
634 : !> \param nspins ...
635 : !> \param max_iso_not0 ...
636 : !> \param bfactor ...
637 : ! **************************************************************************************************
638 25888 : SUBROUTINE Vh_1c_atom_potential(rho_atom, vrrad_0, &
639 : grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, &
640 : nchan_0, nspins, max_iso_not0, bfactor)
641 :
642 : TYPE(rho_atom_type), POINTER :: rho_atom
643 : REAL(dp), DIMENSION(:, :), POINTER :: vrrad_0
644 : TYPE(grid_atom_type), POINTER :: grid_atom
645 : LOGICAL, INTENT(IN) :: core_charge
646 : REAL(dp), DIMENSION(:), POINTER :: vrrad_z
647 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
648 : INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
649 : REAL(dp), INTENT(IN) :: bfactor
650 :
651 : INTEGER :: ir, iso, ispin, nr
652 25888 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s
653 :
654 25888 : nr = grid_atom%nr
655 :
656 25888 : NULLIFY (vr_h, vr_s)
657 25888 : CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
658 :
659 18602668 : Vh1_h = 0.0_dp
660 18640204 : Vh1_s = 0.0_dp
661 :
662 2615052 : IF (core_charge) Vh1_h(:, 1) = vrrad_z(:)
663 :
664 218376 : DO iso = 1, nchan_0
665 21277784 : Vh1_s(:, iso) = vrrad_0(:, iso)
666 : END DO
667 :
668 56373 : DO ispin = 1, nspins
669 458882 : DO iso = 1, max_iso_not0
670 23049579 : Vh1_h(:, iso) = Vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
671 23080064 : Vh1_s(:, iso) = Vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
672 : END DO
673 : END DO
674 :
675 25888 : IF (bfactor /= 0._dp) THEN
676 53896 : DO ir = 1, nr
677 53000 : Vh1_h(ir, 1) = Vh1_h(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
678 53896 : Vh1_s(ir, 1) = Vh1_s(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
679 : END DO
680 : END IF
681 :
682 25888 : END SUBROUTINE Vh_1c_atom_potential
683 :
684 : ! **************************************************************************************************
685 :
686 : ! **************************************************************************************************
687 : !> \brief ...
688 : !> \param energy_hartree_1c ...
689 : !> \param ecoul_1c ...
690 : !> \param rho_atom ...
691 : !> \param rrad_0 ...
692 : !> \param grid_atom ...
693 : !> \param iatom ...
694 : !> \param core_charge ...
695 : !> \param rrad_z ...
696 : !> \param Vh1_h ...
697 : !> \param Vh1_s ...
698 : !> \param nchan_0 ...
699 : !> \param nspins ...
700 : !> \param max_iso_not0 ...
701 : ! **************************************************************************************************
702 25888 : SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
703 : grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
704 : nchan_0, nspins, max_iso_not0)
705 :
706 : REAL(dp), INTENT(INOUT) :: energy_hartree_1c
707 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
708 : TYPE(rho_atom_type), POINTER :: rho_atom
709 : REAL(dp), DIMENSION(:, :), POINTER :: rrad_0
710 : TYPE(grid_atom_type), POINTER :: grid_atom
711 : INTEGER, INTENT(IN) :: iatom
712 : LOGICAL, INTENT(IN) :: core_charge
713 : REAL(dp), DIMENSION(:), POINTER :: rrad_z
714 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
715 : INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
716 :
717 : INTEGER :: iso, ispin, nr
718 : REAL(dp) :: ecoul_1_0, ecoul_1_h, ecoul_1_s, &
719 : ecoul_1_z
720 25888 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
721 :
722 25888 : nr = grid_atom%nr
723 :
724 25888 : NULLIFY (r_h, r_s)
725 25888 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
726 :
727 : ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz
728 25888 : ecoul_1_z = 0.0_dp
729 25888 : IF (core_charge) THEN
730 1300571 : ecoul_1_z = 0.5_dp*SUM(Vh1_h(:, 1)*rrad_z(:)*grid_atom%wr(:))
731 : END IF
732 :
733 : ! Calculate the contributions to Ecoul coming from Vh1_s*rho0
734 25888 : ecoul_1_0 = 0.0_dp
735 218376 : DO iso = 1, nchan_0
736 10651836 : ecoul_1_0 = ecoul_1_0 + 0.5_dp*SUM(Vh1_s(:, iso)*rrad_0(:, iso)*grid_atom%wr(:))
737 : END DO
738 :
739 : ! Calculate the contributions to Ecoul coming from Vh1_h*rho1_h and Vh1_s*rho1_s
740 25888 : ecoul_1_s = 0.0_dp
741 25888 : ecoul_1_h = 0.0_dp
742 56373 : DO ispin = 1, nspins
743 458882 : DO iso = 1, max_iso_not0
744 23049579 : ecoul_1_s = ecoul_1_s + 0.5_dp*SUM(Vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
745 23080064 : ecoul_1_h = ecoul_1_h + 0.5_dp*SUM(Vh1_h(:, iso)*r_h(ispin)%r_coef(:, iso)*grid_atom%wr(:))
746 : END DO
747 : END DO
748 :
749 25888 : CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
750 25888 : CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
751 :
752 25888 : energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
753 25888 : energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
754 :
755 25888 : END SUBROUTINE Vh_1c_atom_energy
756 :
757 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
758 :
759 : ! **************************************************************************************************
760 : !> \brief ...
761 : !> \param rho_atom ...
762 : !> \param aVh1b_hh ...
763 : !> \param aVh1b_ss ...
764 : !> \param aVh1b_00 ...
765 : !> \param Vh1_h ...
766 : !> \param Vh1_s ...
767 : !> \param max_iso_not0 ...
768 : !> \param max_s_harm ...
769 : !> \param llmax ...
770 : !> \param cg_list ...
771 : !> \param cg_n_list ...
772 : !> \param nset ...
773 : !> \param npgf ...
774 : !> \param lmin ...
775 : !> \param lmax ...
776 : !> \param nsotot ...
777 : !> \param maxso ...
778 : !> \param nspins ...
779 : !> \param nchan_0 ...
780 : !> \param gsph ...
781 : !> \param g0_h_w ...
782 : !> \param my_CG ...
783 : !> \param Qlm_gg ...
784 : ! **************************************************************************************************
785 25888 : SUBROUTINE Vh_1c_atom_integrals(rho_atom, &
786 25888 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
787 25888 : max_s_harm, llmax, cg_list, cg_n_list, &
788 : nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
789 25888 : g0_h_w, my_CG, Qlm_gg)
790 :
791 : TYPE(rho_atom_type), POINTER :: rho_atom
792 : REAL(dp), DIMENSION(:, :) :: aVh1b_hh, aVh1b_ss, aVh1b_00
793 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
794 : INTEGER, INTENT(IN) :: max_iso_not0, max_s_harm, llmax
795 : INTEGER, DIMENSION(:, :, :) :: cg_list
796 : INTEGER, DIMENSION(:) :: cg_n_list
797 : INTEGER, INTENT(IN) :: nset
798 : INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax
799 : INTEGER, INTENT(IN) :: nsotot, maxso, nspins, nchan_0
800 : REAL(dp), DIMENSION(:, :), POINTER :: gsph
801 : REAL(dp), DIMENSION(:, 0:) :: g0_h_w
802 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg
803 :
804 : INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, &
805 : iset2, iso, iso1, iso2, ispin, l_ang, &
806 : m1, m2, max_iso_not0_local, n1, n2, nr
807 : REAL(dp) :: gVg_0, gVg_h, gVg_s
808 25888 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
809 :
810 25888 : NULLIFY (int_local_h, int_local_s)
811 : CALL get_rho_atom(rho_atom=rho_atom, &
812 : ga_Vlocal_gb_h=int_local_h, &
813 25888 : ga_Vlocal_gb_s=int_local_s)
814 :
815 : ! Calculate the integrals of the potential with 2 primitives
816 75190324 : aVh1b_hh = 0.0_dp
817 75190324 : aVh1b_ss = 0.0_dp
818 75190324 : aVh1b_00 = 0.0_dp
819 :
820 25888 : nr = SIZE(gsph, 1)
821 :
822 25888 : m1 = 0
823 90953 : DO iset1 = 1, nset
824 : m2 = 0
825 282618 : DO iset2 = 1, nset
826 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
827 217553 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
828 :
829 217553 : n1 = nsoset(lmax(iset1))
830 739586 : DO ipgf1 = 1, npgf(iset1)
831 522033 : n2 = nsoset(lmax(iset2))
832 2217342 : DO ipgf2 = 1, npgf(iset2)
833 : ! with contributions to V1_s*rho0
834 14523544 : DO iso = 1, MIN(nchan_0, max_iso_not0)
835 13045788 : l_ang = indso(1, iso)
836 721928798 : gVg_0 = SUM(Vh1_s(:, iso)*g0_h_w(:, l_ang))
837 26946020 : DO icg = 1, cg_n_list(iso)
838 12422476 : is1 = cg_list(1, icg, iso)
839 12422476 : is2 = cg_list(2, icg, iso)
840 :
841 12422476 : iso1 = is1 + n1*(ipgf1 - 1) + m1
842 12422476 : iso2 = is2 + n2*(ipgf2 - 1) + m2
843 12422476 : gVg_h = 0.0_dp
844 12422476 : gVg_s = 0.0_dp
845 :
846 683459236 : DO ir = 1, nr
847 671036760 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
848 683459236 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
849 : END DO ! ir
850 :
851 12422476 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
852 12422476 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
853 25468264 : aVh1b_00(iso1, iso2) = aVh1b_00(iso1, iso2) + gVg_0*Qlm_gg(iso1, iso2, iso)
854 :
855 : END DO !icg
856 : END DO ! iso
857 : ! without contributions to V1_s*rho0
858 21584661 : DO iso = nchan_0 + 1, max_iso_not0
859 27229661 : DO icg = 1, cg_n_list(iso)
860 6167033 : is1 = cg_list(1, icg, iso)
861 6167033 : is2 = cg_list(2, icg, iso)
862 :
863 6167033 : iso1 = is1 + n1*(ipgf1 - 1) + m1
864 6167033 : iso2 = is2 + n2*(ipgf2 - 1) + m2
865 6167033 : gVg_h = 0.0_dp
866 6167033 : gVg_s = 0.0_dp
867 :
868 323427913 : DO ir = 1, nr
869 317260880 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
870 323427913 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
871 : END DO ! ir
872 :
873 6167033 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
874 25751905 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
875 :
876 : END DO !icg
877 : END DO ! iso
878 : END DO ! ipgf2
879 : END DO ! ipgf1
880 282618 : m2 = m2 + maxso
881 : END DO ! iset2
882 90953 : m1 = m1 + maxso
883 : END DO !iset1
884 56373 : DO ispin = 1, nspins
885 30485 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
886 30485 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
887 30485 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_h(ispin)%r_coef, 1)
888 56373 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_s(ispin)%r_coef, 1)
889 : END DO ! ispin
890 :
891 25888 : END SUBROUTINE Vh_1c_atom_integrals
892 :
893 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
894 :
895 : END MODULE hartree_local_methods
|
