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

more modifications in the final print for eDFT_UKS

This commit is contained in:
Clotilde Marut 2020-07-02 22:15:29 +02:00
parent 09d3b39e63
commit e8d4938fe5
12 changed files with 209 additions and 82 deletions

View File

@ -1,9 +1,25 @@
1 3
S 3
1 38.3600000 0.0238090
2 5.7700000 0.1548910
3 1.2400000 0.4699870
1 10
S 4
1 528.5000000 0.0009400
2 79.3100000 0.0072140
3 18.0500000 0.0359750
4 5.0850000 0.1277820
S 1
1 0.2976000 1.0000000
1 1.6090000 1.0000000
S 1
1 0.5363000 1.0000000
S 1
1 0.1833000 1.0000000
P 1
1 1.2750000 1.0000000
1 5.9940000 1.0000000
P 1
1 1.7450000 1.0000000
P 1
1 0.5600000 1.0000000
D 1
1 4.2990000 1.0000000
D 1
1 1.2230000 1.0000000
F 1
1 2.6800000 1.0000000

View File

@ -6,7 +6,7 @@
# GGA = 2: RB88
# Hybrid = 4
# Hartree-Fock = 666
1 US51
1 UCC
# correlation rung:
# Hartree = 0
# LDA = 1: RVWN5,RMFL20

View File

@ -1,9 +1,25 @@
1 3
S 3
1 38.3600000 0.0238090
2 5.7700000 0.1548910
3 1.2400000 0.4699870
1 10
S 4
1 528.5000000 0.0009400
2 79.3100000 0.0072140
3 18.0500000 0.0359750
4 5.0850000 0.1277820
S 1
1 0.2976000 1.0000000
1 1.6090000 1.0000000
S 1
1 0.5363000 1.0000000
S 1
1 0.1833000 1.0000000
P 1
1 1.2750000 1.0000000
1 5.9940000 1.0000000
P 1
1 1.7450000 1.0000000
P 1
1 0.5600000 1.0000000
D 1
1 4.2990000 1.0000000
D 1
1 1.2230000 1.0000000
F 1
1 2.6800000 1.0000000

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD)
subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD)
! Compute the restricted version of the curvature-corrected exchange ensemble derivative
@ -9,6 +9,8 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(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)
double precision,intent(in) :: rhow(nGrid)
@ -35,9 +37,9 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow
! Single excitation parameters
a1 = 0.0d0
b1 = 0.0d0
c1 = 0.0d0
! a1 = 0.0d0
! b1 = 0.0d0
! c1 = 0.0d0
! Parameters for H2 at equilibrium
@ -47,9 +49,9 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow
! Parameters for stretch H2
a2 = + 0.01922622507087411d0
b2 = - 0.01799647558018601d0
c2 = - 0.022945430666782573d0
! a2 = + 0.01922622507087411d0
! b2 = - 0.01799647558018601d0
! c2 = - 0.022945430666782573d0
! Parameters for He
@ -57,6 +59,14 @@ subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow
! b2 = 2.715266992840757d0
! c2 = 2.1634223380633086d0
a1 = aCC_w1(1)
b1 = aCC_w1(2)
c1 = aCC_w1(3)
a2 = aCC_w2(1)
b2 = aCC_w2(2)
c2 = aCC_w2(3)
w1 = wEns(2)
w2 = wEns(3)

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex)
! Compute the restricted version of the curvature-corrected exchange functional
@ -9,6 +9,8 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(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)
double precision,intent(in) :: rho(nGrid)
@ -28,9 +30,9 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Single excitation parameter
a1 = 0.0d0
b1 = 0.0d0
c1 = 0.0d0
! a1 = 0.0d0
! b1 = 0.0d0
! c1 = 0.0d0
! Parameters for H2 at equilibrium
@ -40,9 +42,9 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Parameters for stretch H2
a2 = + 0.01922622507087411d0
b2 = - 0.01799647558018601d0
c2 = - 0.022945430666782573d0
! a2 = + 0.01922622507087411d0
! b2 = - 0.01799647558018601d0
! c2 = - 0.022945430666782573d0
! Parameters for He
@ -50,6 +52,15 @@ subroutine RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! b2 = 2.715266992840757d0
! c2 = 2.1634223380633086d0
a1 = aCC_w1(1)
b1 = aCC_w1(2)
c1 = aCC_w1(3)
a2 = aCC_w2(1)
b2 = aCC_w2(2)
c2 = aCC_w2(3)
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex)
subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex)
! Compute the restricted version of the curvature-corrected exchange functional
@ -9,6 +9,8 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(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)
double precision,intent(in) :: rhow(nGrid)
@ -30,9 +32,9 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex
! Single excitation parameter
a1 = 0.0d0
b1 = 0.0d0
c1 = 0.0d0
! a1 = 0.0d0
! b1 = 0.0d0
! c1 = 0.0d0
! Parameters for H2 at equilibrium
@ -42,9 +44,9 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex
! Parameters for stretch H2
a2 = + 0.01922622507087411d0
b2 = - 0.01799647558018601d0
c2 = - 0.022945430666782573d0
! a2 = + 0.01922622507087411d0
! b2 = - 0.01799647558018601d0
! c2 = - 0.022945430666782573d0
! Parameters for He
@ -52,6 +54,19 @@ subroutine RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex
! b2 = 2.715266992840757d0
! c2 = 2.1634223380633086d0
a1 = aCC_w1(1)
b1 = aCC_w1(2)
c1 = aCC_w1(3)
a2 = aCC_w2(1)
b2 = aCC_w2(2)
c2 = aCC_w2(3)
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)

View File

@ -1,4 +1,4 @@
subroutine RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx)
! Compute the restricted version of the curvature-corrected exchange potential
@ -9,6 +9,8 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(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) :: nBas
@ -30,9 +32,9 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! Single excitation parameter
a1 = 0.0d0
b1 = 0.0d0
c1 = 0.0d0
! a1 = 0.0d0
! b1 = 0.0d0
! c1 = 0.0d0
! Parameters for H2 at equilibrium
@ -42,9 +44,9 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! Parameters for stretch H2
a2 = + 0.01922622507087411d0
b2 = - 0.01799647558018601d0
c2 = - 0.022945430666782573d0
! a2 = + 0.01922622507087411d0
! b2 = - 0.01799647558018601d0
! c2 = - 0.022945430666782573d0
! Parameters for He
@ -52,6 +54,18 @@ subroutine RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! b2 = 2.715266992840757d0
! c2 = 2.1634223380633086d0
! Parameters for He N -> N-1
a1 = aCC_w1(1)
b1 = aCC_w1(2)
c1 = aCC_w1(3)
! Parameters for He N -> N+1
a2 = aCC_w2(1)
b2 = aCC_w2(2)
c2 = aCC_w2(3)
w1 = wEns(2)
Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)

View File

@ -9,7 +9,8 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(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)
double precision,intent(in) :: rhow(nGrid)

View File

@ -238,8 +238,9 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig
! Compute exchange potential
do ispin=1,nspin
call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,Pw(:,:,ispin),ERI(:,:,:,:), &
AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin),Fx(:,:,ispin),FxHF(:,:,ispin))
call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), &
Fx(:,:,ispin),FxHF(:,:,ispin))
end do
! Compute correlation potential

View File

@ -40,7 +40,7 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_
case ('RCC')
call RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
call RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),rho(:),Ex)
case ('UCC')

View File

@ -133,72 +133,103 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,
write(*,*)
!------------------------------------------------------------------------
! Total and Excitation energies
! Total Energy and IP and EA
!------------------------------------------------------------------------
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' EXCITATION ENERGIES FROM AUXILIARY ENERGIES '
write(*,'(A60)') ' IP and EA FROM AUXILIARY ENERGIES '
write(*,'(A60)') '-------------------------------------------------'
do iEns=2,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Omaux(iEns)+OmxcDD(iEns),' au'
write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',Omaux(2)+OmxcDD(2),' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au'
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2), ' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2),' au'
write(*,*)
end do
write(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3 :',Omaux(3)+OmxcDD(3),' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(3), ' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3),' au'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
do iEns=2,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',(Omaux(iEns)+OmxcDD(iEns))*HaToeV,' eV'
write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',(Omaux(2)+OmxcDD(2))*HaToeV,' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV'
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)*HaToeV,' eV'
write(*,*)
end do
write(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3 :',(Omaux(3)+OmxcDD(3))*HaToeV,' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3)*HaToeV,' eV'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' EXCITATION ENERGIES FROM INDIVIDUAL ENERGIES '
write(*,'(A60)') ' IP and EA FROM INDIVIDUAL ENERGIES '
write(*,'(A60)') '-------------------------------------------------'
do iEns=1,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au'
end do
write(*,'(A60)') '-------------------------------------------------'
do iEns=2,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au'
write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',Om(2), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2), ' au'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2), ' au'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2),' au'
write(*,*)
end do
write(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3 :',Om(3), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(3), ' au'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(3), ' au'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(3), ' au'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3),' au'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
do iEns=2,nEns
write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV, ' eV'
write(*,'(A40,F16.10,A3)') ' Ionization Potential 1 -> 2 :',Om(2)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(2)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV'
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(2)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(2)*HaToeV,' eV'
write(*,*)
end do
write(*,'(A40,F16.10,A3)') ' Electronic Affinity 1 -> 3 :',Om(3)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(3)*HaToeV, ' eV'
write(*,*)
write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(3)*HaToeV, ' eV'
write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(3)*HaToeV,' eV'
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,*)
end subroutine print_unrestricted_individual_energy

View File

@ -108,6 +108,18 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC
read(1,*) (aCC_w1(I),I=1,3)
read(1,*) (aCC_w2(I),I=1,3)
write(*,*)'----------------------------------------------------------'
write(*,*)' parameters for w1-dependant exchange functional coefficient '
write(*,*)'----------------------------------------------------------'
call matout(3,1,aCC_w1)
write(*,*)
write(*,*)'----------------------------------------------------------'
write(*,*)' parameters for w2-dependant exchange functional coefficient '
write(*,*)'----------------------------------------------------------'
call matout(3,1,aCC_w2)
write(*,*)
write(*,*)'----------------------------------------------------------'
write(*,*)' Ensemble weights '
write(*,*)'----------------------------------------------------------'