Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 1612 : 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 1612 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
76 :
77 1612 : CALL timeset(routineN, handle)
78 :
79 1612 : NULLIFY (ecoul_1c)
80 : ! Allocate and Initialize 1-center Potentials and Integrals
81 1612 : CALL allocate_ecoul_1center(ecoul_1c, natom)
82 1612 : hartree_local%ecoul_1c => ecoul_1c
83 :
84 1612 : CALL timestop(handle)
85 :
86 1612 : 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 14954 : 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 14954 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
228 14954 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
229 14954 : 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 14954 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr
234 14954 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_hh, aVh1b_ss, g0_h_w
235 14954 : REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
236 14954 : REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, rrad_0, Vh1_h, &
237 14954 : Vh1_s, vrrad_0, zet
238 14954 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg, Qlm_gg_2nd
239 14954 : 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 14954 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
249 14954 : 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 14954 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
252 : TYPE(rho_atom_type), POINTER :: rho_atom
253 14954 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
254 :
255 14954 : CALL timeset(routineN, handle)
256 :
257 14954 : NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
258 14954 : NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
259 14954 : NULLIFY (rho0_mpole, rhoz_set)
260 14954 : NULLIFY (atom_list, grid_atom, harmonics)
261 14954 : NULLIFY (basis_1c, lmin, lmax, npgf, zet)
262 14954 : 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 14954 : pw_env=pw_env, qs_charges=qs_charges)
269 :
270 14954 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
271 14954 : my_periodic = (poisson_env%method == pw_poisson_periodic)
272 :
273 14954 : back_ch = qs_charges%background*cell%deth
274 :
275 : ! rhoz_set is not accessed in TDDFT
276 14954 : 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 14954 : l_2nd_local_rho = .FALSE.
280 14954 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
281 124 : l_2nd_local_rho = .TRUE.
282 124 : NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
283 124 : 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 14954 : nkind = SIZE(atomic_kind_set, 1)
287 14954 : nspins = dft_control%nspins
288 :
289 14954 : core_charge = .NOT. tddft ! for forces mixed version
290 14954 : my_core_2nd = .TRUE.
291 14954 : 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 14954 : CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
303 14954 : 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 14954 : energy_hartree_1c = 0.0_dp
307 :
308 : ! Here starts the loop over all the atoms
309 44950 : DO ikind = 1, nkind
310 :
311 29996 : 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 29996 : max_iso_not0=max_iso_not0, paw_atom=paw_atom)
316 : CALL get_qs_kind(qs_kind_set(ikind), &
317 29996 : basis_set=basis_1c, basis_type="GAPW_1C")
318 :
319 74946 : 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 28498 : nset=nset, zet=zet)
324 :
325 28498 : max_s_harm = harmonics%max_s_harm
326 28498 : llmax = harmonics%llmax
327 :
328 28498 : nsotot = maxso*nset
329 113992 : ALLOCATE (gsph(nr, nsotot))
330 85494 : ALLOCATE (gexp(nr))
331 142490 : ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0))
332 :
333 : NULLIFY (Vh1_h, Vh1_s)
334 113992 : ALLOCATE (Vh1_h(nr, max_iso_not0))
335 85494 : ALLOCATE (Vh1_s(nr, max_iso_not0))
336 :
337 113992 : ALLOCATE (aVh1b_hh(nsotot, nsotot))
338 85494 : ALLOCATE (aVh1b_ss(nsotot, nsotot))
339 85494 : ALLOCATE (aVh1b_00(nsotot, nsotot))
340 170988 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
341 :
342 28498 : NULLIFY (Qlm_gg, g0_h)
343 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
344 : l0_ikind=lmax0, &
345 28498 : Qlm_gg=Qlm_gg, g0_h=g0_h) ! Qlm_gg of density
346 :
347 28498 : IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
348 276 : 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 276 : Qlm_gg=Qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
352 : END IF
353 28498 : nchan_0 = nsoset(lmax0)
354 :
355 28498 : IF (nchan_0 > max_iso_not0) CPABORT("channels for rho0 > # max of spherical harmonics")
356 :
357 28498 : NULLIFY (rrad_z, my_CG)
358 28498 : my_CG => harmonics%my_CG
359 :
360 : ! set to zero temporary arrays
361 1679718 : sqrtwr = 0.0_dp
362 5045528 : g0_h_w = 0.0_dp
363 1679718 : gexp = 0.0_dp
364 52819922 : gsph = 0.0_dp
365 :
366 1679718 : sqrtwr(1:nr) = SQRT(grid_atom%wr(1:nr))
367 105604 : DO l_ang = 0, lmax0
368 4494104 : 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 101340 : DO iset1 = 1, nset
373 72842 : n1 = nsoset(lmax(iset1))
374 267778 : DO ipgf1 = 1, npgf(iset1)
375 11456436 : gexp(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
376 774568 : DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1))
377 506790 : iso = is1 + (ipgf1 - 1)*n1 + m1
378 506790 : l_ang = indso(1, is1)
379 58773406 : gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
380 : END DO ! is1
381 : END DO ! ipgf1
382 101340 : m1 = m1 + maxso
383 : END DO ! iset1
384 :
385 : ! Distribute the atoms of this kind
386 28498 : num_pe = para_env%num_pe
387 28498 : mepos = para_env%mepos
388 28498 : bo = get_limit(nat, num_pe, mepos)
389 :
390 51253 : DO iat = bo(1), bo(2) !1,nat
391 22755 : iatom = atom_list(iat)
392 22755 : rho_atom => rho_atom_set(iatom)
393 :
394 22755 : NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
395 22755 : IF (core_charge) THEN
396 20043 : rrad_z => rhoz_set(ikind)%r_coef ! for density
397 : END IF
398 22755 : IF (my_core_2nd) THEN
399 20143 : IF (l_2nd_local_rho) THEN
400 100 : vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
401 : ELSE
402 20043 : vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
403 : END IF
404 : END IF
405 22755 : rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
406 22755 : vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
407 22755 : IF (l_2nd_local_rho) THEN
408 200 : rho_atom => rho_atom_set_2nd(iatom)
409 200 : vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
410 : END IF
411 22755 : IF (my_periodic .AND. back_ch .GT. 1.E-3_dp) THEN
412 808 : factor = -2.0_dp*pi/3.0_dp*SQRT(fourpi)*qs_charges%background
413 : ELSE
414 21947 : 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 22755 : nchan_0, nspins, max_iso_not0, factor)
420 :
421 22755 : 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 22755 : nchan_0, nspins, max_iso_not0)
426 :
427 22755 : 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 51253 : g0_h_w, my_CG, Qlm_gg) ! Qlm_gg for density from local_rho_set
436 :
437 : END DO ! iat
438 :
439 28498 : DEALLOCATE (aVh1b_hh)
440 28498 : DEALLOCATE (aVh1b_ss)
441 28498 : DEALLOCATE (aVh1b_00)
442 28498 : DEALLOCATE (Vh1_h, Vh1_s)
443 28498 : DEALLOCATE (cg_list, cg_n_list)
444 28498 : DEALLOCATE (gsph)
445 28498 : DEALLOCATE (gexp)
446 85494 : 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 1498 : CYCLE
452 : END IF ! paw
453 : END DO ! ikind
454 :
455 14954 : CALL para_env%sum(energy_hartree_1c)
456 :
457 14954 : CALL timestop(handle)
458 :
459 29908 : 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 22755 : 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 22755 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s
492 :
493 22755 : nr = grid_atom%nr
494 :
495 22755 : NULLIFY (vr_h, vr_s)
496 22755 : CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
497 :
498 16323840 : Vh1_h = 0.0_dp
499 16323840 : Vh1_s = 0.0_dp
500 :
501 2373478 : IF (core_charge) Vh1_h(:, 1) = vrrad_z(:)
502 :
503 192766 : DO iso = 1, nchan_0
504 18963997 : Vh1_s(:, iso) = vrrad_0(:, iso)
505 : END DO
506 :
507 49536 : DO ispin = 1, nspins
508 402685 : DO iso = 1, max_iso_not0
509 20514219 : Vh1_h(:, iso) = Vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
510 20541000 : Vh1_s(:, iso) = Vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
511 : END DO
512 : END DO
513 :
514 22755 : IF (bfactor /= 0._dp) THEN
515 48408 : DO ir = 1, nr
516 47600 : Vh1_h(ir, 1) = Vh1_h(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
517 48408 : 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 22755 : 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 22755 : 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 22755 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
560 :
561 22755 : nr = grid_atom%nr
562 :
563 22755 : NULLIFY (r_h, r_s)
564 22755 : 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 22755 : ecoul_1_z = 0.0_dp
568 22755 : IF (core_charge) THEN
569 1180333 : 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 22755 : ecoul_1_0 = 0.0_dp
574 192766 : DO iso = 1, nchan_0
575 9493376 : 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 22755 : ecoul_1_s = 0.0_dp
580 22755 : ecoul_1_h = 0.0_dp
581 49536 : DO ispin = 1, nspins
582 402685 : DO iso = 1, max_iso_not0
583 20514219 : ecoul_1_s = ecoul_1_s + 0.5_dp*SUM(Vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
584 20541000 : 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 22755 : CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
589 22755 : CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
590 :
591 22755 : energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
592 22755 : energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
593 :
594 22755 : 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 22755 : SUBROUTINE Vh_1c_atom_integrals(rho_atom, &
625 22755 : aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
626 22755 : max_s_harm, llmax, cg_list, cg_n_list, &
627 : nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
628 22755 : 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 22755 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
648 :
649 22755 : 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 22755 : ga_Vlocal_gb_s=int_local_s)
653 :
654 : ! Calculate the integrals of the potential with 2 primitives
655 65060861 : aVh1b_hh = 0.0_dp
656 65060861 : aVh1b_ss = 0.0_dp
657 65060861 : aVh1b_00 = 0.0_dp
658 :
659 22755 : nr = SIZE(gsph, 1)
660 :
661 22755 : m1 = 0
662 78597 : DO iset1 = 1, nset
663 : m2 = 0
664 236134 : DO iset2 = 1, nset
665 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
666 180292 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
667 :
668 180292 : n1 = nsoset(lmax(iset1))
669 624200 : DO ipgf1 = 1, npgf(iset1)
670 443908 : n2 = nsoset(lmax(iset2))
671 1918049 : DO ipgf2 = 1, npgf(iset2)
672 : ! with contributions to V1_s*rho0
673 12793450 : DO iso = 1, nchan_0
674 11499601 : l_ang = indso(1, iso)
675 642632261 : gVg_0 = SUM(Vh1_s(:, iso)*g0_h_w(:, l_ang))
676 23254925 : DO icg = 1, cg_n_list(iso)
677 10461475 : is1 = cg_list(1, icg, iso)
678 10461475 : is2 = cg_list(2, icg, iso)
679 :
680 10461475 : iso1 = is1 + n1*(ipgf1 - 1) + m1
681 10461475 : iso2 = is2 + n2*(ipgf2 - 1) + m2
682 10461475 : gVg_h = 0.0_dp
683 10461475 : gVg_s = 0.0_dp
684 :
685 583362185 : DO ir = 1, nr
686 572900710 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
687 583362185 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
688 : END DO ! ir
689 :
690 10461475 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
691 10461475 : aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso)
692 21961076 : 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 19052517 : DO iso = nchan_0 + 1, max_iso_not0
698 23584455 : DO icg = 1, cg_n_list(iso)
699 4975846 : is1 = cg_list(1, icg, iso)
700 4975846 : is2 = cg_list(2, icg, iso)
701 :
702 4975846 : iso1 = is1 + n1*(ipgf1 - 1) + m1
703 4975846 : iso2 = is2 + n2*(ipgf2 - 1) + m2
704 4975846 : gVg_h = 0.0_dp
705 4975846 : gVg_s = 0.0_dp
706 :
707 262677376 : DO ir = 1, nr
708 257701530 : gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso)
709 262677376 : gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso)
710 : END DO ! ir
711 :
712 4975846 : aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso)
713 22290606 : 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 236134 : m2 = m2 + maxso
720 : END DO ! iset2
721 78597 : m1 = m1 + maxso
722 : END DO !iset1
723 49536 : DO ispin = 1, nspins
724 26781 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
725 26781 : CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
726 26781 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_h(ispin)%r_coef, 1)
727 49536 : CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_s(ispin)%r_coef, 1)
728 : END DO ! ispin
729 :
730 22755 : END SUBROUTINE Vh_1c_atom_integrals
731 :
732 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
733 :
734 : END MODULE hartree_local_methods
735 :
|
