comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Energy decomposition: correct detailed analysis (#3459)

Browse source
This commit is contained in:
Juerg Hutter 2024-06-11 14:18:47 +02:00 • committed by GitHub
parent 52714c9400
commit d2eb2ba5be
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
2 changed files with 44 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

9
src/ed_analysis.F
View file

@ -582,6 +582,8 @@ CONTAINS
atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
END DO
atmul(1:natom, 2) = atmul(1:natom, 2) - atmul(1:natom, 1)
atmul(1:natom, 3) = atmul(1:natom, 3) + 0.5_dp*atmul(1:natom, 2)
atmul(1:natom, 2) = 0.5_dp*atmul(1:natom, 2) + 0.5_dp*atcore(1:natom)
END IF
!
CALL dbcsr_release(dkmat)
@ -623,7 +625,7 @@ CONTAINS
! KS energy
atener(1:natom) = ateks(1:natom)
! 1/2 of VPP contribution Tr[VPP(K)*P]
CALL group%sum(atcore) ! maybe not needed
CALL group%sum(atcore)
atener(1:natom) = atener(1:natom) + 0.5_dp*atcore(1:natom)
! core energy corrections
CALL group%sum(atecc)
@ -673,7 +675,10 @@ CONTAINS
ekts = energy%kts/REAL(natom, KIND=dp)
atener(1:natom) = atener(1:natom) + ekts
! 0.5 Vpp(at)*D + 0.5 * Vpp*D(at)
IF (detailed_ener) atdet(1:natom, 2) = 0.5_dp*atdet(1:natom, 2) + 0.5_dp*atcore(1:natom)
IF (detailed_ener) THEN
atdet(1:natom, 3) = atdet(1:natom, 3) + 0.5_dp*atdet(1:natom, 2)
atdet(1:natom, 2) = 0.5_dp*atdet(1:natom, 2) + 0.5_dp*atcore(1:natom)
END IF
!
IF (detailed_ener) THEN
IF (unit_nr > 0) THEN

44
src/qs_energy_utils.F
View file

@ -40,6 +40,7 @@ MODULE qs_energy_utils
pw_scale
USE pw_pool_types, ONLY: pw_pool_type
USE pw_types, ONLY: pw_r3d_rs_type
USE qs_core_hamiltonian, ONLY: core_matrices
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
@ -227,25 +228,54 @@ CONTAINS
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER :: ispin
INTEGER :: ispin, natom, nspin
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: atcore
TYPE(atprop_type), POINTER :: atprop
TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), &
TARGET :: core_mat
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_ks, rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: math, matp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_rho_type), POINTER :: rho
NULLIFY (atprop, matrix_h, matrix_ks, rho, rho_ao)
CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_h=matrix_h, &
rho=rho, atprop=atprop)
CALL qs_rho_get(rho, rho_ao=rho_ao)
CALL get_qs_env(qs_env=qs_env, atprop=atprop)
IF (atprop%energy) THEN
CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_h=matrix_h, rho=rho)
CALL qs_rho_get(rho, rho_ao=rho_ao)
! E = 0.5*Tr(H*P+F*P)
atprop%atener = 0._dp
DO ispin = 1, SIZE(rho_ao)
nspin = SIZE(rho_ao)
DO ispin = 1, nspin
CALL atom_trace(matrix_h(1)%matrix, rho_ao(ispin)%matrix, &
0.5_dp, atprop%atener)
CALL atom_trace(matrix_ks(ispin)%matrix, rho_ao(ispin)%matrix, &
0.5_dp, atprop%atener)
END DO
!
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
IF (.NOT. dft_control%qs_control%semi_empirical .AND. &
.NOT. dft_control%qs_control%xtb .AND. &
.NOT. dft_control%qs_control%dftb) THEN
CALL get_qs_env(qs_env=qs_env, natom=natom)
ALLOCATE (atcore(natom))
ALLOCATE (core_mat(1))
ALLOCATE (core_mat(1)%matrix)
CALL dbcsr_create(core_mat(1)%matrix, template=matrix_h(1)%matrix)
CALL dbcsr_copy(core_mat(1)%matrix, matrix_h(1)%matrix)
CALL dbcsr_set(core_mat(1)%matrix, 0.0_dp)
math(1:1, 1:1) => core_mat(1:1)
matp(1:nspin, 1:1) => rho_ao(1:nspin)
CALL core_matrices(qs_env, math, matp, .FALSE., 0, atcore=atcore)
atprop%atener = atprop%atener + 0.5_dp*atcore
DO ispin = 1, nspin
CALL atom_trace(core_mat(1)%matrix, rho_ao(ispin)%matrix, &
-0.5_dp, atprop%atener)
END DO
DEALLOCATE (atcore)
CALL dbcsr_release(core_mat(1)%matrix)
DEALLOCATE (core_mat(1)%matrix)
DEALLOCATE (core_mat)
END IF
END IF
END SUBROUTINE qs_energies_mulliken