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