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_environment_types, ONLY: get_qs_env,&
31 : qs_environment_type
32 : USE qs_grid_atom, ONLY: grid_atom_type
33 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
34 : harmonics_atom_type
35 : USE qs_kind_types, ONLY: get_qs_kind,&
36 : get_qs_kind_set,&
37 : qs_kind_type
38 : USE qs_local_rho_types, ONLY: get_local_rho,&
39 : local_rho_type,&
40 : rhoz_type
41 : USE qs_rho0_types, ONLY: get_rho0_mpole,&
42 : rho0_atom_type,&
43 : rho0_mpole_type
44 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
45 : rho_atom_coeff,&
46 : rho_atom_type
47 : USE util, ONLY: get_limit
48 : #include "./base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 :
52 : PRIVATE
53 :
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hartree_local_methods'
55 :
56 : ! Public Subroutine
57 :
58 : PUBLIC :: init_coulomb_local, calculate_Vh_1center, Vh_1c_gg_integrals
59 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief ...
64 : !> \param hartree_local ...
65 : !> \param natom ...
66 : ! **************************************************************************************************
67 1664 : SUBROUTINE init_coulomb_local(hartree_local, natom)
68 :
69 : TYPE(hartree_local_type), POINTER :: hartree_local
70 : INTEGER, INTENT(IN) :: natom
71 :
72 : CHARACTER(len=*), PARAMETER :: routineN = 'init_coulomb_local'
73 :
74 : INTEGER :: handle
75 1664 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
76 :
77 1664 : CALL timeset(routineN, handle)
78 :
79 1664 : NULLIFY (ecoul_1c)
80 : ! Allocate and Initialize 1-center Potentials and Integrals
81 1664 : CALL allocate_ecoul_1center(ecoul_1c, natom)
82 1664 : hartree_local%ecoul_1c => ecoul_1c
83 :
84 1664 : CALL timestop(handle)
85 :
86 1664 : END SUBROUTINE init_coulomb_local
87 :
88 : ! **************************************************************************************************
89 : !> \brief Calculates Hartree potential for hard and soft densities (including
90 : !> nuclear charge and compensation charges) using numerical integration
91 : !> \param vrad_h ...
92 : !> \param vrad_s ...
93 : !> \param rrad_h ...
94 : !> \param rrad_s ...
95 : !> \param rrad_0 ...
96 : !> \param rrad_z ...
97 : !> \param grid_atom ...
98 : !> \par History
99 : !> 05.2012 JGH refactoring
100 : !> \author ??
101 : ! **************************************************************************************************
102 15 : SUBROUTINE calculate_Vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom)
103 :
104 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vrad_h, vrad_s
105 : TYPE(rho_atom_coeff), DIMENSION(:), INTENT(IN) :: rrad_h, rrad_s
106 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: rrad_0
107 : REAL(dp), DIMENSION(:), INTENT(IN) :: rrad_z
108 : TYPE(grid_atom_type), POINTER :: grid_atom
109 :
110 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_Vh_1center'
111 :
112 : INTEGER :: handle, ir, iso, ispin, l_ang, &
113 : max_s_harm, nchannels, nr, nspins
114 : REAL(dp) :: I1_down, I1_up, I2_down, I2_up, prefactor
115 15 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rho_1, rho_2
116 15 : REAL(dp), DIMENSION(:), POINTER :: wr
117 15 : REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l
118 :
119 15 : CALL timeset(routineN, handle)
120 :
121 15 : nr = grid_atom%nr
122 15 : max_s_harm = SIZE(vrad_h, 2)
123 15 : nspins = SIZE(rrad_h, 1)
124 15 : nchannels = SIZE(rrad_0, 2)
125 :
126 15 : r2l => grid_atom%rad2l
127 15 : oor2l => grid_atom%oorad2l
128 15 : wr => grid_atom%wr
129 :
130 90 : ALLOCATE (rho_1(nr, max_s_harm), rho_2(nr, max_s_harm))
131 9348 : rho_1 = 0.0_dp
132 9348 : rho_2 = 0.0_dp
133 :
134 : ! Case lm = 0
135 765 : rho_1(:, 1) = rrad_z(:)
136 765 : rho_2(:, 1) = rrad_0(:, 1)
137 :
138 135 : DO iso = 2, nchannels
139 6135 : rho_2(:, iso) = rrad_0(:, iso)
140 : END DO
141 :
142 198 : DO iso = 1, max_s_harm
143 549 : DO ispin = 1, nspins
144 18666 : rho_1(:, iso) = rho_1(:, iso) + rrad_h(ispin)%r_coef(:, iso)
145 18849 : rho_2(:, iso) = rho_2(:, iso) + 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(:, iso) = rho_1(:, iso)*wr(:)
152 9333 : rho_2(:, iso) = rho_2(:, iso)*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, iso)
160 183 : I2_up = r2l(nr, l_ang)*rho_2(nr, iso)
161 :
162 9150 : DO ir = nr - 1, 1, -1
163 8967 : I1_down = I1_down + oor2l(ir, l_ang + 1)*rho_1(ir, iso)
164 9150 : I2_down = I2_down + oor2l(ir, l_ang + 1)*rho_2(ir, iso)
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, iso)
174 8967 : I1_down = I1_down - oor2l(ir, l_ang + 1)*rho_1(ir, iso)
175 8967 : I2_up = I2_up + r2l(ir, l_ang)*rho_2(ir, iso)
176 8967 : I2_down = I2_down - oor2l(ir, l_ang + 1)*rho_2(ir, iso)
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 15744 : 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, lmax0, &
225 : lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_s_harm, maxl, maxso, mepos, n1, nat, &
226 : nchan_0, nkind, nr, nset, nsotot, nspins, num_pe
227 15744 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
228 15744 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
229 15744 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf
230 : LOGICAL :: core_charge, l_2nd_local_rho, &
231 : my_core_2nd, my_periodic, paw_atom
232 : REAL(dp) :: back_ch, factor
233 15744 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr
234 15744 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_hh, aVh1b_ss, g0_h_w
235 15744 : REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
236 15744 : REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, rrad_0, Vh1_h, &
237 15744 : Vh1_s, vrrad_0, zet
238 15744 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg, Qlm_gg_2nd
239 15744 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
240 : TYPE(cell_type), POINTER :: cell
241 : TYPE(dft_control_type), POINTER :: dft_control
242 : TYPE(grid_atom_type), POINTER :: grid_atom
243 : TYPE(gto_basis_set_type), POINTER :: basis_1c
244 : TYPE(harmonics_atom_type), POINTER :: harmonics
245 : TYPE(pw_env_type), POINTER :: pw_env
246 : TYPE(pw_poisson_type), POINTER :: poisson_env
247 : TYPE(qs_charges_type), POINTER :: qs_charges
248 15744 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
249 15744 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho0_atom_set_2nd
250 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole, rho0_mpole_2nd
251 15744 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
252 : TYPE(rho_atom_type), POINTER :: rho_atom
253 15744 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
254 :
255 15744 : CALL timeset(routineN, handle)
256 :
257 15744 : NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
258 15744 : NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
259 15744 : NULLIFY (rho0_mpole, rhoz_set)
260 15744 : NULLIFY (atom_list, grid_atom, harmonics)
261 15744 : NULLIFY (basis_1c, lmin, lmax, npgf, zet)
262 15744 : NULLIFY (gsph)
263 :
264 : CALL get_qs_env(qs_env=qs_env, &
265 : cell=cell, dft_control=dft_control, &
266 : atomic_kind_set=atomic_kind_set, &
267 : qs_kind_set=qs_kind_set, &
268 15744 : pw_env=pw_env, qs_charges=qs_charges)
269 :
270 15744 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
271 15744 : my_periodic = (poisson_env%method == pw_poisson_periodic)
272 :
273 15744 : back_ch = qs_charges%background*cell%deth
274 :
275 : ! rhoz_set is not accessed in TDDFT
276 15744 : CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set) ! for integral space
277 :
278 : ! for forces we need a second local_rho_set
279 15744 : l_2nd_local_rho = .FALSE.
280 15744 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
281 128 : l_2nd_local_rho = .TRUE.
282 128 : NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
283 128 : CALL get_local_rho(local_rho_set_2nd, rho_atom_set_2nd, rho0_atom_set_2nd, rho0_mpole_2nd, rhoz_set=rhoz_set_2nd)
284 : END IF
285 :
286 15744 : nkind = SIZE(atomic_kind_set, 1)
287 15744 : nspins = dft_control%nspins
288 :
289 15744 : core_charge = .NOT. tddft ! for forces mixed version
290 15744 : my_core_2nd = .TRUE.
291 15744 : IF (PRESENT(core_2nd)) my_core_2nd = .NOT. core_2nd ! if my_core_2nd true, include core charge
292 :
293 : ! The aim of the following code was to return immediately if the subroutine
294 : ! was called for triplet excited states in spin-restricted case. This check
295 : ! is also performed before invocation of this subroutine. It should be save
296 : ! to remove the optional argument 'do_triplet' from the subroutine interface.
297 : !IF (tddft) THEN
298 : ! CPASSERT(PRESENT(do_triplet))
299 : ! IF (nspins == 1 .AND. do_triplet) RETURN
300 : !END IF
301 :
302 15744 : CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
303 15744 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax_0)
304 :
305 : ! Put to 0 the local hartree energy contribution from 1 center integrals
306 15744 : energy_hartree_1c = 0.0_dp
307 :
308 : ! Here starts the loop over all the atoms
309 47110 : DO ikind = 1, nkind
310 :
311 31366 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
312 : CALL get_qs_kind(qs_kind_set(ikind), &
313 : grid_atom=grid_atom, &
314 : harmonics=harmonics, ngrid_rad=nr, &
315 31366 : max_iso_not0=max_iso_not0, paw_atom=paw_atom)
316 : CALL get_qs_kind(qs_kind_set(ikind), &
317 31366 : basis_set=basis_1c, basis_type="GAPW_1C")
318 :
319 78476 : IF (paw_atom) THEN
320 : !=========== PAW ===============
321 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
322 : maxso=maxso, npgf=npgf, maxl=maxl, &
323 29782 : nset=nset, zet=zet)
324 :
325 29782 : max_s_harm = harmonics%max_s_harm
326 29782 : llmax = harmonics%llmax
327 :
328 29782 : nsotot = maxso*nset
329 119128 : ALLOCATE (gsph(nr, nsotot))
330 89346 : ALLOCATE (gexp(nr))
331 148910 : ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0))
332 :
333 : NULLIFY (Vh1_h, Vh1_s)
334 119128 : ALLOCATE (Vh1_h(nr, max_iso_not0))
335 89346 : ALLOCATE (Vh1_s(nr, max_iso_not0))
336 :
337 119128 : ALLOCATE (aVh1b_hh(nsotot, nsotot))
338 89346 : ALLOCATE (aVh1b_ss(nsotot, nsotot))
339 89346 : ALLOCATE (aVh1b_00(nsotot, nsotot))
340 178692 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
341 :
342 29782 : NULLIFY (Qlm_gg, g0_h)
343 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
344 : l0_ikind=lmax0, &
345 29782 : Qlm_gg=Qlm_gg, g0_h=g0_h) ! Qlm_gg of density
346 :
347 29782 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
348 284 : NULLIFY (Qlm_gg_2nd, g0_h_2nd)
349 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole_2nd, ikind=ikind, &
350 : l0_ikind=lmax0_2nd, &
351 284 : Qlm_gg=Qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
352 : END IF
353 29782 : nchan_0 = nsoset(lmax0)
354 :
355 29782 : IF (nchan_0 > max_iso_not0) CPABORT("channels for rho0 > # max of spherical harmonics")
356 :
357 29782 : NULLIFY (rrad_z, my_CG)
358 29782 : my_CG => harmonics%my_CG
359 :
360 : ! set to zero temporary arrays
361 1746202 : sqrtwr = 0.0_dp
362 5244428 : g0_h_w = 0.0_dp
363 1746202 : gexp = 0.0_dp
364 56714306 : gsph = 0.0_dp
365 :
366 1746202 : sqrtwr(1:nr) = SQRT(grid_atom%wr(1:nr))
367 110516 : DO l_ang = 0, lmax0
368 4683416 : g0_h_w(1:nr, l_ang) = g0_h(1:nr, l_ang)*grid_atom%wr(1:nr)
369 : END DO
370 :
371 : m1 = 0
372 107306 : DO iset1 = 1, nset
373 77524 : n1 = nsoset(lmax(iset1))
374 283716 : DO ipgf1 = 1, npgf(iset1)
375 12037492 : gexp(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
376 823082 : DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1))
377 539366 : iso = is1 + (ipgf1 - 1)*n1 + m1
378 539366 : l_ang = indso(1, is1)
379 62092838 : gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
380 : END DO ! is1
381 : END DO ! ipgf1
382 107306 : m1 = m1 + maxso
383 : END DO ! iset1
384 :
385 : ! Distribute the atoms of this kind
386 29782 : num_pe = para_env%num_pe
387 29782 : mepos = para_env%mepos
388 29782 : bo = get_limit(nat, num_pe, mepos)
389 :
390 53867 : DO iat = bo(1), bo(2) !1,nat
391 24085 : iatom = atom_list(iat)
392 24085 : rho_atom => rho_atom_set(iatom)
393 :
394 24085 : NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
395 24085 : IF (core_charge) THEN
396 21185 : rrad_z => rhoz_set(ikind)%r_coef ! for density
397 : END IF
398 24085 : IF (my_core_2nd) THEN
399 21288 : IF (l_2nd_local_rho) THEN
400 103 : vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
401 : ELSE
402 21185 : vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
403 : END IF
404 : END IF
405 24085 : rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
406 24085 : vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
407 24085 : IF (l_2nd_local_rho) THEN
408 206 : rho_atom => rho_atom_set_2nd(iatom)
409 206 : vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
410 : END IF
411 24085 : IF (my_periodic .AND. back_ch .GT. 1.E-3_dp) THEN
412 896 : factor = -2.0_dp*pi/3.0_dp*SQRT(fourpi)*qs_charges%background
413 : ELSE
414 23189 : factor = 0._dp
415 : END IF
416 :
417 : CALL Vh_1c_atom_potential(rho_atom, vrrad_0, &
418 : grid_atom, my_core_2nd, vrrad_z, Vh1_h, Vh1_s, & ! core charge for potential (2nd)
419 24085 : nchan_0, nspins, max_iso_not0, factor)
420 :
421 24085 : IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom) ! rho_atom for density
422 :
423 : CALL Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
424 : grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, & ! core charge for density
425 24085 : nchan_0, nspins, max_iso_not0)
426 :
427 24085 : IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)
428 :
429 : CALL Vh_1c_atom_integrals(rho_atom, & ! results (int_local_h and int_local_s) written on rho_atom_2nd
430 : ! int_local_h and int_local_s are used in update_ks_atom
431 : ! on int_local_h mixed core / non-core
432 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
433 : max_s_harm, llmax, cg_list, cg_n_list, &
434 : nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
435 53867 : g0_h_w, my_CG, Qlm_gg) ! Qlm_gg for density from local_rho_set
436 :
437 : END DO ! iat
438 :
439 29782 : DEALLOCATE (aVh1b_hh)
440 29782 : DEALLOCATE (aVh1b_ss)
441 29782 : DEALLOCATE (aVh1b_00)
442 29782 : DEALLOCATE (Vh1_h, Vh1_s)
443 29782 : DEALLOCATE (cg_list, cg_n_list)
444 29782 : DEALLOCATE (gsph)
445 29782 : DEALLOCATE (gexp)
446 89346 : DEALLOCATE (sqrtwr, g0_h_w)
447 :
448 : ELSE
449 : !=========== NO PAW ===============
450 : ! This term is taken care of using the core density as in GPW
451 1584 : CYCLE
452 : END IF ! paw
453 : END DO ! ikind
454 :
455 15744 : CALL para_env%sum(energy_hartree_1c)
456 :
457 15744 : CALL timestop(handle)
458 :
459 31488 : END SUBROUTINE Vh_1c_gg_integrals
460 :
461 : ! **************************************************************************************************
462 :
463 : ! **************************************************************************************************
464 : !> \brief ...
465 : !> \param rho_atom ...
466 : !> \param vrrad_0 ...
467 : !> \param grid_atom ...
468 : !> \param core_charge ...
469 : !> \param vrrad_z ...
470 : !> \param Vh1_h ...
471 : !> \param Vh1_s ...
472 : !> \param nchan_0 ...
473 : !> \param nspins ...
474 : !> \param max_iso_not0 ...
475 : !> \param bfactor ...
476 : ! **************************************************************************************************
477 24085 : SUBROUTINE Vh_1c_atom_potential(rho_atom, vrrad_0, &
478 : grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, &
479 : nchan_0, nspins, max_iso_not0, bfactor)
480 :
481 : TYPE(rho_atom_type), POINTER :: rho_atom
482 : REAL(dp), DIMENSION(:, :), POINTER :: vrrad_0
483 : TYPE(grid_atom_type), POINTER :: grid_atom
484 : LOGICAL, INTENT(IN) :: core_charge
485 : REAL(dp), DIMENSION(:), POINTER :: vrrad_z
486 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
487 : INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
488 : REAL(dp), INTENT(IN) :: bfactor
489 :
490 : INTEGER :: ir, iso, ispin, nr
491 24085 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s
492 :
493 24085 : nr = grid_atom%nr
494 :
495 24085 : NULLIFY (vr_h, vr_s)
496 24085 : CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
497 :
498 17237584 : Vh1_h = 0.0_dp
499 17237584 : Vh1_s = 0.0_dp
500 :
501 2492453 : IF (core_charge) Vh1_h(:, 1) = vrrad_z(:)
502 :
503 204850 : DO iso = 1, nchan_0
504 20080235 : Vh1_s(:, iso) = vrrad_0(:, iso)
505 : END DO
506 :
507 52426 : DO ispin = 1, nspins
508 425831 : DO iso = 1, max_iso_not0
509 21565275 : Vh1_h(:, iso) = Vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
510 21593616 : Vh1_s(:, iso) = Vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
511 : END DO
512 : END DO
513 :
514 24085 : IF (bfactor /= 0._dp) THEN
515 53896 : DO ir = 1, nr
516 53000 : Vh1_h(ir, 1) = Vh1_h(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
517 53896 : Vh1_s(ir, 1) = Vh1_s(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
518 : END DO
519 : END IF
520 :
521 24085 : END SUBROUTINE Vh_1c_atom_potential
522 :
523 : ! **************************************************************************************************
524 :
525 : ! **************************************************************************************************
526 : !> \brief ...
527 : !> \param energy_hartree_1c ...
528 : !> \param ecoul_1c ...
529 : !> \param rho_atom ...
530 : !> \param rrad_0 ...
531 : !> \param grid_atom ...
532 : !> \param iatom ...
533 : !> \param core_charge ...
534 : !> \param rrad_z ...
535 : !> \param Vh1_h ...
536 : !> \param Vh1_s ...
537 : !> \param nchan_0 ...
538 : !> \param nspins ...
539 : !> \param max_iso_not0 ...
540 : ! **************************************************************************************************
541 24085 : SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
542 : grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
543 : nchan_0, nspins, max_iso_not0)
544 :
545 : REAL(dp), INTENT(INOUT) :: energy_hartree_1c
546 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
547 : TYPE(rho_atom_type), POINTER :: rho_atom
548 : REAL(dp), DIMENSION(:, :), POINTER :: rrad_0
549 : TYPE(grid_atom_type), POINTER :: grid_atom
550 : INTEGER, INTENT(IN) :: iatom
551 : LOGICAL, INTENT(IN) :: core_charge
552 : REAL(dp), DIMENSION(:), POINTER :: rrad_z
553 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
554 : INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
555 :
556 : INTEGER :: iso, ispin, nr
557 : REAL(dp) :: ecoul_1_0, ecoul_1_h, ecoul_1_s, &
558 : ecoul_1_z
559 24085 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
560 :
561 24085 : nr = grid_atom%nr
562 :
563 24085 : NULLIFY (r_h, r_s)
564 24085 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
565 :
566 : ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz
567 24085 : ecoul_1_z = 0.0_dp
568 24085 : IF (core_charge) THEN
569 1239575 : ecoul_1_z = 0.5_dp*SUM(Vh1_h(:, 1)*rrad_z(:)*grid_atom%wr(:))
570 : END IF
571 :
572 : ! Calculate the contributions to Ecoul coming from Vh1_s*rho0
573 24085 : ecoul_1_0 = 0.0_dp
574 204850 : DO iso = 1, nchan_0
575 10052160 : ecoul_1_0 = ecoul_1_0 + 0.5_dp*SUM(Vh1_s(:, iso)*rrad_0(:, iso)*grid_atom%wr(:))
576 : END DO
577 :
578 : ! Calculate the contributions to Ecoul coming from Vh1_h*rho1_h and Vh1_s*rho1_s
579 24085 : ecoul_1_s = 0.0_dp
580 24085 : ecoul_1_h = 0.0_dp
581 52426 : DO ispin = 1, nspins
582 425831 : DO iso = 1, max_iso_not0
583 21565275 : ecoul_1_s = ecoul_1_s + 0.5_dp*SUM(Vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
584 21593616 : ecoul_1_h = ecoul_1_h + 0.5_dp*SUM(Vh1_h(:, iso)*r_h(ispin)%r_coef(:, iso)*grid_atom%wr(:))
585 : END DO
586 : END DO
587 :
588 24085 : CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
589 24085 : CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
590 :
591 24085 : energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
592 24085 : energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
593 :
594 24085 : END SUBROUTINE Vh_1c_atom_energy
595 :
596 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
597 :
598 : ! **************************************************************************************************
599 : !> \brief ...
600 : !> \param rho_atom ...
601 : !> \param aVh1b_hh ...
602 : !> \param aVh1b_ss ...
603 : !> \param aVh1b_00 ...
604 : !> \param Vh1_h ...
605 : !> \param Vh1_s ...
606 : !> \param max_iso_not0 ...
607 : !> \param max_s_harm ...
608 : !> \param llmax ...
609 : !> \param cg_list ...
610 : !> \param cg_n_list ...
611 : !> \param nset ...
612 : !> \param npgf ...
613 : !> \param lmin ...
614 : !> \param lmax ...
615 : !> \param nsotot ...
616 : !> \param maxso ...
617 : !> \param nspins ...
618 : !> \param nchan_0 ...
619 : !> \param gsph ...
620 : !> \param g0_h_w ...
621 : !> \param my_CG ...
622 : !> \param Qlm_gg ...
623 : ! **************************************************************************************************
624 24085 : SUBROUTINE Vh_1c_atom_integrals(rho_atom, &
625 24085 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
626 24085 : max_s_harm, llmax, cg_list, cg_n_list, &
627 : nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
628 24085 : g0_h_w, my_CG, Qlm_gg)
629 :
630 : TYPE(rho_atom_type), POINTER :: rho_atom
631 : REAL(dp), DIMENSION(:, :) :: aVh1b_hh, aVh1b_ss, aVh1b_00
632 : REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s
633 : INTEGER, INTENT(IN) :: max_iso_not0, max_s_harm, llmax
634 : INTEGER, DIMENSION(:, :, :) :: cg_list
635 : INTEGER, DIMENSION(:) :: cg_n_list
636 : INTEGER, INTENT(IN) :: nset
637 : INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax
638 : INTEGER, INTENT(IN) :: nsotot, maxso, nspins, nchan_0
639 : REAL(dp), DIMENSION(:, :), POINTER :: gsph
640 : REAL(dp), DIMENSION(:, 0:) :: g0_h_w
641 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg
642 :
643 : INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, &
644 : iset2, iso, iso1, iso2, ispin, l_ang, &
645 : m1, m2, max_iso_not0_local, n1, n2, nr
646 : REAL(dp) :: gVg_0, gVg_h, gVg_s
647 24085 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
648 :
649 24085 : NULLIFY (int_local_h, int_local_s)
650 : CALL get_rho_atom(rho_atom=rho_atom, &
651 : ga_Vlocal_gb_h=int_local_h, &
652 24085 : ga_Vlocal_gb_s=int_local_s)
653 :
654 : ! Calculate the integrals of the potential with 2 primitives
655 72529203 : aVh1b_hh = 0.0_dp
656 72529203 : aVh1b_ss = 0.0_dp
657 72529203 : aVh1b_00 = 0.0_dp
658 :
659 24085 : nr = SIZE(gsph, 1)
660 :
661 24085 : m1 = 0
662 84052 : DO iset1 = 1, nset
663 : m2 = 0
664 259120 : DO iset2 = 1, nset
665 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
666 199153 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
667 :
668 199153 : n1 = nsoset(lmax(iset1))
669 683058 : DO ipgf1 = 1, npgf(iset1)
670 483905 : n2 = nsoset(lmax(iset2))
671 2073385 : DO ipgf2 = 1, npgf(iset2)
672 : ! with contributions to V1_s*rho0
673 13721318 : DO iso = 1, nchan_0
674 12330991 : l_ang = indso(1, iso)
675 685474151 : gVg_0 = SUM(Vh1_s(:, iso)*g0_h_w(:, l_ang))
676 25123130 : DO icg = 1, cg_n_list(iso)
677 11401812 : is1 = cg_list(1, icg, iso)
678 11401812 : is2 = cg_list(2, icg, iso)
679 :
680 11401812 : iso1 = is1 + n1*(ipgf1 - 1) + m1
681 11401812 : iso2 = is2 + n2*(ipgf2 - 1) + m2
682 11401812 : gVg_h = 0.0_dp
683 11401812 : gVg_s = 0.0_dp
684 :
685 631405372 : DO ir = 1, nr
686 620003560 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
687 631405372 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
688 : END DO ! ir
689 :
690 11401812 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
691 11401812 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
692 23732803 : aVh1b_00(iso1, iso2) = aVh1b_00(iso1, iso2) + gVg_0*Qlm_gg(iso1, iso2, iso)
693 :
694 : END DO !icg
695 : END DO ! iso
696 : ! without contributions to V1_s*rho0
697 20361776 : DO iso = nchan_0 + 1, max_iso_not0
698 25338324 : DO icg = 1, cg_n_list(iso)
699 5460453 : is1 = cg_list(1, icg, iso)
700 5460453 : is2 = cg_list(2, icg, iso)
701 :
702 5460453 : iso1 = is1 + n1*(ipgf1 - 1) + m1
703 5460453 : iso2 = is2 + n2*(ipgf2 - 1) + m2
704 5460453 : gVg_h = 0.0_dp
705 5460453 : gVg_s = 0.0_dp
706 :
707 287392333 : DO ir = 1, nr
708 281931880 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
709 287392333 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
710 : END DO ! ir
711 :
712 5460453 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
713 23947997 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
714 :
715 : END DO !icg
716 : END DO ! iso
717 : END DO ! ipgf2
718 : END DO ! ipgf1
719 259120 : m2 = m2 + maxso
720 : END DO ! iset2
721 84052 : m1 = m1 + maxso
722 : END DO !iset1
723 52426 : DO ispin = 1, nspins
724 28341 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
725 28341 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
726 28341 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_h(ispin)%r_coef, 1)
727 52426 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_s(ispin)%r_coef, 1)
728 : END DO ! ispin
729 :
730 24085 : END SUBROUTINE Vh_1c_atom_integrals
731 :
732 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
733 :
734 : END MODULE hartree_local_methods
735 :
|
