diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 deleted file mode 100644 index 263fb51..0000000 --- a/src/eDFT/LIM_RKS.f90 +++ /dev/null @@ -1,150 +0,0 @@ -subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGrid,weight, & - maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,& - c,occnum,Cx_choice) - -! Perform restricted Kohn-Sham calculation for ensembles - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: x_rung,c_rung - character(len=12),intent(in) :: x_DFA,c_DFA - logical,intent(in) :: LDA_centered - integer,intent(in) :: nEns - double precision,intent(in) :: aCC_w1(3) - double precision,intent(in) :: aCC_w2(3) - integer,intent(in) :: nGrid - double precision,intent(in) :: weight(nGrid) - integer,intent(in) :: maxSCF,max_diis,guess_type - double precision,intent(in) :: thresh - integer,intent(in) :: nBas - double precision,intent(in) :: AO(nBas,nGrid) - double precision,intent(in) :: dAO(ncart,nBas,nGrid) - - integer,intent(in) :: nO,nV - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ENuc - double precision,intent(in) :: occnum(nBas,nspin,nEns) - integer,intent(in) :: Cx_choice - - - double precision,intent(out) :: c(nBas,nBas) - -! Local variables - - integer :: iEns - double precision :: Ew(nEnS) - double precision :: wLIM(nEns) - double precision :: Om(nEns) - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'* Linear-interpolation method *' - write(*,*)'* for excitation energies *' - write(*,*)'************************************************' - write(*,*) - -! Initializatio - - Ew(:) = 0d0 - Om(:) = 0d0 - -!------------------------------------------------------------------------ -! Zero-weight calculation -!------------------------------------------------------------------------ - - write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' ZERO-WEIGHT CALCULATION ' - write(*,'(A40)') '*************************************************' - - wLIM(1) = 1d0 - wLIM(2) = 0d0 - wLIM(3) = 0d0 - - do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns) - end do - write(*,'(A40)') '*************************************************' - write(*,*) - - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum,Cx_choice) - -!------------------------------------------------------------------------ -! Equiensemble calculation -!------------------------------------------------------------------------ - - write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' TWO-STATE EQUI-WEIGHT CALCULATION ' - write(*,'(A40)') '*************************************************' - - wLIM(1) = 0.5d0 - wLIM(2) = 0.0d0 - wLIM(3) = 0.5d0 - - do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns) - end do - write(*,'(A40)') '*************************************************' - write(*,*) - - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum,Cx_choice) - -!------------------------------------------------------------------------ -! Equiensemble calculation -!------------------------------------------------------------------------ - - write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' THREE-STATE EQUI-WEIGHT CALCULATION ' - write(*,'(A40)') '*************************************************' - - wLIM(1) = 1d0/3d0 - wLIM(2) = 1d0/3d0 - wLIM(3) = 1d0/3d0 - - do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns) - end do - write(*,'(A40)') '*************************************************' - write(*,*) - -! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & -! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum,Cx_choice) - - -!------------------------------------------------------------------------ -! LIM excitation energies -!------------------------------------------------------------------------ - - Om(2) = 2d0*(Ew(2) - Ew(1)) -! Om(3) = 3d0*(Ew(3) - Ew(2)) + 0.5d0*Om(2) - - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES ' - write(*,'(A60)') '-------------------------------------------------' - - write(*,'(A44,F16.10,A3)') ' Ensemble energy #1 ',Ew(1),' au' - write(*,'(A44,F16.10,A3)') ' Ensemble energy #2 ',Ew(2),' au' - write(*,'(A44,F16.10,A3)') ' Ensemble energy #3 ',Ew(3),' au' - write(*,'(A60)') '-------------------------------------------------' - do iEns=2,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au' - end do - write(*,*) - do iEns=2,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV' - end do - write(*,'(A60)') '-------------------------------------------------' - - -end subroutine LIM_RKS diff --git a/src/eDFT/MOM_RKS.f90 b/src/eDFT/MOM_RKS.f90 deleted file mode 100644 index 35e2c8c..0000000 --- a/src/eDFT/MOM_RKS.f90 +++ /dev/null @@ -1,148 +0,0 @@ -subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & - aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,& - nO,nV,S,T,V,Hc,ERI,X,ENuc,c,occnum,Cx_choice) - -! Perform restricted Kohn-Sham calculation for ensembles - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: x_rung,c_rung - character(len=12),intent(in) :: x_DFA,c_DFA - logical,intent(in) :: LDA_centered - integer,intent(in) :: nEns - integer,intent(in) :: nGrid - double precision,intent(in) :: weight(nGrid) - double precision,intent(in) :: aCC_w1(3) - double precision,intent(in) :: aCC_w2(3) - integer,intent(in) :: maxSCF,max_diis,guess_type - double precision,intent(in) :: thresh - integer,intent(in) :: nBas - double precision,intent(in) :: AO(nBas,nGrid) - double precision,intent(in) :: dAO(ncart,nBas,nGrid) - - integer,intent(in) :: nO,nV - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ENuc - double precision,intent(in) :: occnum(nBas,nspin,nEns) - integer,intent(in) :: Cx_choice - - - double precision,intent(out) :: c(nBas,nBas) - -! Local variables - - integer :: iEns - double precision :: Ew(nEns) - double precision :: wMOM(nEns) - double precision :: Om(nEns) - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'* Maximum Overlao method *' - write(*,*)'* for excitation energies *' - write(*,*)'************************************************' - write(*,*) - -! Initialization - - Ew(:) = 0d0 - Om(:) = 0d0 - -!------------------------------------------------------------------------ -! Zero-weight calculation -!------------------------------------------------------------------------ - - write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' STATE-SPECIFIC CALCULATION #1 ' - write(*,'(A40)') '*************************************************' - - wMOM(1) = 1d0 - wMOM(2) = 0d0 - wMOM(3) = 0d0 - - do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wMOM(iEns) - end do - write(*,'(A40)') '*************************************************' - write(*,*) - - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c,occnum,Cx_choice) - -!------------------------------------------------------------------------ -! Equiensemble calculation -!------------------------------------------------------------------------ - - write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' STATE-SPECIFIC CALCULATION #2 ' - write(*,'(A40)') '*************************************************' - - wMOM(1) = 0d0 - wMOM(2) = 1d0 - wMOM(3) = 0d0 - - do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wMOM(iEns) - end do - write(*,'(A40)') '*************************************************' - write(*,*) - -! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & -! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c,occnum,Cx_choice) - -!------------------------------------------------------------------------ -! Equiensemble calculation -!------------------------------------------------------------------------ - - write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' STATE-SPECIFIC CALCULATION #3 ' - write(*,'(A40)') '*************************************************' - - wMOM(1) = 0d0 - wMOM(2) = 0d0 - wMOM(3) = 1.0d0 - - do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wMOM(iEns) - end do - write(*,'(A40)') '*************************************************' - write(*,*) - - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c,occnum,Cx_choice) - -!------------------------------------------------------------------------ -! MOM excitation energies -!------------------------------------------------------------------------ - - Om(:) = Ew(:) - Ew(1) - - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' MAXIMUM OVERLAP METHOD EXCITATION ENERGIES ' - write(*,'(A60)') '-------------------------------------------------' - - write(*,'(A44,F16.10,A3)') ' Ensemble energy #1 ',Ew(1),' au' - write(*,'(A44,F16.10,A3)') ' Ensemble energy #2 ',Ew(2),' au' - write(*,'(A44,F16.10,A3)') ' Ensemble energy #3 ',Ew(3),' au' - write(*,'(A60)') '-------------------------------------------------' - do iEns=2,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au' - end do - write(*,*) - do iEns=2,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV' - end do - write(*,'(A60)') '-------------------------------------------------' - - -end subroutine MOM_RKS diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 04dcf32..0ce947b 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -147,80 +147,6 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid)) call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO) - LDA_centered = .true. - -!------------------------------------------------------------------------ -! Compute GOK-RKS energy -!------------------------------------------------------------------------ - -! if(method == 'GOK-RKS') then - -! call cpu_time(start_KS) -! call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, & -! maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO(1),nV(1), & -! S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice) -! call cpu_time(end_KS) - -! t_KS = end_KS - start_KS -! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOK-RKS = ',t_KS,' seconds' -! write(*,*) - -!end if - -!------------------------------------------------------------------------ -! Compute LIM excitation energies -!------------------------------------------------------------------------ - -! if(method == 'LIM-RKS') then - -! call cpu_time(start_KS) -! call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & -! aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & -! S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered) -! call cpu_time(end_KS) - -! t_KS = end_KS - start_KS -! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for LIM-RKS = ',t_KS,' seconds' -! write(*,*) - -! end if - -!------------------------------------------------------------------------ -! Compute MOM excitation energies -!------------------------------------------------------------------------ - -! if(method == 'MOM-RKS') then - -! call cpu_time(start_KS) -! call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & -! aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & -! S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:),occnum,Cx_choice,doNcentered) -! call cpu_time(end_KS) - -! t_KS = end_KS - start_KS -! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds' -! write(*,*) - -! end if - -!------------------------------------------------------------------------ -! Compute GOK-UKS energy (BROKEN) -!------------------------------------------------------------------------ - -! if(method == 'GOK-UKS') then - -! call cpu_time(start_KS) -! call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, & -! nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,occnum, & -! Cx_choice,doNcentered) -! call cpu_time(end_KS) - -! t_KS = end_KS - start_KS -! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UKS = ',t_KS,' seconds' -! write(*,*) - -! end if - !------------------------------------------------------------------------ ! Compute UKS energy !------------------------------------------------------------------------ diff --git a/src/eDFT/fock_exchange_energy.f90 b/src/eDFT/fock_exchange_energy.f90 deleted file mode 100644 index 17e7d49..0000000 --- a/src/eDFT/fock_exchange_energy.f90 +++ /dev/null @@ -1,25 +0,0 @@ -subroutine fock_exchange_energy(nBas,P,Fx,Ex) - -! Compute the (exact) Fock exchange energy - - implicit none - -! Input variables - - integer,intent(in) :: nBas - double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: Fx(nBas,nBas) - -! Local variables - - double precision,external :: trace_matrix - -! Output variables - - double precision,intent(out) :: Ex - -! Compute HF exchange energy - - Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx)) - -end subroutine fock_exchange_energy diff --git a/src/eDFT/print_RKS.f90 b/src/eDFT/print_RKS.f90 deleted file mode 100644 index 509559d..0000000 --- a/src/eDFT/print_RKS.f90 +++ /dev/null @@ -1,75 +0,0 @@ -subroutine print_RKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) - -! Print one- and two-electron energies and other stuff for KS calculation - - implicit none - include 'parameters.h' - - integer,intent(in) :: nBas - integer,intent(in) :: nO - double precision,intent(in) :: eps(nBas) - double precision,intent(in) :: c(nBas,nBas) - double precision,intent(in) :: ENuc - double precision,intent(in) :: ET - double precision,intent(in) :: EV - double precision,intent(in) :: EJ - double precision,intent(in) :: Ex - double precision,intent(in) :: Ec - double precision,intent(in) :: Ew - - integer :: HOMO - integer :: LUMO - double precision :: Gap - -! HOMO, LUMO, and Gap - - HOMO = nO - - LUMO = HOMO + 1 - - Gap = eps(LUMO) - eps(HOMO) - -! Dump results - - - write(*,*) - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40)') ' Summary ' - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au' - write(*,'(A40,1X,F16.10,A3)') ' Kinetic energy: ',ET,' au' - write(*,'(A40,1X,F16.10,A3)') ' Potential energy: ',EV,' au' - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex + Ec,' au' - write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',EJ,' au' - write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au' - write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',Ec,' au' - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',Ew ,' au' - write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ', ENuc,' au' - write(*,'(A40,1X,F16.10,A3)') ' Kohn-Sham energy: ',Ew + ENuc,' au' - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F13.6,A3)') ' KS HOMO energy:',eps(HOMO),' au' - write(*,'(A40,F13.6,A3)') ' KS LUMO energy:',eps(LUMO),' au' - write(*,'(A40,F13.6,A3)') ' KS HOMO-LUMO gap:',Gap ,' au' - write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F13.6,A3)') ' KS HOMO energy:',eps(HOMO)*HatoeV,' eV' - write(*,'(A40,F13.6,A3)') ' KS LUMO energy:',eps(LUMO)*HatoeV,' eV' - write(*,'(A40,F13.6,A3)') ' KS HOMO-LUMO gap:',Gap *HatoeV,' eV' - write(*,'(A60)') '-------------------------------------------------' - write(*,*) - -! Print results - - write(*,'(A50)') '-----------------------------------------' - write(*,'(A50)') ' Kohn-Sham orbital coefficients ' - write(*,'(A50)') '-----------------------------------------' - call matout(nBas,nBas,c(:,:)) - write(*,*) - write(*,'(A50)') '---------------------------------------' - write(*,'(A50)') ' Kohn-Sham orbital energies ' - write(*,'(A50)') '---------------------------------------' - call matout(nBas,1,eps(:)) - write(*,*) - -end subroutine print_RKS