From 192a6345de7472f54758ae6d7e1be38c2906fc9c Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 6 Apr 2020 19:40:42 +0200 Subject: [PATCH] GIC functional --- include/parameters.h | 1 - input/basis | 48 ++- input/dft | 4 +- input/weight | 48 ++- src/QuAcK/G0T0.f90 | 69 +--- src/QuAcK/excitation_density_Tmatrix.f90 | 302 +++++++++--------- src/QuAcK/excitation_density_Tmatrix_so.f90 | 82 +++++ src/QuAcK/linear_response_A_matrix.f90 | 2 - src/QuAcK/linear_response_B_matrix.f90 | 2 - src/QuAcK/renormalization_factor_Tmatrix.f90 | 8 +- src/QuAcK/self_energy_Tmatrix_diag.f90 | 12 +- src/QuAcK/soG0T0.f90 | 4 +- src/eDFT/Makefile | 2 +- ..._lda_exchange_derivative_discontinuity.f90 | 65 ++++ src/eDFT/RGIC_lda_exchange_energy.f90 | 52 +++ .../RGIC_lda_exchange_individual_energy.f90 | 61 ++++ src/eDFT/RGIC_lda_exchange_potential.f90 | 61 ++++ ..._lda_exchange_derivative_discontinuity.f90 | 10 +- src/eDFT/RMFL20_lda_exchange_energy.f90 | 2 +- .../RMFL20_lda_exchange_individual_energy.f90 | 6 +- src/eDFT/fock_exchange_individual_energy.f90 | 6 +- src/eDFT/hartree_individual_energy.f90 | 10 +- .../lda_exchange_derivative_discontinuity.f90 | 4 + src/eDFT/lda_exchange_energy.f90 | 10 +- src/eDFT/lda_exchange_individual_energy.f90 | 4 + src/eDFT/lda_exchange_potential.f90 | 10 +- 26 files changed, 602 insertions(+), 283 deletions(-) create mode 100644 src/QuAcK/excitation_density_Tmatrix_so.f90 create mode 100644 src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 create mode 100644 src/eDFT/RGIC_lda_exchange_energy.f90 create mode 100644 src/eDFT/RGIC_lda_exchange_individual_energy.f90 create mode 100644 src/eDFT/RGIC_lda_exchange_potential.f90 diff --git a/include/parameters.h b/include/parameters.h index ee9572e..279fe1e 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -18,6 +18,5 @@ double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) -! double precision,parameter :: Cx1 = - 0.913d0*(4d0/3d0)*(1d0/pi)**(1d0/3d0) double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) diff --git a/input/basis b/input/basis index a1f7a8d..dc1936c 100644 --- a/input/basis +++ b/input/basis @@ -1,27 +1,43 @@ -1 5 +1 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 +2 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 diff --git a/input/dft b/input/dft index ca731a0..7679b2d 100644 --- a/input/dft +++ b/input/dft @@ -2,11 +2,11 @@ GOK-RKS # exchange rung: # Hartree = 0 -# LDA = 1: RS51,RMFL20 +# LDA = 1: RS51,RMFL20,RGIC # GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 666 HF + 1 RGIC # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 diff --git a/input/weight b/input/weight index a1f7a8d..dc1936c 100644 --- a/input/weight +++ b/input/weight @@ -1,27 +1,43 @@ -1 5 +1 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 +2 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index c74aca7..c0e8086 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -45,6 +45,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m double precision,allocatable :: X2s(:,:),X2t(:,:) double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) + double precision,allocatable :: rho1st(:,:,:),rho2st(:,:,:) double precision,allocatable :: SigT(:) double precision,allocatable :: Z(:) @@ -81,6 +82,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), & + rho1st(nBas,nO,nVVt),rho2st(nBas,nV,nOOt), & SigT(nBas),Z(nBas),eG0T0(nBas)) !---------------------------------------------- @@ -100,14 +102,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s(:)) call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s(:)) -! Compute excitation densities for the T-matrix - - rho1s(:,:,:) = 0d0 - rho2s(:,:,:) = 0d0 - - call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & - X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) - !---------------------------------------------- ! Triplet manifold !---------------------------------------------- @@ -125,72 +119,33 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t(:)) call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t(:)) -! Compute excitation densities for the T-matrix +!----------------------------------------------------------- +! Compute excitation densities for the T-matrix self-energy +!----------------------------------------------------------- - rho1t(:,:,:) = 0d0 - rho2t(:,:,:) = 0d0 - - call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & + call excitation_density_Tmatrix(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,ERI(:,:,:,:),rho1st(:,:,:),rho2st(:,:,:), & + X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:), & X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- - SigT(:) = 0d0 - rho2s(:,:,:) = 0d0 - rho2t(:,:,:) = 0d0 +! rho2s(:,:,:) = 0d0 +! rho2t(:,:,:) = 0d0 +! rho2st(:,:,:) = 0d0 call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & - SigT(:)) -!---------------------------------------------- -! Singlet manifold -!---------------------------------------------- - - ispin = 1 - -! Compute excitation densities for the T-matrix - - rho1s(:,:,:) = 0d0 - rho2s(:,:,:) = 0d0 - - call excitation_density_Tmatrix(ispin,0d0,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & - X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) - -!---------------------------------------------- -! Triplet manifold -!---------------------------------------------- - - ispin = 2 - -! Compute excitation densities for the T-matrix - - rho1t(:,:,:) = 0d0 - rho2t(:,:,:) = 0d0 - - call excitation_density_Tmatrix(ispin,0d0,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & - X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) - -!---------------------------------------------- -! Compute T-matrix version of the self-energy -!---------------------------------------------- - - rho2s(:,:,:) = 0d0 - rho2t(:,:,:) = 0d0 - - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & - Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & - Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & - SigT(:)) + rho1st(:,:,:),rho2st(:,:,:),SigT(:)) ! Compute renormalization factor for T-matrix self-energy call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & - Z(:)) + rho1st(:,:,:),rho2st(:,:,:),Z(:)) !---------------------------------------------- ! Solve the quasi-particle equation diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index fd900fc..f5a9e26 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -1,4 +1,5 @@ -subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) +subroutine excitation_density_Tmatrix(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,ERI,rho1st,rho2st, & + X1s,Y1s,rho1s,X2s,Y2s,rho2s,X1t,Y1t,rho1t,X2t,Y2t,rho2t) ! Compute excitation densities for T-matrix self-energy @@ -6,14 +7,24 @@ subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y ! Input variables - integer,intent(in) :: ispin - double precision,intent(in) :: db - integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: X1(nVV,nVV) - double precision,intent(in) :: Y1(nOO,nVV) - double precision,intent(in) :: X2(nVV,nOO) - double precision,intent(in) :: Y2(nOO,nOO) + integer,intent(in) :: nOOs + integer,intent(in) :: nOOt + double precision,intent(in) :: X1s(nVVs,nVVs) + double precision,intent(in) :: Y1s(nOOs,nVVs) + double precision,intent(in) :: X2s(nVVs,nOOs) + double precision,intent(in) :: Y2s(nOOs,nOOs) + integer,intent(in) :: nVVs + integer,intent(in) :: nVVt + double precision,intent(in) :: X1t(nVVt,nVVt) + double precision,intent(in) :: Y1t(nOOt,nVVt) + double precision,intent(in) :: X2t(nVVt,nOOt) + double precision,intent(in) :: Y2t(nOOt,nOOt) ! Local variables @@ -25,187 +36,176 @@ subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y ! Output variables - double precision,intent(out) :: rho1(nBas,nO,nVV) - double precision,intent(out) :: rho2(nBas,nV,nOO) + double precision,intent(out) :: rho1s(nBas,nO,nVVs) + double precision,intent(out) :: rho2s(nBas,nV,nOOs) + double precision,intent(out) :: rho1t(nBas,nO,nVVt) + double precision,intent(out) :: rho2t(nBas,nV,nOOt) + double precision,intent(out) :: rho1st(nBas,nO,nVVt) + double precision,intent(out) :: rho2st(nBas,nV,nOOt) + +! Initialization + + rho1s(:,:,:) = 0d0 + rho2s(:,:,:) = 0d0 + rho1t(:,:,:) = 0d0 + rho2t(:,:,:) = 0d0 + rho1st(:,:,:) = 0d0 + rho2st(:,:,:) = 0d0 !---------------------------------------------- ! Singlet manifold !---------------------------------------------- - if(ispin == 1) then + do p=nC+1,nBas-nR - do p=nC+1,nBas-nR + do i=nC+1,nO + do ab=1,nVVs - do i=nC+1,nO - do ab=1,nVV - - cd = 0 - do c=nO+1,nBas-nR - do d=c,nBas-nR - cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) & -! + db*(ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) & -! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) & - + 0d0*db*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) & - + (1d0 - db)*ERI(p,i,c,d)*X1(cd,ab) - end do + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 + rho1s(p,i,ab) = rho1s(p,i,ab) + 2.0d0*ERI(p,i,c,d)*X1s(cd,ab) end do - - kl = 0 - do k=nC+1,nO - do l=k,nO - kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) & -! + db*(ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) & -! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) & - + 0d0*db*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) & - + (1d0 - db)*ERI(p,i,k,l)*Y1(kl,ab) - end do - end do - end do - end do - do a=1,nV-nR - do ij=1,nOO - - cd = 0 - do c=nO+1,nBas-nR - do d=c,nBas-nR - cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) + db*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) - end do + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + rho1s(p,i,ab) = rho1s(p,i,ab) + 2.0d0*ERI(p,i,k,l)*Y1s(kl,ab) end do - - kl = 0 - do k=nC+1,nO - do l=k,nO - kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) + db*(ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) - end do - end do - end do - end do + end do end do - end if + do a=1,nV-nR + do ij=1,nOOs + + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 + rho2s(p,a,ij) = rho2s(p,a,ij) + 2.0d0*ERI(p,nO+a,c,d)*X2s(cd,ij) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + rho2s(p,a,ij) = rho2s(p,a,ij) + 2.0d0*ERI(p,nO+a,k,l)*Y2s(kl,ij) + end do + end do + + end do + end do + + end do !---------------------------------------------- ! Triplet manifold !---------------------------------------------- - if(ispin == 2) then + do p=nC+1,nBas-nR - do p=nC+1,nBas-nR + do i=nC+1,nO + do ab=1,nVVt - do i=nC+1,nO - do ab=1,nVV - - cd = 0 - do c=nO+1,nBas-nR - do d=c+1,nBas-nR - cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + 1.0d0*db*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) & - + (1d0-db)*0d0*(ERI(p,i,c,d))*X1(cd,ab) - end do + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho1t(p,i,ab) = rho1t(p,i,ab) + 1d0*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1t(cd,ab) end do - - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + 1.0d0*db*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) & - + (1d0-db)*0d0*(ERI(p,i,k,l))*Y1(kl,ab) - end do - end do - end do - end do - do a=1,nV-nR - do ij=1,nOO - - cd = 0 - do c=nO+1,nBas-nR - do d=c+1,nBas-nR - cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) + 1d0*(ERI(p,nO+a,c,d) - db*ERI(p,nO+a,d,c))*X2(cd,ij) - end do + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho1t(p,i,ab) = rho1t(p,i,ab) + 1d0*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1t(kl,ab) end do - - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) + 1d0*(ERI(p,nO+a,k,l) - db*ERI(p,nO+a,l,k))*Y2(kl,ij) - end do - end do - end do - end do + end do end do - end if + do a=1,nV-nR + do ij=1,nOOt -!---------------------------------------------- -! Spinorbital basis -!---------------------------------------------- - - if(ispin == 3) then - - do p=nC+1,nBas-nR - - do i=nC+1,nO - do ab=1,nVV - - cd = 0 - do c=nO+1,nBas-nR - do d=c+1,nBas-nR - cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) - end do + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho2t(p,a,ij) = rho2t(p,a,ij) + 1d0*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2t(cd,ij) end do - - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) - end do - end do - end do - end do - do a=1,nV-nR - do ij=1,nOO - - cd = 0 - do c=nO+1,nBas-nR - do d=c+1,nBas-nR - cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) - end do + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho2t(p,a,ij) = rho2t(p,a,ij) + 1d0*(ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2t(kl,ij) end do - - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) - end do - end do - end do - end do + end do end do - end if + end do + +!---------------------------------------------- +! Singlet-triplet crossed term +!---------------------------------------------- + + do p=nC+1,nBas-nR + + do i=nC+1,nO + do ab=1,nVVt + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho1st(p,i,ab) = rho1st(p,i,ab) + 2d0*ERI(p,i,c,d)*X1t(cd,ab) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho1st(p,i,ab) = rho1st(p,i,ab) + 2d0*ERI(p,i,k,l)*Y1t(kl,ab) + end do + end do + + end do + end do + + do a=1,nV-nR + do ij=1,nOOt + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho2st(p,a,ij) = rho2st(p,a,ij) + 2d0*ERI(p,nO+a,c,d)*X2t(cd,ij) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho2st(p,a,ij) = rho2st(p,a,ij) + 2d0*ERI(p,nO+a,k,l)*Y2t(kl,ij) + end do + end do + + end do + end do + + end do end subroutine excitation_density_Tmatrix diff --git a/src/QuAcK/excitation_density_Tmatrix_so.f90 b/src/QuAcK/excitation_density_Tmatrix_so.f90 new file mode 100644 index 0000000..0d18299 --- /dev/null +++ b/src/QuAcK/excitation_density_Tmatrix_so.f90 @@ -0,0 +1,82 @@ +subroutine excitation_density_Tmatrix_so(nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) + +! Compute excitation densities for T-matrix self-energy + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + integer :: p + integer :: ab,cd,ij,kl + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: rho1(nBas,nO,nVV) + double precision,intent(out) :: rho2(nBas,nV,nOO) + +! Initialization + + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + do p=nC+1,nBas-nR + + do i=nC+1,nO + do ab=1,nVV + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) + end do + end do + + end do + end do + + do a=1,nV-nR + do ij=1,nOO + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) + end do + end do + + end do + end do + + end do + +end subroutine excitation_density_Tmatrix_so diff --git a/src/QuAcK/linear_response_A_matrix.f90 b/src/QuAcK/linear_response_A_matrix.f90 index d2a7c9b..1a1d391 100644 --- a/src/QuAcK/linear_response_A_matrix.f90 +++ b/src/QuAcK/linear_response_A_matrix.f90 @@ -49,8 +49,6 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + (1d0 + delta_spin)*lambda*ERI(i,b,a,j) & - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a) -! + (1d0 + delta_spin)*lambda*ERI(i,j,a,b) & -! - (1d0 - delta_dRPA)*lambda*ERI(i,a,j,b) enddo enddo diff --git a/src/QuAcK/linear_response_B_matrix.f90 b/src/QuAcK/linear_response_B_matrix.f90 index 30e9c2b..5226059 100644 --- a/src/QuAcK/linear_response_B_matrix.f90 +++ b/src/QuAcK/linear_response_B_matrix.f90 @@ -46,8 +46,6 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_ B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) & - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a) -! B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,b,a,j) & -! - (1d0 - delta_dRPA)*lambda*ERI(i,a,b,j) enddo enddo diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90 index ce4be02..d6326c3 100644 --- a/src/QuAcK/renormalization_factor_Tmatrix.f90 +++ b/src/QuAcK/renormalization_factor_Tmatrix.f90 @@ -1,6 +1,6 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, & Omega1s,rho1s,Omega2s,rho2s,Omega1t,rho1t,Omega2t,rho2t, & - Z) + rho1st,rho2st,Z) ! Compute renormalization factor of the T-matrix self-energy @@ -15,9 +15,9 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV integer,intent(in) :: nVVs,nVVt double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt) - double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt) + double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt),rho1st(nBas,nO,nVVt) double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt) - double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt) + double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt),rho2st(nBas,nV,nOOt) ! Local variables @@ -69,6 +69,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV do cd=1,nVVt eps = e(p) + e(i) - Omega1t(cd) Z(p) = Z(p) + (rho1t(p,i,cd)/eps)**2 + Z(p) = Z(p) + (rho1st(p,i,cd)/eps)**2 enddo enddo enddo @@ -80,6 +81,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV do kl=1,nOOt eps = e(p) + e(nO+a) - Omega2t(kl) Z(p) = Z(p) + (rho2t(p,a,kl)/eps)**2 + Z(p) = Z(p) + (rho2st(p,a,kl)/eps)**2 enddo enddo enddo diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90 index 3882dd1..73e4ad1 100644 --- a/src/QuAcK/self_energy_Tmatrix_diag.f90 +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -1,6 +1,6 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, & Omega1s,rho1s,Omega2s,rho2s,Omega1t,rho1t,Omega2t,rho2t, & - SigT) + rho1st,rho2st,SigT) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -19,9 +19,9 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, integer,intent(in) :: nVVs,nVVt double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt) - double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt) + double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt),rho1st(nBas,nO,nVVt) double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt) - double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt) + double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt),rho2st(nBas,nV,nOOt) ! Local variables @@ -32,6 +32,10 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, double precision,intent(out) :: SigT(nBas) +! Initialization + + SigT(:) = 0d0 + !---------------------------------------------- ! Singlet part of the T-matrix self-energy !---------------------------------------------- @@ -69,6 +73,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do cd=1,nVVt eps = e(p) + e(i) - Omega1t(cd) SigT(p) = SigT(p) + rho1t(p,i,cd)**2/eps + SigT(p) = SigT(p) + rho1st(p,i,cd)**2/eps enddo enddo enddo @@ -80,6 +85,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do kl=1,nOOt eps = e(p) + e(nO+a) - Omega2t(kl) SigT(p) = SigT(p) + rho2t(p,a,kl)**2/eps + SigT(p) = SigT(p) + rho2st(p,a,kl)**2/eps enddo enddo enddo diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index 88793f9..79d2c27 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -91,14 +91,14 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) rho1(:,:,:) = 0d0 rho2(:,:,:) = 0d0 - call excitation_density_Tmatrix(ispin,1d0,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), & + call excitation_density_Tmatrix_so(nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), & X1(:,:),Y1(:,:),rho1(:,:,:),X2(:,:),Y2(:,:),rho2(:,:,:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- - rho2(:,:,:) = 0d0 +! rho2(:,:,:) = 0d0 call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & diff --git a/src/eDFT/Makefile b/src/eDFT/Makefile index 1024aba..3efafe3 100644 --- a/src/eDFT/Makefile +++ b/src/eDFT/Makefile @@ -3,7 +3,7 @@ BDIR =../../bin ODIR = obj OODIR = ../IntPak/obj SDIR =. -FC = gfortran -I$(IDIR) -Wall +FC = gfortran -I$(IDIR) ifeq ($(DEBUG),1) FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant else diff --git a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..cd0fa3f --- /dev/null +++ b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 @@ -0,0 +1,65 @@ +subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD) + +! Compute the restricted version of the GIC exchange individual energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + +! Local variables + + integer :: iEns,jEns + integer :: iG + double precision :: r + double precision :: dExdw(nEns) + double precision,external :: Kronecker_delta + + double precision :: a,b,c,w + double precision :: dCxGICdw + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + +! Compute correlation energy for ground- and doubly-excited states + + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 + + w = wEns(2) + dCxGICdw = CxLDA*(0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0))) + + dExdw(:) = 0d0 + + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + + if(r > threshold) then + + dExdw(1) = 0d0 + dExdw(2) = dExdw(2) + weight(iG)*dCxGICdw*r**(4d0/3d0) + + end if + + end do + + ExDD(:) = 0d0 + + do iEns=1,nEns + do jEns=2,nEns + + ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) + + end do + end do + +end subroutine RGIC_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/RGIC_lda_exchange_energy.f90 b/src/eDFT/RGIC_lda_exchange_energy.f90 new file mode 100644 index 0000000..dafc17b --- /dev/null +++ b/src/eDFT/RGIC_lda_exchange_energy.f90 @@ -0,0 +1,52 @@ +subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + +! Compute the restricted version of the GIC exchange functional + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: r + + double precision :: a,b,c,w + double precision :: CxGIC + +! Output variables + + double precision :: Ex + +! Weight-denepdent Cx coefficient + + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 + + w = wEns(2) + CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + +! Compute GIC-LDA exchange energy + + Ex = 0d0 + + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0) + Ex = Ex + weight(iG)*CxGIC*r**(4d0/3d0) + endif + + enddo + +end subroutine RGIC_lda_exchange_energy diff --git a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 new file mode 100644 index 0000000..68c553a --- /dev/null +++ b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 @@ -0,0 +1,61 @@ +subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex) + +! Compute the restricted version of the GIC exchange functional + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: CxGIC + double precision :: r,rI + double precision :: e_p,dedr + + double precision :: a,b,c,w + +! Output variables + + double precision,intent(out) :: Ex + +! Weight-dependent Cx coefficient for RMFL20 exchange functional + + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 + + w = wEns(2) + CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + +! Compute LDA exchange matrix in the AO basis + + Ex = 0d0 + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + rI = max(0d0,rho(iG)) + + if(r > threshold .and. rI > threshold) then + + e_p = CxLDA*r**(1d0/3d0) + dedr = 1d0/3d0*CxLDA*r**(-2d0/3d0) + Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) + + e_p = CxGIC*r**(1d0/3d0) + dedr = 1d0/3d0*CxGIC*r**(-2d0/3d0) + Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) + + endif + + enddo + +end subroutine RGIC_lda_exchange_individual_energy diff --git a/src/eDFT/RGIC_lda_exchange_potential.f90 b/src/eDFT/RGIC_lda_exchange_potential.f90 new file mode 100644 index 0000000..0b99764 --- /dev/null +++ b/src/eDFT/RGIC_lda_exchange_potential.f90 @@ -0,0 +1,61 @@ +subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + +! Compute the restricted version of the GIC exchange potential + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: mu,nu,iG + double precision :: r,vAO + + double precision :: a,b,c,w + double precision :: CxGIC + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Weight-dependent Cx coefficient for RMFL20 exchange functional + + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 + + w = wEns(2) + CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) + +! Compute LDA exchange matrix in the AO basis + + Fx(:,:) = 0d0 + + do mu=1,nBas + do nu=1,nBas + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) + Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxLDA*r**(1d0/3d0) + Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxGIC*r**(1d0/3d0) + + endif + + enddo + enddo + enddo + +end subroutine RGIC_lda_exchange_potential diff --git a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 index 3a8414c..15e1c3b 100644 --- a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 @@ -34,8 +34,10 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r r = max(0d0,rhow(iG)) if(r > threshold) then - dExdw(1) = dExdw(1) + weight(iG)*Cx0*r**(4d0/3d0) - dExdw(2) = dExdw(2) + weight(iG)*Cx1*r**(4d0/3d0) + + dExdw(1) = 0d0 + dExdw(2) = dExdw(2) + weight(iG)*(Cx1 - Cx0)*r**(4d0/3d0) + end if end do @@ -43,9 +45,9 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r ExDD(:) = 0d0 do iEns=1,nEns - do jEns=1,nEns + do jEns=2,nEns - ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*(dExdw(jEns) - dExdw(1)) + ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) end do end do diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 index 5e9cb1e..c9ca09f 100644 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_energy.f90 @@ -28,7 +28,7 @@ subroutine RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex ! Weight-denepdent Cx coefficient if(LDA_centered) then - Cxw = CxLDA + (Cx1 - Cx0)*wEns(2) + Cxw = CxLDA + (Cx1 - Cx0)*wEns(2) else Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 end if diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 index 6551b39..e582ac1 100644 --- a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 @@ -20,7 +20,7 @@ subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,we integer :: iG double precision :: Cxw double precision :: r,rI - double precision :: e,dedr + double precision :: e_p,dedr ! Output variables @@ -44,9 +44,9 @@ subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,we if(r > threshold .and. rI > threshold) then - e = Cxw*r**(1d0/3d0) + e_p = Cxw*r**(1d0/3d0) dedr = 1d0/3d0*Cxw*r**(-2d0/3d0) - Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) + Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) endif diff --git a/src/eDFT/fock_exchange_individual_energy.f90 b/src/eDFT/fock_exchange_individual_energy.f90 index e1839eb..9b6d3ae 100644 --- a/src/eDFT/fock_exchange_individual_energy.f90 +++ b/src/eDFT/fock_exchange_individual_energy.f90 @@ -24,8 +24,8 @@ subroutine fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) allocate(Fx(nBas,nBas)) - call fock_exchange_potential(nBas,P(:,:),ERI(:,:,:,:),Fx(:,:)) - - Ex = 0.5d0*trace_matrix(nBas,matmul(P(:,:),Fx(:,:))) !& + call fock_exchange_potential(nBas,Pw(:,:),ERI(:,:,:,:),Fx(:,:)) + Ex = trace_matrix(nBas,matmul(P(:,:),Fx(:,:))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),Fx(:,:))) end subroutine fock_exchange_individual_energy diff --git a/src/eDFT/hartree_individual_energy.f90 b/src/eDFT/hartree_individual_energy.f90 index 248c120..cf1dd3f 100644 --- a/src/eDFT/hartree_individual_energy.f90 +++ b/src/eDFT/hartree_individual_energy.f90 @@ -28,8 +28,9 @@ subroutine hartree_individual_energy(rung,nBas,ERI,J,Pw,P,EJ) case(0) - call hartree_coulomb(nBas,P(:,:),ERI(:,:,:,:),J(:,:)) - EJ = 0.5d0*trace_matrix(nBas,matmul(P(:,:),J(:,:))) + call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:)) + EJ = trace_matrix(nBas,matmul(P(:,:),J(:,:))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:))) ! LDA functionals @@ -57,8 +58,9 @@ subroutine hartree_individual_energy(rung,nBas,ERI,J,Pw,P,EJ) case(666) - call hartree_coulomb(nBas,P(:,:),ERI(:,:,:,:),J(:,:)) - EJ = 0.5d0*trace_matrix(nBas,matmul(P(:,:),J(:,:))) + call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:)) + EJ = trace_matrix(nBas,matmul(P(:,:),J(:,:))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:))) end select diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 index cb4527c..618b6ba 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -37,6 +37,10 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow call RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) + case ('RGIC') + + call RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) + case default call print_warning('!!! LDA exchange functional not available !!!') diff --git a/src/eDFT/lda_exchange_energy.f90 b/src/eDFT/lda_exchange_energy.f90 index 4a153df..d2f1616 100644 --- a/src/eDFT/lda_exchange_energy.f90 +++ b/src/eDFT/lda_exchange_energy.f90 @@ -23,24 +23,22 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) select case (DFA) -! Slater's LDA correlation functional - case ('S51') call S51_lda_exchange_energy(nGrid,weight,rho,Ex) -! Restricted version of Slater's LDA correlation functional - case ('RS51') call RS51_lda_exchange_energy(nGrid,weight,rho,Ex) -! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 exchange functional - case ('RMFL20') call RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) + case ('RGIC') + + call RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + case default call print_warning('!!! LDA exchange functional not available !!!') diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 index 0d120f1..1a90522 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -32,6 +32,10 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weigh call RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) + case ('RGIC') + + call RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) + case default call print_warning('!!! LDA correlation functional not available !!!') diff --git a/src/eDFT/lda_exchange_potential.f90 b/src/eDFT/lda_exchange_potential.f90 index d3e5bc9..0d56ed8 100644 --- a/src/eDFT/lda_exchange_potential.f90 +++ b/src/eDFT/lda_exchange_potential.f90 @@ -26,24 +26,22 @@ subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,A select case (DFA) -! Restricted version of Slater's LDA correlation functional - case ('RS51') call RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) -! Slater's LDA correlation functional - case ('S51') call S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) -! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 functional - case ('RMFL20') call RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + case ('RGIC') + + call RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + case default call print_warning('!!! LDA exchange functional not available !!!')