PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-04 13:13:51 +01:00
Code Releases Activity

remove GIC print

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-04-09 21:58:10 +02:00
parent 4587d3ff5c
commit ca2186bc6f
12 changed files with 104 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
input/basis
View File

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

4
input/dft
View File

@ -6,7 +6,7 @@
# GGA = 2: RB88
# Hybrid = 4
# Hartree-Fock = 666
1 RGIC
666 HF
# correlation rung:
# Hartree = 0
# LDA = 1: RVWN5,RMFL20
@ -19,6 +19,6 @@
# Number of states in ensemble (nEns)
2
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.5
1.0
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type
32 0.00001 T 5 1 1

48
input/weight
View File

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

5
src/eDFT/GOK_RKS.f90
View File

@ -1,5 +1,5 @@
subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,EwGIC,c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c)
! Perform restricted Kohn-Sham calculation for ensembles
@ -80,7 +80,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri
! Output variables
double precision,intent(out) :: Ew
double precision,intent(out) :: EwGIC
! Hello world
@ -343,6 +342,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri
call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), &
nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), &
J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:))
J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:))
end subroutine GOK_RKS

27
src/eDFT/LIM_RKS.f90
View File

@ -35,10 +35,10 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
! Local variables
integer :: iEns
double precision :: EwZW,EwGICZW
double precision :: EwEW,EwGICEW
double precision :: EwZW
double precision :: EwEW
double precision :: wLIM(nEns)
double precision :: Om(nEns),OmGIC(nEns)
double precision :: Om(nEns)
! Hello world
@ -67,7 +67,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
write(*,*)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,EwGICZW,c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,c)
!------------------------------------------------------------------------
! Equiensemble calculation
@ -84,20 +84,15 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
write(*,*)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,EwGICEW,c)
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,c)
!------------------------------------------------------------------------
! LIM excitation energies
!------------------------------------------------------------------------
Om(1) = 0d0
OmGIC(1) = 0d0
Om(2) = 0d0
OmGIC(2) = 0d0
Om(:) = 0d0
if(wEns(2) > 10d-3) then
Om(2) = (EwEW - EwZW)/wEns(2)
OmGIC(2) = (EwGICEW - EwGICZW)/wEns(2)
Om(2) = (EwEW - EwZW)/wEns(2)
end if
write(*,'(A60)') '-------------------------------------------------'
@ -105,20 +100,12 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A44,F16.10,A3)') ' Zero-weight ensemble energy',EwZW, ' au'
write(*,'(A44,F16.10,A3)') ' Zero-weight GIC ensemble energy',EwGICZW,' au'
write(*,*)
write(*,'(A44,F16.10,A3)') ' Equi-weight ensemble energy',EwEW, ' au'
write(*,'(A44,F16.10,A3)') ' Equi-weight GIC ensemble energy',EwGICEW, ' au'
write(*,'(A60)') '-------------------------------------------------'
do iEns=2,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au'
write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV'
end do
write(*,*)
do iEns=2,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' GIC Excitation energy 1 ->',iEns,': ',OmGIC(iEns), ' au'
write(*,'(A40,I2,A2,F16.10,A3)') ' GIC Excitation energy 1 ->',iEns,': ',OmGIC(iEns)*HaToeV,' eV'
end do
write(*,'(A60)') '-------------------------------------------------'

12
src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90
View File

@ -43,15 +43,15 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho
! Parameters for stretch H2
a = + 0.01922622507087411d0
b = - 0.01799647558018601d0
c = - 0.022945430666782573d0
! a = + 0.01922622507087411d0
! b = - 0.01799647558018601d0
! c = - 0.022945430666782573d0
! Parameters for He
! a = 1.9125735895875828d0
! b = 2.715266992840757d0
! c = 2.1634223380633086d0
a = 1.9125735895875828d0
b = 2.715266992840757d0
c = 2.1634223380633086d0
w = wEns(2)
dCxGICdw = (0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0)))

12
src/eDFT/RGIC_lda_exchange_energy.f90
View File

@ -35,15 +35,15 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Parameters for stretch H2
a = + 0.01922622507087411d0
b = - 0.01799647558018601d0
c = - 0.022945430666782573d0
! a = + 0.01922622507087411d0
! b = - 0.01799647558018601d0
! c = - 0.022945430666782573d0
! Parameters for He
! a = 1.9125735895875828d0
! b = 2.715266992840757d0
! c = 2.1634223380633086d0
a = 1.9125735895875828d0
b = 2.715266992840757d0
c = 2.1634223380633086d0
w = wEns(2)
CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)

12
src/eDFT/RGIC_lda_exchange_individual_energy.f90
View File

@ -38,15 +38,15 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E
! Parameters for stretch H2
a = + 0.01922622507087411d0
b = - 0.01799647558018601d0
c = - 0.022945430666782573d0
! a = + 0.01922622507087411d0
! b = - 0.01799647558018601d0
! c = - 0.022945430666782573d0
! Parameters for He
! a = 1.9125735895875828d0
! b = 2.715266992840757d0
! c = 2.1634223380633086d0
a = 1.9125735895875828d0
b = 2.715266992840757d0
c = 2.1634223380633086d0
w = wEns(2)
CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)

12
src/eDFT/RGIC_lda_exchange_potential.f90
View File

@ -38,15 +38,15 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! Parameters for stretch H2
a = + 0.01922622507087411d0
b = - 0.01799647558018601d0
c = - 0.022945430666782573d0
! a = + 0.01922622507087411d0
! b = - 0.01799647558018601d0
! c = - 0.022945430666782573d0
! Parameters for He
! a = 1.9125735895875828d0
! b = 2.715266992840757d0
! c = 2.1634223380633086d0
a = 1.9125735895875828d0
b = 2.715266992840757d0
c = 2.1634223380633086d0
w = wEns(2)
CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)

4
src/eDFT/eDFT.f90
View File

@ -7,7 +7,7 @@ program eDFT
integer :: nNuc,nBas
integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin)
double precision :: ENuc,Ew,EwGIC
double precision :: ENuc,Ew
double precision,allocatable :: ZNuc(:),rNuc(:,:)
@ -164,7 +164,7 @@ program eDFT
call cpu_time(start_KS)
call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), &
maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), &
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,EwGIC,c(:,:))
S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,c(:,:))
call cpu_time(end_KS)
t_KS = end_KS - start_KS

4
src/eDFT/print_restricted_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, &
subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, &
Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD)
! Print individual energies for eDFT calculation
@ -11,7 +11,6 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,
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)
double precision,intent(in) :: EV(nEns)
double precision,intent(in) :: EJ(nEns)
@ -36,7 +35,6 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,
write(*,'(A60)') ' ENSEMBLE ENERGIES'
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew + ENuc,' au'
write(*,'(A44,F16.10,A3)') ' GIC Ensemble energy: ',EwGIC + ENuc,' au'
write(*,'(A60)') '-------------------------------------------------'
write(*,*)

14
src/eDFT/restricted_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,EwGIC,E,Om)
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om)
! Compute individual energies as well as excitation energies
@ -55,7 +55,6 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
! Output variables
double precision,intent(out) :: EwGIC
double precision,intent(out) :: E(nEns)
double precision,intent(out) :: Om(nEns)
@ -129,15 +128,6 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
+ Ex(iEns) + Ec(iEns) + ExcDD(iEns)
end do
!------------------------------------------------------------------------
! Total energy with ghost-interaction correction
!------------------------------------------------------------------------
EwGIC = 0d0
do iEns=1,nEns
EwGIC = EwGIC + wEns(iEns)*E(iEns)
end do
!------------------------------------------------------------------------
! Excitation energies
!------------------------------------------------------------------------
@ -162,7 +152,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n
! Dump results
!------------------------------------------------------------------------
call print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), &
call print_restricted_individual_energy(nEns,ENuc,Ew,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), &
Eaux(:),ExDD(:),EcDD(:),ExcDD(:),E(:), &
Om(:),Omx(:),Omc(:),Omxc(:),Omaux,OmxDD(:),OmcDD(:),OmxcDD(:))