(git:374b731)
Loading...
Searching...
No Matches
hartree_local_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
13 USE cell_types, ONLY: cell_type
19 USE kinds, ONLY: dp
20 USE mathconstants, ONLY: fourpi,&
21 pi
23 USE orbital_pointers, ONLY: indso,&
24 nsoset
25 USE pw_env_types, ONLY: pw_env_get,&
35 USE qs_kind_types, ONLY: get_qs_kind,&
41 USE qs_rho0_types, ONLY: get_rho0_mpole,&
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
59
60CONTAINS
61
62! **************************************************************************************************
63!> \brief ...
64!> \param hartree_local ...
65!> \param natom ...
66! **************************************************************************************************
67 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 TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
76
77 CALL timeset(routinen, handle)
78
79 NULLIFY (ecoul_1c)
80 ! Allocate and Initialize 1-center Potentials and Integrals
81 CALL allocate_ecoul_1center(ecoul_1c, natom)
82 hartree_local%ecoul_1c => ecoul_1c
83
84 CALL timestop(handle)
85
86 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rho_1, rho_2
116 REAL(dp), DIMENSION(:), POINTER :: wr
117 REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l
118
119 CALL timeset(routinen, handle)
120
121 nr = grid_atom%nr
122 max_s_harm = SIZE(vrad_h, 2)
123 nspins = SIZE(rrad_h, 1)
124 nchannels = SIZE(rrad_0, 2)
125
126 r2l => grid_atom%rad2l
127 oor2l => grid_atom%oorad2l
128 wr => grid_atom%wr
129
130 ALLOCATE (rho_1(nr, max_s_harm), rho_2(nr, max_s_harm))
131 rho_1 = 0.0_dp
132 rho_2 = 0.0_dp
133
134 ! Case lm = 0
135 rho_1(:, 1) = rrad_z(:)
136 rho_2(:, 1) = rrad_0(:, 1)
137
138 DO iso = 2, nchannels
139 rho_2(:, iso) = rrad_0(:, iso)
140 END DO
141
142 DO iso = 1, max_s_harm
143 DO ispin = 1, nspins
144 rho_1(:, iso) = rho_1(:, iso) + rrad_h(ispin)%r_coef(:, iso)
145 rho_2(:, iso) = rho_2(:, iso) + rrad_s(ispin)%r_coef(:, iso)
146 END DO
147
148 l_ang = indso(1, iso)
149 prefactor = fourpi/(2._dp*l_ang + 1._dp)
150
151 rho_1(:, iso) = rho_1(:, iso)*wr(:)
152 rho_2(:, iso) = rho_2(:, iso)*wr(:)
153
154 i1_up = 0.0_dp
155 i1_down = 0.0_dp
156 i2_up = 0.0_dp
157 i2_down = 0.0_dp
158
159 i1_up = r2l(nr, l_ang)*rho_1(nr, iso)
160 i2_up = r2l(nr, l_ang)*rho_2(nr, iso)
161
162 DO ir = nr - 1, 1, -1
163 i1_down = i1_down + oor2l(ir, l_ang + 1)*rho_1(ir, iso)
164 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 (oor2l(nr, l_ang + 1)*i1_up + r2l(nr, l_ang)*i1_down)
169 vrad_s(nr, iso) = vrad_s(nr, iso) + prefactor* &
170 (oor2l(nr, l_ang + 1)*i2_up + r2l(nr, l_ang)*i2_down)
171
172 DO ir = nr - 1, 1, -1
173 i1_up = i1_up + r2l(ir, l_ang)*rho_1(ir, iso)
174 i1_down = i1_down - oor2l(ir, l_ang + 1)*rho_1(ir, iso)
175 i2_up = i2_up + r2l(ir, l_ang)*rho_2(ir, iso)
176 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 (oor2l(ir, l_ang + 1)*i1_up + r2l(ir, l_ang)*i1_down)
180 vrad_s(ir, iso) = vrad_s(ir, iso) + prefactor* &
181 (oor2l(ir, l_ang + 1)*i2_up + r2l(ir, l_ang)*i2_down)
182
183 END DO
184
185 END DO
186
187 DEALLOCATE (rho_1, rho_2)
188
189 CALL timestop(handle)
190
191 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
228 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
229 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 REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr
234 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: avh1b_00, avh1b_hh, avh1b_ss, g0_h_w
235 REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
236 REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, rrad_0, vh1_h, &
237 vh1_s, vrrad_0, zet
238 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg, qlm_gg, qlm_gg_2nd
239 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
249 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 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
252 TYPE(rho_atom_type), POINTER :: rho_atom
253 TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
254
255 CALL timeset(routinen, handle)
256
257 NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
258 NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
259 NULLIFY (rho0_mpole, rhoz_set)
260 NULLIFY (atom_list, grid_atom, harmonics)
261 NULLIFY (basis_1c, lmin, lmax, npgf, zet)
262 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 pw_env=pw_env, qs_charges=qs_charges)
269
270 CALL pw_env_get(pw_env, poisson_env=poisson_env)
271 my_periodic = (poisson_env%method == pw_poisson_periodic)
272
273 back_ch = qs_charges%background*cell%deth
274
275 ! rhoz_set is not accessed in TDDFT
276 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 l_2nd_local_rho = .false.
280 IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
281 l_2nd_local_rho = .true.
282 NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
283 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 nkind = SIZE(atomic_kind_set, 1)
287 nspins = dft_control%nspins
288
289 core_charge = .NOT. tddft ! for forces mixed version
290 my_core_2nd = .true.
291 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 CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
303 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 energy_hartree_1c = 0.0_dp
307
308 ! Here starts the loop over all the atoms
309 DO ikind = 1, nkind
310
311 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 max_iso_not0=max_iso_not0, paw_atom=paw_atom)
316 CALL get_qs_kind(qs_kind_set(ikind), &
317 basis_set=basis_1c, basis_type="GAPW_1C")
318
319 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 nset=nset, zet=zet)
324
325 max_s_harm = harmonics%max_s_harm
326 llmax = harmonics%llmax
327
328 nsotot = maxso*nset
329 ALLOCATE (gsph(nr, nsotot))
330 ALLOCATE (gexp(nr))
331 ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0))
332
333 NULLIFY (vh1_h, vh1_s)
334 ALLOCATE (vh1_h(nr, max_iso_not0))
335 ALLOCATE (vh1_s(nr, max_iso_not0))
336
337 ALLOCATE (avh1b_hh(nsotot, nsotot))
338 ALLOCATE (avh1b_ss(nsotot, nsotot))
339 ALLOCATE (avh1b_00(nsotot, nsotot))
340 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
341
342 NULLIFY (qlm_gg, g0_h)
343 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
344 l0_ikind=lmax0, &
345 qlm_gg=qlm_gg, g0_h=g0_h) ! Qlm_gg of density
346
347 IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
348 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 qlm_gg=qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
352 END IF
353 nchan_0 = nsoset(lmax0)
354
355 IF (nchan_0 > max_iso_not0) cpabort("channels for rho0 > # max of spherical harmonics")
356
357 NULLIFY (rrad_z, my_cg)
358 my_cg => harmonics%my_CG
359
360 ! set to zero temporary arrays
361 sqrtwr = 0.0_dp
362 g0_h_w = 0.0_dp
363 gexp = 0.0_dp
364 gsph = 0.0_dp
365
366 sqrtwr(1:nr) = sqrt(grid_atom%wr(1:nr))
367 DO l_ang = 0, lmax0
368 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 DO iset1 = 1, nset
373 n1 = nsoset(lmax(iset1))
374 DO ipgf1 = 1, npgf(iset1)
375 gexp(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
376 DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1))
377 iso = is1 + (ipgf1 - 1)*n1 + m1
378 l_ang = indso(1, is1)
379 gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
380 END DO ! is1
381 END DO ! ipgf1
382 m1 = m1 + maxso
383 END DO ! iset1
384
385 ! Distribute the atoms of this kind
386 num_pe = para_env%num_pe
387 mepos = para_env%mepos
388 bo = get_limit(nat, num_pe, mepos)
389
390 DO iat = bo(1), bo(2) !1,nat
391 iatom = atom_list(iat)
392 rho_atom => rho_atom_set(iatom)
393
394 NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
395 IF (core_charge) THEN
396 rrad_z => rhoz_set(ikind)%r_coef ! for density
397 END IF
398 IF (my_core_2nd) THEN
399 IF (l_2nd_local_rho) THEN
400 vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
401 ELSE
402 vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
403 END IF
404 END IF
405 rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
406 vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
407 IF (l_2nd_local_rho) THEN
408 rho_atom => rho_atom_set_2nd(iatom)
409 vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
410 END IF
411 IF (my_periodic .AND. back_ch .GT. 1.e-3_dp) THEN
412 factor = -2.0_dp*pi/3.0_dp*sqrt(fourpi)*qs_charges%background
413 ELSE
414 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 nchan_0, nspins, max_iso_not0, factor)
420
421 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 nchan_0, nspins, max_iso_not0)
426
427 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 g0_h_w, my_cg, qlm_gg) ! Qlm_gg for density from local_rho_set
436
437 END DO ! iat
438
439 DEALLOCATE (avh1b_hh)
440 DEALLOCATE (avh1b_ss)
441 DEALLOCATE (avh1b_00)
442 DEALLOCATE (vh1_h, vh1_s)
443 DEALLOCATE (cg_list, cg_n_list)
444 DEALLOCATE (gsph)
445 DEALLOCATE (gexp)
446 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 cycle
452 END IF ! paw
453 END DO ! ikind
454
455 CALL para_env%sum(energy_hartree_1c)
456
457 CALL timestop(handle)
458
459 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 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 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s
492
493 nr = grid_atom%nr
494
495 NULLIFY (vr_h, vr_s)
496 CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
497
498 vh1_h = 0.0_dp
499 vh1_s = 0.0_dp
500
501 IF (core_charge) vh1_h(:, 1) = vrrad_z(:)
502
503 DO iso = 1, nchan_0
504 vh1_s(:, iso) = vrrad_0(:, iso)
505 END DO
506
507 DO ispin = 1, nspins
508 DO iso = 1, max_iso_not0
509 vh1_h(:, iso) = vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
510 vh1_s(:, iso) = vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
511 END DO
512 END DO
513
514 IF (bfactor /= 0._dp) THEN
515 DO ir = 1, nr
516 vh1_h(ir, 1) = vh1_h(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
517 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 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 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 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
560
561 nr = grid_atom%nr
562
563 NULLIFY (r_h, r_s)
564 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 ecoul_1_z = 0.0_dp
568 IF (core_charge) THEN
569 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 ecoul_1_0 = 0.0_dp
574 DO iso = 1, nchan_0
575 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 ecoul_1_s = 0.0_dp
580 ecoul_1_h = 0.0_dp
581 DO ispin = 1, nspins
582 DO iso = 1, max_iso_not0
583 ecoul_1_s = ecoul_1_s + 0.5_dp*sum(vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
584 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 CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
589 CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
590
591 energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
592 energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
593
594 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 SUBROUTINE vh_1c_atom_integrals(rho_atom, &
625 aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
626 max_s_harm, llmax, cg_list, cg_n_list, &
627 nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
628 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 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
648
649 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 ga_vlocal_gb_s=int_local_s)
653
654 ! Calculate the integrals of the potential with 2 primitives
655 avh1b_hh = 0.0_dp
656 avh1b_ss = 0.0_dp
657 avh1b_00 = 0.0_dp
658
659 nr = SIZE(gsph, 1)
660
661 m1 = 0
662 DO iset1 = 1, nset
663 m2 = 0
664 DO iset2 = 1, nset
665 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
666 max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
667
668 n1 = nsoset(lmax(iset1))
669 DO ipgf1 = 1, npgf(iset1)
670 n2 = nsoset(lmax(iset2))
671 DO ipgf2 = 1, npgf(iset2)
672 ! with contributions to V1_s*rho0
673 DO iso = 1, nchan_0
674 l_ang = indso(1, iso)
675 gvg_0 = sum(vh1_s(:, iso)*g0_h_w(:, l_ang))
676 DO icg = 1, cg_n_list(iso)
677 is1 = cg_list(1, icg, iso)
678 is2 = cg_list(2, icg, iso)
679
680 iso1 = is1 + n1*(ipgf1 - 1) + m1
681 iso2 = is2 + n2*(ipgf2 - 1) + m2
682 gvg_h = 0.0_dp
683 gvg_s = 0.0_dp
684
685 DO ir = 1, nr
686 gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
687 gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
688 END DO ! ir
689
690 avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
691 avh1b_ss(iso1, iso2) = avh1b_ss(iso1, iso2) + gvg_s*my_cg(is1, is2, iso)
692 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 DO iso = nchan_0 + 1, max_iso_not0
698 DO icg = 1, cg_n_list(iso)
699 is1 = cg_list(1, icg, iso)
700 is2 = cg_list(2, icg, iso)
701
702 iso1 = is1 + n1*(ipgf1 - 1) + m1
703 iso2 = is2 + n2*(ipgf2 - 1) + m2
704 gvg_h = 0.0_dp
705 gvg_s = 0.0_dp
706
707 DO ir = 1, nr
708 gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
709 gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
710 END DO ! ir
711
712 avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
713 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 m2 = m2 + maxso
720 END DO ! iset2
721 m1 = m1 + maxso
722 END DO !iset1
723 DO ispin = 1, nspins
724 CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
725 CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
726 CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_h(ispin)%r_coef, 1)
727 CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_s(ispin)%r_coef, 1)
728 END DO ! ispin
729
730 END SUBROUTINE vh_1c_atom_integrals
731
732!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
733
734END MODULE hartree_local_methods
735
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Handles all functions related to the CELL.
Definition cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public init_coulomb_local(hartree_local, natom)
...
subroutine, public vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, core_2nd)
Calculates one center GAPW Hartree energies and matrix elements Hartree potentials are input Takes po...
subroutine, public calculate_vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom)
Calculates Hartree potential for hard and soft densities (including nuclear charge and compensation c...
subroutine, public allocate_ecoul_1center(ecoul_1c, natom)
...
subroutine, public set_ecoul_1c(ecoul_1c, iatom, ecoul_1_h, ecoul_1_s, ecoul_1_z, ecoul_1_0)
...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public fourpi
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nsoset
integer, dimension(:, :), allocatable, public indso
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_periodic
container for information about total charges on the grids
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
subroutine, public get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set)
...
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, qlm_gg, qlm_car, qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs)
...
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
All kind of helpful little routines.
Definition util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
Container for information about total charges on the grids.
Provides all information about a quickstep kind.