From 431eab634da0c98b82ca3a320393d932f343b7d2 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 5 May 2020 16:45:10 +0200 Subject: [PATCH] saving eDFT --- input/basis | 69 +++++++++++---- input/dft | 4 +- input/methods | 4 +- input/options | 4 +- input/weight | 69 +++++++++++---- src/eDFT/GOK_RKS.f90 | 2 +- src/eDFT/MOM_RKS.f90 | 6 +- ..._lda_exchange_derivative_discontinuity.f90 | 85 ------------------- src/eDFT/RGIC_lda_exchange_energy.f90 | 68 --------------- .../RGIC_lda_exchange_individual_energy.f90 | 72 ---------------- src/eDFT/RGIC_lda_exchange_potential.f90 | 75 ---------------- .../lda_exchange_derivative_discontinuity.f90 | 4 +- src/eDFT/lda_exchange_energy.f90 | 4 +- src/eDFT/lda_exchange_individual_energy.f90 | 4 +- src/eDFT/lda_exchange_potential.f90 | 4 +- src/eDFT/restricted_auxiliary_energy.f90 | 2 +- src/eDFT/restricted_density_matrix.f90 | 2 +- 17 files changed, 124 insertions(+), 354 deletions(-) delete mode 100644 src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 delete mode 100644 src/eDFT/RGIC_lda_exchange_energy.f90 delete mode 100644 src/eDFT/RGIC_lda_exchange_individual_energy.f90 delete mode 100644 src/eDFT/RGIC_lda_exchange_potential.f90 diff --git a/input/basis b/input/basis index a1f7a8d..1ce59ba 100644 --- a/input/basis +++ b/input/basis @@ -1,27 +1,62 @@ -1 5 +1 14 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 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 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 - + 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 diff --git a/input/dft b/input/dft index 79fc2fe..3bb1b92 100644 --- a/input/dft +++ b/input/dft @@ -1,12 +1,12 @@ # Restricted or unrestricted KS calculation - MOM-RKS + LIM-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 # GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 1 RGIC + 1 RCC # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 diff --git a/input/methods b/input/methods index 8afbf57..582b287 100644 --- a/input/methods +++ b/input/methods @@ -7,7 +7,7 @@ # drCCD rCCD lCCD pCCD F F F F # CIS CID CISD - T F F + F F F # RPA RPAx ppRPA F F F # G0F2 evGF2 G0F3 evGF3 @@ -15,6 +15,6 @@ # G0W0 evGW qsGW T F F # G0T0 evGT qsGT - F F F + T F F # MCMP2 F diff --git a/input/options b/input/options index 51ab7a5..417c173 100644 --- a/input/options +++ b/input/options @@ -8,8 +8,8 @@ T T # GF: maxSCF thresh DIIS n_diis lin renorm 256 0.00001 T 5 T 3 -# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta - 256 0.00001 T 5 F F T F F F T 0.00367493 +# GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta + 256 0.00001 T 5 F F T F F F T 0.0 # ACFDT: AC Kx XBS F F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/input/weight b/input/weight index a1f7a8d..1ce59ba 100644 --- a/input/weight +++ b/input/weight @@ -1,27 +1,62 @@ -1 5 +1 14 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 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 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 - + 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 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index 7df8b16..b8814b8 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -328,7 +328,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - stop +! stop end if diff --git a/src/eDFT/MOM_RKS.f90 b/src/eDFT/MOM_RKS.f90 index 5f41a01..24b716d 100644 --- a/src/eDFT/MOM_RKS.f90 +++ b/src/eDFT/MOM_RKS.f90 @@ -91,8 +91,8 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c) +! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & +! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c) !------------------------------------------------------------------------ ! Equiensemble calculation @@ -122,7 +122,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & Om(:) = Ew(:) - Ew(1) write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES ' + write(*,'(A60)') ' MAXIMUM OVERLAP METHOD EXCITATION ENERGIES ' write(*,'(A60)') '-------------------------------------------------' write(*,'(A44,F16.10,A3)') ' Ensemble energy #1 ',Ew(1),' au' diff --git a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 deleted file mode 100644 index 9437c10..0000000 --- a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 +++ /dev/null @@ -1,85 +0,0 @@ -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,allocatable :: dExdw(:) - double precision,external :: Kronecker_delta - - double precision :: a,b,c,w - double precision :: dCxGICdw - -! Output variables - - double precision,intent(out) :: ExDD(nEns) - -! Memory allocation - - allocate(dExdw(nEns)) - -! Weight-dependent Cx coefficient for RMFL20 exchange functional - -! Parameters for H2 at equilibrium - - a = + 0.5739189000851961d0 - b = - 0.0003469882157336496d0 - c = - 0.2676338054343272d0 - -! Parameters for stretch H2 - -! a = + 0.01918229168254928d0 -! b = - 0.01545313842512261d0 -! c = - 0.012720073519142448d0 - -! Parameters for He - -! a = 1.9015719148496788d0 -! b = 2.5236598782764412d0 -! c = 1.6652282199359842d0 - - w = 0.5d0*wEns(2) + wEns(3) - dCxGICdw = (0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0))) - dCxGICdw = CxLDA*dCxGICdw - - dExdw(:) = 0d0 - - do iG=1,nGrid - - r = max(0d0,rhow(iG)) - - if(r > threshold) then - - dExdw(1) = 0d0 - dExdw(2) = dExdw(2) + 0.5d0*weight(iG)*dCxGICdw*r**(4d0/3d0) - dExdw(3) = dExdw(3) + 1.0d0*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 deleted file mode 100644 index d8a5844..0000000 --- a/src/eDFT/RGIC_lda_exchange_energy.f90 +++ /dev/null @@ -1,68 +0,0 @@ -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-dependent Cx coefficient - -! Weight-dependent Cx coefficient for RMFL20 exchange functional - -! Parameters for H2 at equilibrium - - a = + 0.5739189000851961d0 - b = - 0.0003469882157336496d0 - c = - 0.2676338054343272d0 - -! Parameters for stretch H2 - -! a = + 0.01918229168254928d0 -! b = - 0.01545313842512261d0 -! c = - 0.012720073519142448d0 - -! Parameters for He - -! a = 1.9015719148496788d0 -! b = 2.5236598782764412d0 -! c = 1.6652282199359842d0 - - w = 0.5d0*wEns(2) + wEns(3) - CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) - CxGIC = CxLDA*CxGIC - -! Compute GIC-LDA exchange energy - - Ex = 0d0 - - do iG=1,nGrid - - r = max(0d0,rho(iG)) - - if(r > threshold) then - 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 deleted file mode 100644 index 42a53f4..0000000 --- a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 +++ /dev/null @@ -1,72 +0,0 @@ -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 - -! Parameters for H2 at equilibrium - - a = + 0.5739189000851961d0 - b = - 0.0003469882157336496d0 - c = - 0.2676338054343272d0 - -! Parameters for stretch H2 - -! a = + 0.01918229168254928d0 -! b = - 0.01545313842512261d0 -! c = - 0.012720073519142448d0 - -! Parameters for He - -! a = 1.9015719148496788d0 -! b = 2.5236598782764412d0 -! c = 1.6652282199359842d0 - - w = 0.5d0*wEns(2) + wEns(3) - CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) - CxGIC = CxLDA*CxGIC - -! 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 = 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 deleted file mode 100644 index f8f0426..0000000 --- a/src/eDFT/RGIC_lda_exchange_potential.f90 +++ /dev/null @@ -1,75 +0,0 @@ -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 - -! Parameters for H2 at equilibrium - - a = + 0.5739189000851961d0 - b = - 0.0003469882157336496d0 - c = - 0.2676338054343272d0 - -! Parameters for stretch H2 - -! a = + 0.01918229168254928d0 -! b = - 0.01545313842512261d0 -! c = - 0.012720073519142448d0 - -! Parameters for He - -! a = 1.9015719148496788d0 -! b = 2.5236598782764412d0 -! c = 1.6652282199359842d0 - - w = 0.5d0*wEns(2) + wEns(3) - CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) - CxGIC = CxLDA*CxGIC - -! 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*CxGIC*r**(1d0/3d0) - - endif - - enddo - enddo - enddo - -end subroutine RGIC_lda_exchange_potential diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 index 0f0dac6..75296f7 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -37,9 +37,9 @@ 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') + case ('RCC') - call RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) + call RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) case default diff --git a/src/eDFT/lda_exchange_energy.f90 b/src/eDFT/lda_exchange_energy.f90 index d2f1616..5b652b4 100644 --- a/src/eDFT/lda_exchange_energy.f90 +++ b/src/eDFT/lda_exchange_energy.f90 @@ -35,9 +35,9 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) call RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) - case ('RGIC') + case ('RCC') - call RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + call RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) case default diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 index 624c6eb..bf8caa6 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -32,9 +32,9 @@ 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') + case ('RCC') - call RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) + call RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) case default diff --git a/src/eDFT/lda_exchange_potential.f90 b/src/eDFT/lda_exchange_potential.f90 index 0d56ed8..38d77ab 100644 --- a/src/eDFT/lda_exchange_potential.f90 +++ b/src/eDFT/lda_exchange_potential.f90 @@ -38,9 +38,9 @@ subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,A call RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) - case ('RGIC') + case ('RCC') - call RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + call RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) case default diff --git a/src/eDFT/restricted_auxiliary_energy.f90 b/src/eDFT/restricted_auxiliary_energy.f90 index 71b1c8f..16490f5 100644 --- a/src/eDFT/restricted_auxiliary_energy.f90 +++ b/src/eDFT/restricted_auxiliary_energy.f90 @@ -37,7 +37,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux) Eaux(iEns) = 0d0 end if - Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+1) + Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+2) ! Doubly-excited state density matrix diff --git a/src/eDFT/restricted_density_matrix.f90 b/src/eDFT/restricted_density_matrix.f90 index 907cc31..93ac4b4 100644 --- a/src/eDFT/restricted_density_matrix.f90 +++ b/src/eDFT/restricted_density_matrix.f90 @@ -38,7 +38,7 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P) end if P(:,:,iEns) = P(:,:,iEns) + 1d0*matmul(c(:,nO :nO ),transpose(c(:,nO :nO ))) & - + 1d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1))) + + 1d0*matmul(c(:,nO+2:nO+2),transpose(c(:,nO+2:nO+2))) ! Doubly-excited state density matrix