33 USE dbcsr_api,
ONLY: &
34 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_dot, dbcsr_get_info, &
35 dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_scale, dbcsr_set, &
36 dbcsr_type, dbcsr_type_symmetric
82 integrate_v_gaussian_rspace,&
104#include "./base/base_uses.f90"
109 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ed_analysis'
126 INTEGER,
INTENT(IN) :: unit_nr
128 CHARACTER(len=*),
PARAMETER :: routinen =
'edmf_analysis'
130 INTEGER :: handle, iatom, ikind, iorb, ispin, jorb, &
131 nao, natom, nimages, nkind, no, norb, &
133 INTEGER,
DIMENSION(2) :: nocc
134 INTEGER,
DIMENSION(:),
POINTER :: refbas_blk_sizes
135 LOGICAL :: detailed_ener, do_hfx, ewald_correction, explicit, gapw, gapw_xc, &
136 ref_orb_canonical, skip_localize, uniform_occupation, uocc
137 REAL(kind=
dp) :: ateps, checksum, e1, e2, e_pot, ealpha, &
138 ecc, egcp, ehfx, ekts, evdw, focc, &
140 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: amval, ate1h, ate1xc, atecc, ateks, &
141 atener, atewald, odiag
142 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: atdet, atmul, mcharge, mweight
143 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: bcenter
144 REAL(kind=
dp),
DIMENSION(:),
POINTER :: occupation_numbers
148 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: c_iao_coef, ciao, fij_mat, orb_weight, &
151 TYPE(dbcsr_distribution_type) :: dbcsr_dist
152 TYPE(dbcsr_p_type) :: dve_mat
153 TYPE(dbcsr_p_type),
ALLOCATABLE,
DIMENSION(:) :: exc_mat, ks_mat, vhxc_mat
154 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_h, matrix_hfx, matrix_ks, &
155 matrix_p, matrix_s, matrix_t
156 TYPE(dbcsr_type) :: dkmat, dmat
157 TYPE(dbcsr_type),
POINTER :: smat
162 TYPE(
mo_set_type),
DIMENSION(:),
POINTER :: mos, mos_loc
169 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
175 IF (.NOT. explicit)
RETURN
177 CALL timeset(routinen, handle)
179 IF (unit_nr > 0)
THEN
180 WRITE (unit=unit_nr, fmt=
"(/,4(T2,A))") &
181 "!-----------------------------------------------------------------------------!", &
182 "! ENERGY DECOMPOSITION ANALYSIS !", &
183 "! Janus J Eriksen, JCP 153 214109 (2020) !", &
184 "!-----------------------------------------------------------------------------!"
193 ewald_correction = (ealpha /= 0.0_dp)
195 CALL get_qs_env(qs_env, dft_control=dft_control)
196 nimages = dft_control%nimages
197 IF (nimages > 1)
THEN
198 IF (unit_nr > 0)
THEN
199 WRITE (unit=unit_nr, fmt=
"(T2,A)") &
200 "K-Points: MAO's determined and analyzed using Gamma-Point only."
203 gapw = dft_control%qs_control%gapw
204 gapw_xc = dft_control%qs_control%gapw_xc
205 IF (ewald_correction)
THEN
206 IF (gapw .OR. gapw_xc)
THEN
207 CALL cp_warn(__location__,
"Ewald Correction for GAPW and GAPW_XC not available.")
208 cpabort(
"Ewald Correction")
216 nspin = dft_control%nspins
217 ALLOCATE (mos_loc(nspin))
222 CALL get_mo_set(mos(ispin), uniform_occupation=uniform_occupation)
223 IF (.NOT. uniform_occupation) uocc = .false.
225 IF (unit_nr > 0)
THEN
227 WRITE (unit=unit_nr, fmt=
"(T2,A)") &
228 "MO's have a uniform occupation pattern"
230 WRITE (unit=unit_nr, fmt=
"(T2,A)") &
231 "MO's have varying occupations"
237 iao_env%do_iao = .true.
238 iao_env%do_charges = .true.
239 IF (skip_localize)
THEN
240 iao_env%do_bondorbitals = .false.
241 iao_env%do_center = .false.
243 iao_env%do_bondorbitals = .true.
244 iao_env%do_center = .true.
246 iao_env%eps_occ = 1.0e-4_dp
248 iao_env%pos_periodic = .NOT. all(cell%perd == 0)
256 IF (ref_orb_canonical)
THEN
257 CALL get_mo_set(mos_loc(ispin), mo_coeff=moref)
263 ALLOCATE (bcenter(5, no, nspin))
266 c_iao_coef=c_iao_coef, mos=mos_loc, bond_centers=bcenter)
269 IF (iao_env%do_bondorbitals .AND. iao_env%do_center)
THEN
275 smat => matrix_s(1)%matrix
276 ALLOCATE (rotmat(nspin))
278 CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
284 template_fmstruct=cloc%matrix_struct)
289 CALL parallel_gemm(
'T',
'N', norb, norb, nao, 1.0_dp, moref, &
290 smo, 0.0_dp, rotmat(ispin))
296 ALLOCATE (fij_mat(nspin))
298 CALL cp_fm_create(fij_mat(ispin), rotmat(ispin)%matrix_struct)
302 CALL get_mo_set(mos(ispin), nmo=norb, occupation_numbers=occupation_numbers)
307 CALL parallel_gemm(
'N',
'N', norb, norb, norb, 1.0_dp, fij_mat(ispin), &
308 rotmat(ispin), 0.0_dp, smo)
309 CALL parallel_gemm(
'T',
'N', norb, norb, norb, 1.0_dp, rotmat(ispin), &
310 smo, 0.0_dp, fij_mat(ispin))
316 ALLOCATE (ciao(nspin))
318 CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
324 template_fmstruct=cloc%matrix_struct)
329 CALL parallel_gemm(
'T',
'N', nref, norb, nao, 1.0_dp, c_iao_coef(ispin), &
330 smo, 0.0_dp, ciao(ispin))
335 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
336 nkind =
SIZE(qs_kind_set)
337 ALLOCATE (ref_basis_set_list(nkind))
339 qs_kind => qs_kind_set(ikind)
340 NULLIFY (ref_basis_set_list(ikind)%gto_basis_set)
342 CALL get_qs_kind(qs_kind=qs_kind, basis_set=refbasis, basis_type=
"MIN")
343 IF (
ASSOCIATED(refbasis)) ref_basis_set_list(ikind)%gto_basis_set => refbasis
345 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, natom=natom)
346 ALLOCATE (refbas_blk_sizes(natom))
348 basis=ref_basis_set_list)
351 ALLOCATE (orb_weight(nspin))
352 ALLOCATE (mcharge(natom, 1))
354 CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
357 template_fmstruct=cloc%matrix_struct)
363 CALL dbcsr_get_info(smat, distribution=dbcsr_dist)
364 CALL dbcsr_create(matrix=dmat, name=
"DMAT", &
365 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
366 row_blk_size=refbas_blk_sizes, col_blk_size=refbas_blk_sizes, &
368 CALL dbcsr_reserve_diag_blocks(dmat)
372 template_fmstruct=ciao(1)%matrix_struct)
379 mo_coeff=cloc, nmo=norb, &
380 occupation_numbers=occupation_numbers)
384 CALL cp_fm_to_fm(ciao(ispin), cvec, ncol=1, source_start=iorb, target_start=1)
385 focc = occupation_numbers(iorb)
386 CALL dbcsr_set(dmat, 0.0_dp)
390 start_col=iorb, n_rows=natom, n_cols=1)
391 checksum = sum(mcharge(:, 1))
392 IF (abs(focc - checksum) > 1.e-6_dp)
THEN
393 CALL cp_warn(__location__,
"Sum of atomic orbital weights is incorrect")
394 IF (unit_nr > 0)
THEN
395 WRITE (unit=unit_nr, fmt=
"(T2,A,F10.6,T40,A,F10.6)") &
396 "Orbital occupation:", focc, &
397 "Sum of atomic orbital weights:", checksum
403 ALLOCATE (odiag(norb))
406 IF (odiag(iorb) < 1.e-8_dp) cycle
407 CALL dbcsr_set(dmat, 0.0_dp)
408 CALL cp_fm_to_fm(ciao(ispin), cvec, ncol=1, source_start=iorb, target_start=1)
411 IF (focc < 1.e-12_dp) cycle
412 CALL cp_fm_to_fm(ciao(ispin), cvec2, ncol=1, source_start=jorb, target_start=1)
416 checksum = sum(mcharge(:, 1))
418 IF (abs(focc - checksum) > 1.e-6_dp)
THEN
419 CALL cp_warn(__location__,
"Sum of atomic orbital weights is incorrect")
420 IF (unit_nr > 0)
THEN
421 WRITE (unit=unit_nr, fmt=
"(T2,A,F10.6,T40,A,F10.6)") &
422 "Orbital occupation:", focc, &
423 "Sum of atomic orbital weights:", checksum
426 mcharge(:, 1) = mcharge(:, 1)/focc
428 start_col=iorb, n_rows=natom, n_cols=1)
434 CALL dbcsr_release(dmat)
439 ALLOCATE (atener(natom), ateks(natom), atecc(natom))
440 ALLOCATE (ate1xc(natom), ate1h(natom), atewald(natom))
447 IF (detailed_ener)
THEN
448 ALLOCATE (atdet(natom, 10), atmul(natom, 10), amval(natom))
453 CALL dbcsr_create(dkmat, template=smat)
454 CALL dbcsr_copy(dkmat, smat)
455 CALL dbcsr_set(dkmat, 0.0_dp)
457 CALL get_qs_env(qs_env, matrix_h=matrix_h, matrix_ks=matrix_ks, kinetic=matrix_t)
458 ALLOCATE (ks_mat(nspin), vhxc_mat(nspin), exc_mat(1))
460 ALLOCATE (ks_mat(ispin)%matrix)
461 CALL dbcsr_create(ks_mat(ispin)%matrix, template=matrix_h(1)%matrix)
462 CALL dbcsr_copy(ks_mat(ispin)%matrix, matrix_h(1)%matrix)
463 CALL dbcsr_add(ks_mat(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp)
464 CALL dbcsr_scale(ks_mat(ispin)%matrix, 0.5_dp)
466 ALLOCATE (vhxc_mat(ispin)%matrix)
467 CALL dbcsr_create(vhxc_mat(ispin)%matrix, template=smat)
468 CALL dbcsr_copy(vhxc_mat(ispin)%matrix, smat)
469 CALL dbcsr_set(vhxc_mat(ispin)%matrix, 0.0_dp)
471 ALLOCATE (exc_mat(1)%matrix)
472 CALL dbcsr_create(exc_mat(1)%matrix, template=smat)
473 CALL dbcsr_copy(exc_mat(1)%matrix, smat)
474 CALL dbcsr_set(exc_mat(1)%matrix, 0.0_dp)
476 CALL vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
478 IF (ewald_correction)
THEN
479 ALLOCATE (dve_mat%matrix)
480 CALL dbcsr_create(dve_mat%matrix, template=matrix_h(1)%matrix)
481 CALL dbcsr_copy(dve_mat%matrix, matrix_h(1)%matrix)
482 CALL dbcsr_set(dve_mat%matrix, 0.0_dp)
483 CALL vh_ewald_correction(qs_env, ealpha, dve_mat, atewald)
487 CALL dbcsr_add(ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix, 1.0_dp, 1.0_dp)
490 IF (detailed_ener .AND. do_hfx)
THEN
491 ALLOCATE (matrix_hfx(nspin))
493 ALLOCATE (matrix_hfx(nspin)%matrix)
494 CALL dbcsr_create(matrix_hfx(ispin)%matrix, template=smat)
495 CALL dbcsr_copy(matrix_hfx(ispin)%matrix, smat)
496 CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
498 CALL get_hfx_ks_matrix(qs_env, matrix_hfx)
502 CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc, nmo=norb)
503 ALLOCATE (mweight(1, norb))
506 start_col=1, n_rows=1, n_cols=norb)
512 CALL dbcsr_dot(dkmat, ks_mat(ispin)%matrix, ecc)
513 ateks(iatom) = ateks(iatom) + ecc
515 IF (detailed_ener)
THEN
516 CALL dbcsr_dot(dkmat, matrix_t(1)%matrix, e1)
517 atdet(iatom, 1) = atdet(iatom, 1) + e1
518 CALL dbcsr_dot(dkmat, matrix_h(1)%matrix, e2)
519 atdet(iatom, 2) = atdet(iatom, 2) + (e2 - e1)
521 CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 0.5_dp)
522 CALL dbcsr_dot(dkmat, matrix_hfx(ispin)%matrix, ehfx)
523 atdet(iatom, 5) = atdet(iatom, 5) + ehfx
524 CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 2.0_dp)
528 CALL dbcsr_dot(dkmat, exc_mat(1)%matrix, e1)
529 atdet(iatom, 3) = atdet(iatom, 3) + ecc - e2 - e1 - ehfx
530 atdet(iatom, 4) = atdet(iatom, 4) + e1
532 IF (ewald_correction)
THEN
533 CALL dbcsr_dot(dkmat, dve_mat%matrix, ecc)
534 atewald(iatom) = atewald(iatom) + ecc
541 IF (detailed_ener)
THEN
543 CALL get_qs_env(qs_env, rho=rho, para_env=para_env)
549 atmul(1:natom, 1) = atmul(1:natom, 1) + amval(1:natom)
552 atmul(1:natom, 2) = atmul(1:natom, 2) + amval(1:natom)
553 atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
556 atmul(1:natom, 3) = atmul(1:natom, 3) + amval(1:natom)
560 atmul(1:natom, 5) = atmul(1:natom, 5) + 0.5_dp*amval(1:natom)
561 atmul(1:natom, 3) = atmul(1:natom, 3) - 0.5_dp*amval(1:natom)
565 atmul(1:natom, 4) = atmul(1:natom, 4) + amval(1:natom)
566 atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
568 atmul(1:natom, 2) = atmul(1:natom, 2) - atmul(1:natom, 1)
571 CALL dbcsr_release(dkmat)
573 CALL dbcsr_release(ks_mat(ispin)%matrix)
574 CALL dbcsr_release(vhxc_mat(ispin)%matrix)
575 DEALLOCATE (ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix)
578 CALL dbcsr_release(exc_mat(1)%matrix)
579 DEALLOCATE (exc_mat(1)%matrix)
580 DEALLOCATE (ks_mat, vhxc_mat, exc_mat)
582 DEALLOCATE (refbas_blk_sizes)
583 DEALLOCATE (ref_basis_set_list)
584 IF (detailed_ener .AND. do_hfx)
THEN
586 CALL dbcsr_release(matrix_hfx(ispin)%matrix)
587 DEALLOCATE (matrix_hfx(ispin)%matrix)
589 DEALLOCATE (matrix_hfx)
591 IF (ewald_correction)
THEN
592 CALL dbcsr_release(dve_mat%matrix)
593 DEALLOCATE (dve_mat%matrix)
606 atener(1:natom) = ateks(1:natom)
608 CALL group%sum(atecc)
609 atener(1:natom) = atener(1:natom) + atecc(1:natom)
610 IF (detailed_ener) atdet(1:natom, 6) = atdet(1:natom, 6) + atecc(1:natom)
612 CALL group%sum(ate1xc)
613 CALL group%sum(ate1h)
614 atener(1:natom) = atener(1:natom) + ate1xc(1:natom) + ate1h(1:natom)
615 IF (detailed_ener)
THEN
616 IF (gapw .OR. gapw_xc) atdet(1:natom, 8) = atdet(1:natom, 8) + ate1xc(1:natom)
617 IF (gapw) atdet(1:natom, 9) = atdet(1:natom, 9) + ate1h(1:natom)
620 atecc(1:natom) = 0.0_dp
622 CALL group%sum(atecc)
623 atener(1:natom) = atener(1:natom) + atecc(1:natom)
624 IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
625 IF (ewald_correction)
THEN
626 atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
629 atecc(1:natom) = 0.0_dp
631 CALL group%sum(atecc)
632 atener(1:natom) = atener(1:natom) + atecc(1:natom)
633 IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
634 IF (ewald_correction)
THEN
635 atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
639 atecc(1:natom) = 0.0_dp
640 CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
643 CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env)
644 IF (
ASSOCIATED(gcp_env))
THEN
647 CALL group%sum(atecc)
648 atener(1:natom) = atener(1:natom) + atecc(1:natom)
649 IF (detailed_ener) atdet(1:natom, 10) = atdet(1:natom, 10) + atecc(1:natom)
652 ekts = energy%kts/real(natom, kind=
dp)
653 atener(1:natom) = atener(1:natom) + ekts
655 IF (detailed_ener)
THEN
656 IF (unit_nr > 0)
THEN
657 WRITE (unit_nr, fmt=
"(/,T2,A)")
"Detailed IAO Atomic Energy Components "
658 CALL write_atdet(unit_nr, atdet(:, 1),
"Kinetic Energy")
659 CALL write_atdet(unit_nr, atdet(:, 2),
"Short-Range Core and/or PP Energy")
660 CALL write_atdet(unit_nr, atdet(:, 3),
"Hartree Energy (long-ranged part)")
661 CALL write_atdet(unit_nr, atdet(:, 5),
"Exact Exchange Energy")
662 CALL write_atdet(unit_nr, atdet(:, 4),
"Exchange-Correlation Energy")
663 CALL write_atdet(unit_nr, atdet(:, 6),
"Atomic Core Hartree: Int(nc V(n+nc) dx")
664 CALL write_atdet(unit_nr, atdet(:, 7),
"Core Self Energy Corrections")
666 CALL write_atdet(unit_nr, atdet(:, 9),
"GAPW: 1-center Hartree Energy")
668 IF (gapw .OR. gapw_xc)
THEN
669 CALL write_atdet(unit_nr, atdet(:, 8),
"GAPW: 1-center XC Energy")
671 IF (abs(evdw) > 1.e-14_dp)
THEN
672 CALL write_atdet(unit_nr, atdet(:, 10),
"vdW Pairpotential Energy")
675 atecc(iatom) = sum(atdet(iatom, 1:10)) - atener(iatom)
677 CALL write_atdet(unit_nr, atecc(:),
"Missing Energy")
679 WRITE (unit_nr, fmt=
"(/,T2,A)")
"Detailed Mulliken Atomic Energy Components "
680 CALL write_atdet(unit_nr, atmul(:, 1),
"Kinetic Energy")
681 CALL write_atdet(unit_nr, atmul(:, 2),
"Short-Range Core and/or PP Energy")
682 CALL write_atdet(unit_nr, atmul(:, 3),
"Hartree Energy (long-ranged part)")
683 CALL write_atdet(unit_nr, atmul(:, 5),
"Exact Exchange Energy")
684 CALL write_atdet(unit_nr, atmul(:, 4),
"Exchange-Correlation Energy")
685 CALL write_atdet(unit_nr, atdet(:, 6),
"Atomic Core Hartree: Int(nc V(n+nc) dx")
686 CALL write_atdet(unit_nr, atdet(:, 7),
"Core Self Energy Corrections")
688 CALL write_atdet(unit_nr, atdet(:, 9),
"GAPW: 1-center Hartree Energy")
690 IF (gapw .OR. gapw_xc)
THEN
691 CALL write_atdet(unit_nr, atdet(:, 8),
"GAPW: 1-center XC Energy")
693 IF (abs(evdw) > 1.e-14_dp)
THEN
694 CALL write_atdet(unit_nr, atdet(:, 10),
"vdW Pairpotential Energy")
699 IF (unit_nr > 0)
THEN
702 CALL write_atener(unit_nr, particle_set, atener,
"Atomic Energy Decomposition")
703 sum_energy = sum(atener(:))
704 checksum = abs(e_pot - sum_energy)
705 WRITE (unit=unit_nr, fmt=
"((T2,A,T56,F25.13))") &
706 "Potential energy (Atomic):", sum_energy, &
707 "Potential energy (Total) :", e_pot, &
708 "Difference :", checksum
709 cpassert((checksum < ateps*abs(e_pot)))
711 IF (ewald_correction)
THEN
712 WRITE (unit=unit_nr, fmt=
"(/,(T2,A,F10.3,A))")
"Universal Ewald Parameter: ", ealpha,
" [Angstrom]"
713 CALL write_atener(unit_nr, particle_set, atewald,
"Change in Atomic Energy Decomposition")
714 sum_energy = sum(atewald(:))
715 WRITE (unit=unit_nr, fmt=
"((T2,A,T56,F25.13))") &
716 "Total Change in Potential energy:", sum_energy
720 IF (detailed_ener)
THEN
721 DEALLOCATE (atdet, atmul, amval)
724 IF (unit_nr > 0)
THEN
725 WRITE (unit=unit_nr, fmt=
"(/,T2,A)") &
726 "!--------------------------- END OF ED ANALYSIS ------------------------------!"
729 DEALLOCATE (atener, ateks, atecc, ate1xc, ate1h, atewald)
731 CALL timestop(handle)
744 SUBROUTINE vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
746 TYPE(dbcsr_p_type),
DIMENSION(:) :: vhxc_mat, exc_mat
747 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: atecc, ate1xc, ate1h
749 CHARACTER(len=*),
PARAMETER :: routinen =
'vhxc_correction'
751 INTEGER :: handle, iatom, ispin, natom, nspins
752 LOGICAL :: do_admm_corr, gapw, gapw_xc
753 REAL(kind=
dp) :: eh1, exc1
756 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p
767 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
774 CALL timeset(routinen, handle)
776 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env)
778 nspins = dft_control%nspins
779 gapw = dft_control%qs_control%gapw
780 gapw_xc = dft_control%qs_control%gapw_xc
781 do_admm_corr = .false.
782 IF (dft_control%do_admm)
THEN
785 IF (do_admm_corr)
THEN
787 xc_section => admm_env%xc_section_aux
794 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
795 CALL auxbas_pw_pool%create_pw(xc_den)
796 ALLOCATE (vxc(nspins))
798 CALL auxbas_pw_pool%create_pw(vxc(ispin))
800 IF (needs%tau .OR. needs%tau_spin)
THEN
801 ALLOCATE (vtau(nspins))
803 CALL auxbas_pw_pool%create_pw(vtau(ispin))
808 CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
809 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env, atecc=atecc)
810 IF (gapw .OR. gapw_xc)
THEN
811 CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, &
812 rho_atom_set=rho_atom_set, ecoul_1c=ecoul_1c, &
813 natom=natom, para_env=para_env)
817 CALL vh_1c_gg_integrals(qs_env, eh1, ecoul_1c, local_rho_set, para_env, tddft=.false.)
818 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
820 local_rho_set=local_rho_set, atener=atecc)
825 CALL get_qs_env(qs_env, rho_xc=rho_struct, dispersion_env=dispersion_env)
827 CALL get_qs_env(qs_env, rho=rho_struct, dispersion_env=dispersion_env)
830 cpassert(.NOT. do_admm_corr)
832 IF (needs%tau .OR. needs%tau_spin)
THEN
833 CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
834 xc_den=xc_den, vxc=vxc, vtau=vtau)
836 CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
837 xc_den=xc_den, vxc=vxc)
841 CALL pw_axpy(xc_den, vxc(ispin))
842 CALL pw_scale(vxc(ispin), vxc(ispin)%pw_grid%dvol)
843 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vxc(ispin), &
844 hmat=vhxc_mat(ispin), calculate_forces=.false., &
845 gapw=(gapw .OR. gapw_xc))
846 IF (needs%tau .OR. needs%tau_spin)
THEN
847 CALL pw_scale(vtau(ispin), -0.5_dp*vtau(ispin)%pw_grid%dvol)
848 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vtau(ispin), &
849 hmat=vhxc_mat(ispin), calculate_forces=.false., &
850 compute_tau=.true., gapw=(gapw .OR. gapw_xc))
853 CALL pw_scale(xc_den, xc_den%pw_grid%dvol)
854 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=xc_den, &
855 hmat=exc_mat(1), calculate_forces=.false., &
856 gapw=(gapw .OR. gapw_xc))
858 IF (gapw .OR. gapw_xc)
THEN
860 CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
861 CALL update_ks_atom(qs_env, vhxc_mat, matrix_p, forces=.false., kscale=-0.5_dp)
864 ate1xc(iatom) = ate1xc(iatom) + &
865 rho_atom_set(iatom)%exc_h - rho_atom_set(iatom)%exc_s
869 ate1h(iatom) = ate1h(iatom) + &
870 ecoul_1c(iatom)%ecoul_1_h - ecoul_1c(iatom)%ecoul_1_s + &
871 ecoul_1c(iatom)%ecoul_1_z - ecoul_1c(iatom)%ecoul_1_0
876 CALL auxbas_pw_pool%give_back_pw(xc_den)
878 CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
880 IF (needs%tau .OR. needs%tau_spin)
THEN
882 CALL auxbas_pw_pool%give_back_pw(vtau(ispin))
886 CALL timestop(handle)
888 END SUBROUTINE vhxc_correction
897 SUBROUTINE vh_ewald_correction(qs_env, ealpha, vh_mat, atewd)
899 REAL(kind=
dp),
INTENT(IN) :: ealpha
900 TYPE(dbcsr_p_type) :: vh_mat
901 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: atewd
903 CHARACTER(len=*),
PARAMETER :: routinen =
'vh_ewald_correction'
905 INTEGER :: handle, ikind, ispin, natom, nkind
906 REAL(kind=
dp) :: eps_core_charge, rhotot, zeff
907 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: alpha, atecc, ccore
911 POINTER :: sab_orb, sac_ae
920 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
923 CALL timeset(routinen, handle)
926 ALLOCATE (atecc(natom))
927 CALL get_qs_env(qs_env=qs_env, nkind=nkind, qs_kind_set=qs_kind_set, &
928 dft_control=dft_control)
929 ALLOCATE (alpha(nkind), ccore(nkind))
932 alpha(ikind) = ealpha
933 ccore(ikind) = zeff*(ealpha/
pi)**1.5_dp
937 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
938 poisson_env=poisson_env)
940 CALL auxbas_pw_pool%create_pw(v_hewd_gspace)
941 CALL auxbas_pw_pool%create_pw(v_hewd_rspace)
942 CALL auxbas_pw_pool%create_pw(rho_tot_ewd_g)
943 CALL auxbas_pw_pool%create_pw(rho_tot_ewd_r)
949 DO ispin = 1, dft_control%nspins
950 CALL pw_axpy(rho_r(ispin), rho_tot_ewd_r)
954 CALL pw_poisson_solve(poisson_env, density=rho_tot_ewd_g, vhartree=v_hewd_gspace)
956 CALL pw_scale(v_hewd_rspace, v_hewd_rspace%pw_grid%dvol)
958 CALL integrate_v_gaussian_rspace(v_hewd_rspace, qs_env, alpha, ccore, atecc)
959 atewd(1:natom) = atewd(1:natom) + atecc(1:natom)
962 CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
963 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env, atecc=atecc)
964 atewd(1:natom) = atewd(1:natom) - atecc(1:natom)
966 CALL pw_axpy(v_hartree_rspace, v_hewd_rspace, -1.0_dp)
967 CALL dbcsr_set(vh_mat%matrix, 0.0_dp)
968 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hewd_rspace, hmat=vh_mat, &
969 calculate_forces=.false.)
971 CALL dbcsr_scale(vh_mat%matrix, 0.5_dp)
974 eps_core_charge = dft_control%qs_control%eps_core_charge
975 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
976 particle_set=particle_set, sab_orb=sab_orb, sac_ppl=sac_ae)
977 CALL build_erfc(vh_mat, qs_kind_set, atomic_kind_set, particle_set, &
978 alpha, ccore, eps_core_charge, sab_orb, sac_ae)
979 CALL dbcsr_scale(vh_mat%matrix, -1.0_dp)
981 CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha(ikind), &
982 ccore_charge=ccore(ikind))
984 CALL build_erfc(vh_mat, qs_kind_set, atomic_kind_set, particle_set, &
985 alpha, ccore, eps_core_charge, sab_orb, sac_ae)
986 CALL dbcsr_scale(vh_mat%matrix, -1.0_dp)
988 CALL auxbas_pw_pool%give_back_pw(v_hewd_gspace)
989 CALL auxbas_pw_pool%give_back_pw(v_hewd_rspace)
990 CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_g)
991 CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_r)
993 DEALLOCATE (alpha, ccore)
995 CALL timestop(handle)
997 END SUBROUTINE vh_ewald_correction
1004 SUBROUTINE get_hfx_ks_matrix(qs_env, matrix_hfx)
1006 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_hfx
1008 CHARACTER(len=*),
PARAMETER :: routinen =
'get_hfx_ks_matrix'
1011 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p
1014 CALL timeset(routinen, handle)
1018 CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, update_energy=.false., &
1019 recalc_integrals=.false.)
1021 CALL timestop(handle)
1023 END SUBROUTINE get_hfx_ks_matrix
1034 SUBROUTINE write_atener(iounit, particle_set, atener, label)
1036 INTEGER,
INTENT(IN) :: iounit
1038 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: atener
1039 CHARACTER(LEN=*),
INTENT(IN) :: label
1043 IF (iounit > 0)
THEN
1044 WRITE (unit=iounit, fmt=
"(/,T2,A)") trim(label)
1045 WRITE (unit=iounit, fmt=
"(T4,A,T30,A,T42,A,T54,A,T69,A)") &
1046 "Atom Element",
"X",
"Y",
"Z",
"Energy[a.u.]"
1047 natom =
SIZE(atener)
1049 WRITE (unit=iounit, fmt=
"(I6,T12,A2,T24,3F12.6,F21.10)") i, &
1050 trim(adjustl(particle_set(i)%atomic_kind%element_symbol)), &
1051 particle_set(i)%r(1:3)*
angstrom, atener(i)
1053 WRITE (unit=iounit, fmt=
"(A)")
""
1056 END SUBROUTINE write_atener
1067 SUBROUTINE write_atdet(iounit, atener, label)
1068 INTEGER,
INTENT(IN) :: iounit
1069 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: atener
1070 CHARACTER(LEN=*),
INTENT(IN) :: label
1074 IF (iounit > 0)
THEN
1075 natom =
SIZE(atener)
1076 WRITE (unit=iounit, fmt=
"(T2,A)") trim(label)
1077 WRITE (unit=iounit, fmt=
"(5G16.8)") atener(1:natom)
1080 END SUBROUTINE write_atdet
Types and set/get functions for auxiliary density matrix methods.
Define the atomic kind types and their sub types.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public eriksen2020
Handles all functions related to the CELL.
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
subroutine, public build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, calpha, ccore, eps_core_charge, sab_orb, sac_ae)
Integrals = -Z*erfc(a*r)/r.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Calculate Energy Decomposition analysis.
subroutine, public edmf_analysis(qs_env, input_section, unit_nr)
...
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...
Utilities for hfx and admm methods.
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data)
Add the hfx contributions to the Hamiltonian.
Calculate intrinsic atomic orbitals and analyze wavefunctions.
subroutine, public iao_wfn_analysis(qs_env, iao_env, unit_nr, c_iao_coef, mos, bond_centers)
...
subroutine, public print_center_spread(moments, nocc, print_section, iounit)
prints charge center and spreads for all orbitals
subroutine, public iao_charges(p_matrix, charges)
compute the atomic charges (orthogonal basis)
Calculate ntrinsic atomic orbitals and analyze wavefunctions.
subroutine, public iao_set_default(iao_env)
...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
basic linear algebra operations for full matrixes
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public angstrom
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
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Calculation of the energies concerning the core charge distribution.
subroutine, public calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, e_overlap_core, atecc)
Calculate the overlap energy of the core charge distribution.
subroutine, public calculate_ecore_self(qs_env, e_self_core, atecc)
Calculate the self energy of the core charge distribution.
subroutine, public calculate_ecore_alpha(qs_env, alpha, atecc)
Calculate the overlap and self energy of the core charge distribution for a given alpha Use a minimum...
Calculation of dispersion using pair potentials.
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition of disperson types for DFT calculations.
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.
Calculation of gCP pair potentials.
subroutine, public calculate_gcp_pairpot(qs_env, gcp_env, energy, calculate_forces, ategcp)
...
Definition of gCP types for DFT calculations.
Integrate single or product functions over a potential on a RS grid.
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.
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
collects routines that perform operations directly related to MOs
Definition and initialisation of the mo data type.
subroutine, public duplicate_mo_set(mo_set_new, mo_set_old)
allocate a new mo_set, and copy the old data
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Define the neighbor list data types and the corresponding functionality.
subroutine, public integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, local_rho_set_2nd, atener, kforce)
...
subroutine, public zero_rho_atom_integrals(rho_atom_set)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
subroutine, public calculate_vxc_atom(qs_env, energy_only, exc1, gradient_atom_set, adiabatic_rescale_factor, kind_set_external, rho_atom_set_external, xc_section_external)
...
subroutine, public qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, xc_ener, xc_den, vxc, vtau)
calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
stores some data used in wavefunction fitting
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
keeps the information about the structure of a full matrix
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...