4
1
mirror of https://github.com/pfloos/quack synced 2024-06-02 11:25:28 +02:00

more clean up

This commit is contained in:
Pierre-Francois Loos 2021-11-30 09:35:04 +01:00
parent c711a7e216
commit f00215d9ad
5 changed files with 0 additions and 472 deletions

View File

@ -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

View File

@ -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

View File

@ -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
!------------------------------------------------------------------------

View File

@ -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

View File

@ -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