From 59d798afc9ea30feec55e21aacb15106680788cd Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 27 Mar 2020 20:46:13 +0100 Subject: [PATCH] fix bug in print energy eDFT --- input/basis | 69 +++++-------------- input/dft | 4 +- input/weight | 69 +++++-------------- src/QuAcK/excitation_density_Tmatrix.f90 | 24 +++---- ...a_correlation_derivative_discontinuity.f90 | 2 +- src/eDFT/RMFL20_lda_correlation_energy.f90 | 2 +- .../print_restricted_individual_energy.f90 | 9 +-- .../restricted_elda_correlation_energy.f90 | 3 +- src/eDFT/restricted_individual_energy.f90 | 10 +-- 9 files changed, 61 insertions(+), 131 deletions(-) diff --git a/input/basis b/input/basis index 1ce59ba..a1f7a8d 100644 --- a/input/basis +++ b/input/basis @@ -1,62 +1,27 @@ -1 14 +1 5 S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.7977000 1.0000000 + 1 0.1220000 1.0000000 S 1 - 1 0.2581000 1.0000000 -S 1 - 1 0.0898900 1.0000000 -S 1 - 1 0.0236300 1.0000000 + 1 0.0297400 1.0000000 P 1 - 1 2.2920000 1.0000000 + 1 0.7270000 1.0000000 P 1 - 1 0.8380000 1.0000000 -P 1 - 1 0.2920000 1.0000000 -P 1 - 1 0.0848000 1.0000000 -D 1 - 1 2.0620000 1.0000000 -D 1 - 1 0.6620000 1.0000000 -D 1 - 1 0.1900000 1.0000000 -F 1 - 1 1.3970000 1.0000000 -F 1 - 1 0.3600000 1.0000000 -2 14 + 1 0.1410000 1.0000000 +2 5 S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.7977000 1.0000000 + 1 0.1220000 1.0000000 S 1 - 1 0.2581000 1.0000000 -S 1 - 1 0.0898900 1.0000000 -S 1 - 1 0.0236300 1.0000000 + 1 0.0297400 1.0000000 P 1 - 1 2.2920000 1.0000000 + 1 0.7270000 1.0000000 P 1 - 1 0.8380000 1.0000000 -P 1 - 1 0.2920000 1.0000000 -P 1 - 1 0.0848000 1.0000000 -D 1 - 1 2.0620000 1.0000000 -D 1 - 1 0.6620000 1.0000000 -D 1 - 1 0.1900000 1.0000000 -F 1 - 1 1.3970000 1.0000000 -F 1 - 1 0.3600000 1.0000000 + 1 0.1410000 1.0000000 + diff --git a/input/dft b/input/dft index 7d25eb7..1e4c470 100644 --- a/input/dft +++ b/input/dft @@ -6,14 +6,14 @@ # GGA = 2: # Hybrid = 4 # Hartree-Fock = 666 - 1 RMFL20 + 666 HF # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 - 1 RMFL20 + 666 HF # quadrature grid SG-n 1 # Number of states in ensemble (nEns) diff --git a/input/weight b/input/weight index 1ce59ba..a1f7a8d 100644 --- a/input/weight +++ b/input/weight @@ -1,62 +1,27 @@ -1 14 +1 5 S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.7977000 1.0000000 + 1 0.1220000 1.0000000 S 1 - 1 0.2581000 1.0000000 -S 1 - 1 0.0898900 1.0000000 -S 1 - 1 0.0236300 1.0000000 + 1 0.0297400 1.0000000 P 1 - 1 2.2920000 1.0000000 + 1 0.7270000 1.0000000 P 1 - 1 0.8380000 1.0000000 -P 1 - 1 0.2920000 1.0000000 -P 1 - 1 0.0848000 1.0000000 -D 1 - 1 2.0620000 1.0000000 -D 1 - 1 0.6620000 1.0000000 -D 1 - 1 0.1900000 1.0000000 -F 1 - 1 1.3970000 1.0000000 -F 1 - 1 0.3600000 1.0000000 -2 14 + 1 0.1410000 1.0000000 +2 5 S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.7977000 1.0000000 + 1 0.1220000 1.0000000 S 1 - 1 0.2581000 1.0000000 -S 1 - 1 0.0898900 1.0000000 -S 1 - 1 0.0236300 1.0000000 + 1 0.0297400 1.0000000 P 1 - 1 2.2920000 1.0000000 + 1 0.7270000 1.0000000 P 1 - 1 0.8380000 1.0000000 -P 1 - 1 0.2920000 1.0000000 -P 1 - 1 0.0848000 1.0000000 -D 1 - 1 2.0620000 1.0000000 -D 1 - 1 0.6620000 1.0000000 -D 1 - 1 0.1900000 1.0000000 -F 1 - 1 1.3970000 1.0000000 -F 1 - 1 0.3600000 1.0000000 + 1 0.1410000 1.0000000 + diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index b8daed0..f0a8aa3 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -73,10 +73,10 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do c=nO+1,nBas-nR do d=c,nBas-nR cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,c,d)*X2(cd,ij) -! rho2(p,a,ij) = rho2(p,a,ij) & -! + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) & -! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d))) +! rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,c,d)*X2(cd,ij) + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) & + /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -84,10 +84,10 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do k=nC+1,nO do l=k,nO kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,k,l)*Y2(kl,ij) -! rho2(p,a,ij) = rho2(p,a,ij) & -! + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) & -! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l))) +! rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,k,l)*Y2(kl,ij) + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) & + /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -113,7 +113,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do c=nO+1,nBas-nR do d=c+1,nBas-nR cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) + (2d0+1d0/sqrt(2d0))*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) + rho1(p,i,ab) = rho1(p,i,ab) + 1.5d0*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) end do end do @@ -121,7 +121,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do k=nC+1,nO do l=k+1,nO kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) + (2d0+1d0/sqrt(2d0))*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) + rho1(p,i,ab) = rho1(p,i,ab) + 1.5d0*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) end do end do @@ -135,7 +135,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r 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) + rho2(p,a,ij) = rho2(p,a,ij) + 0.5d0*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) end do end do @@ -143,7 +143,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r 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) + rho2(p,a,ij) = rho2(p,a,ij) + 0.5d0*(ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) end do end do diff --git a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 index eb313aa..d8eab4c 100644 --- a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 @@ -42,7 +42,7 @@ subroutine RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weigh do iEns=1,nEns - call restricted_elda_correlation_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rhow(:),dEcdw(iEns)) + call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rhow(:),dEcdw(iEns)) end do diff --git a/src/eDFT/RMFL20_lda_correlation_energy.f90 b/src/eDFT/RMFL20_lda_correlation_energy.f90 index 5ebe526..8819239 100644 --- a/src/eDFT/RMFL20_lda_correlation_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_energy.f90 @@ -43,7 +43,7 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) do iEns=1,nEns - call restricted_elda_correlation_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rho(:),EceLDA(iEns)) + call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rho(:),EceLDA(iEns)) end do diff --git a/src/eDFT/print_restricted_individual_energy.f90 b/src/eDFT/print_restricted_individual_energy.f90 index c560b18..847aa17 100644 --- a/src/eDFT/print_restricted_individual_energy.f90 +++ b/src/eDFT/print_restricted_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine print_restricted_individual_energy(nEns,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E, & +subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E, & Om,Omx,Omc,Omxc,OmxDD,OmcDD,OmxcDD) ! Print individual energies for eDFT calculation @@ -9,6 +9,7 @@ subroutine print_restricted_individual_energy(nEns,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,E ! Input variables integer,intent(in) :: nEns + double precision,intent(in) :: ENuc double precision,intent(in) :: Ew double precision,intent(in) :: EwGIC double precision,intent(in) :: ET(nEns) @@ -32,8 +33,8 @@ subroutine print_restricted_individual_energy(nEns,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,E write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') ' ENSEMBLE ENERGIES' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew, ' au' - write(*,'(A44,F16.10,A3)') ' GIC Ensemble energy: ',EwGIC,' au' + write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew + ENuc,' au' + write(*,'(A44,F16.10,A3)') ' GIC Ensemble energy: ',EwGIC + ENuc,' au' write(*,'(A60)') '-------------------------------------------------' write(*,*) @@ -126,7 +127,7 @@ subroutine print_restricted_individual_energy(nEns,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,E write(*,'(A60)') ' INDIVIDUAL AND EXCITATION ENERGIES' write(*,'(A60)') '-------------------------------------------------' do iEns=1,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns),' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' end do write(*,'(A60)') '-------------------------------------------------' diff --git a/src/eDFT/restricted_elda_correlation_energy.f90 b/src/eDFT/restricted_elda_correlation_energy.f90 index a54d88f..06bda0a 100644 --- a/src/eDFT/restricted_elda_correlation_energy.f90 +++ b/src/eDFT/restricted_elda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine restricted_elda_correlation_energy(nEns,aMFL,nGrid,weight,rho,Ec) +subroutine restricted_elda_correlation_energy(aMFL,nGrid,weight,rho,Ec) ! Compute the restricted LDA correlation energy of 2-glomium for various states @@ -7,7 +7,6 @@ subroutine restricted_elda_correlation_energy(nEns,aMFL,nGrid,weight,rho,Ec) ! Input variables - integer,intent(in) :: nEns double precision,intent(in) :: aMFL(3) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index 69d26cb..0f97e30 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -121,9 +121,9 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri !------------------------------------------------------------------------ do iEns=1,nEns - Exc(iEns) = Ex(iEns) + Ec(iEns) - E(iEns) = ENuc + ET(iEns) + EV(iEns) + EJ(iEns) & - + Ex(iEns) + Ec(iEns) + ExcDD(iEns) + Exc(iEns) = Ex(iEns) + Ec(iEns) + E(iEns) = ET(iEns) + EV(iEns) + EJ(iEns) & + + Ex(iEns) + Ec(iEns) + ExcDD(iEns) end do !------------------------------------------------------------------------ @@ -157,8 +157,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri ! Dump results !------------------------------------------------------------------------ - call print_restricted_individual_energy(nEns,Ew,EwGIC,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & - ExDD(:),EcDD(:),ExcDD(:),E(:), & + call print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & + ExDD(:),EcDD(:),ExcDD(:),E(:), & Om(:),Omx(:),Omc(:),Omxc(:),OmxDD(:),OmcDD(:),OmxcDD(:)) end subroutine restricted_individual_energy