10
1
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00

huckel guess in eDFT

This commit is contained in:
Pierre-Francois Loos 2020-03-29 21:35:40 +02:00
parent 438a36a82c
commit e4dcb5fadf
14 changed files with 162 additions and 89 deletions

View File

@ -1,27 +1,26 @@
1 5 1 8
S 6
1 1264.5857000 0.0019448
2 189.9368100 0.0148351
3 43.1590890 0.0720906
4 12.0986630 0.2371542
5 3.8063232 0.4691987
6 1.2728903 0.3565202
S 3 S 3
1 13.0100000 0.0196850 1 3.1964631 -0.1126487
2 1.9620000 0.1379770 2 0.7478133 -0.2295064
3 0.4446000 0.4781480 3 0.2199663 1.1869167
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1 S 1
1 0.1220000 1.0000000 1 0.0823099 1.0000000
P 1
1 0.0823099 1.0000000
S 1 S 1
1 0.0297400 1.0000000 1 0.0207000 1.0000000
P 1 P 1
1 0.7270000 1.0000000 1 0.0207000 1.0000000
P 1 D 1
1 0.1410000 1.0000000 1 0.4000000 1.0000000
2 5
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.1220000 1.0000000
S 1
1 0.0297400 1.0000000
P 1
1 0.7270000 1.0000000
P 1
1 0.1410000 1.0000000

View File

@ -13,12 +13,12 @@
# GGA = 2: # GGA = 2:
# Hybrid = 4: # Hybrid = 4:
# Hartree-Fock = 666 # Hartree-Fock = 666
1 RVWN5 1 RMFL20
# quadrature grid SG-n # quadrature grid SG-n
1 1
# Number of states in ensemble (nEns) # Number of states in ensemble (nEns)
2 2
# Ensemble weights: wEns(1),...,wEns(nEns-1) # Ensemble weights: wEns(1),...,wEns(nEns-1)
0.5 0.0 1.0 0.0
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type
32 0.00001 T 5 1 1 32 0.00001 T 5 2 1

View File

@ -1,5 +1,4 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
2 1 1 0 0 1 2 2 0 0
# Znuc x y z # Znuc x y z
H 0.0 0.0 0.0 Be 0.0 0.0 0.0
H 0.0 0.0 1.4

View File

@ -1,4 +1,3 @@
2 1
H 0.0000000000 0.0000000000 0.0000000000 Be 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 0.7408481486

View File

@ -1,27 +1,26 @@
1 5 1 8
S 6
1 1264.5857000 0.0019448
2 189.9368100 0.0148351
3 43.1590890 0.0720906
4 12.0986630 0.2371542
5 3.8063232 0.4691987
6 1.2728903 0.3565202
S 3 S 3
1 13.0100000 0.0196850 1 3.1964631 -0.1126487
2 1.9620000 0.1379770 2 0.7478133 -0.2295064
3 0.4446000 0.4781480 3 0.2199663 1.1869167
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1 S 1
1 0.1220000 1.0000000 1 0.0823099 1.0000000
P 1
1 0.0823099 1.0000000
S 1 S 1
1 0.0297400 1.0000000 1 0.0207000 1.0000000
P 1 P 1
1 0.7270000 1.0000000 1 0.0207000 1.0000000
P 1 D 1
1 0.1410000 1.0000000 1 0.4000000 1.0000000
2 5
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
S 1
1 0.1220000 1.0000000
S 1
1 0.0297400 1.0000000
P 1
1 0.7270000 1.0000000
P 1
1 0.1410000 1.0000000

View File

@ -128,17 +128,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS
! Guess coefficients and eigenvalues ! Guess coefficients and eigenvalues
if(.not. restart) then if(.not. restart) call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,Fx,X,cp,F,Fp,eps,c,Pw)
if(guess_type == 1) then
F(:,:) = Hc(:,:)
else if(guess_type == 2) then
call random_number(F(:,:))
end if
end if
! Initialization ! Initialization

View File

@ -16,7 +16,7 @@ subroutine MFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,r
! Local variables ! Local variables
logical :: LDA_centered = .false. logical :: LDA_centered = .true.
integer :: iEns,isp integer :: iEns,isp
double precision :: EcLDA(nsp) double precision :: EcLDA(nsp)
double precision,allocatable :: aMFL(:,:) double precision,allocatable :: aMFL(:,:)

View File

@ -43,9 +43,7 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec)
! Compute correlation energy for ground and doubly-excited states ! Compute correlation energy for ground and doubly-excited states
do iEns=1,nEns do iEns=1,nEns
call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rho(:),EceLDA(iEns)) call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rho(:),EceLDA(iEns))
end do end do
! LDA-centered functional ! LDA-centered functional
@ -61,11 +59,8 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec)
! Weight-denpendent functional for ensembles ! Weight-denpendent functional for ensembles
Ec = 0d0 Ec = 0d0
do iEns=1,nEns do iEns=1,nEns
Ec = Ec + wEns(iEns)*EceLDA(iEns) Ec = Ec + wEns(iEns)*EceLDA(iEns)
end do end do
end subroutine RMFL20_lda_correlation_energy end subroutine RMFL20_lda_correlation_energy

View File

@ -48,8 +48,6 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,
! LDA-centered functional ! LDA-centered functional
EcLDA = 0d0
call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA) call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA)
if(LDA_centered) then if(LDA_centered) then
@ -61,7 +59,6 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,
! Weight-denpendent functional for ensembles ! Weight-denpendent functional for ensembles
Ec = 0d0 Ec = 0d0
do iEns=1,nEns do iEns=1,nEns
Ec = Ec + wEns(iEns)*EceLDA(iEns) Ec = Ec + wEns(iEns)*EceLDA(iEns)
enddo enddo

View File

@ -44,15 +44,11 @@ include 'parameters.h'
! Compute correlation energy for ground, singly-excited and doubly-excited states ! Compute correlation energy for ground, singly-excited and doubly-excited states
do iEns=1,nEns do iEns=1,nEns
call restricted_elda_correlation_potential(nEns,aMFL(:,iEns),nGrid,weight,nBas,AO,rho,FceLDA(:,:,iEns)) call restricted_elda_correlation_potential(nEns,aMFL(:,iEns),nGrid,weight,nBas,AO,rho,FceLDA(:,:,iEns))
end do end do
! LDA-centered functional ! LDA-centered functional
FcLDA(:,:) = 0d0
call RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,FcLDA) call RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,FcLDA)
if(LDA_centered) then if(LDA_centered) then
@ -64,11 +60,8 @@ include 'parameters.h'
! Weight-denpendent functional for ensembles ! Weight-denpendent functional for ensembles
Fc(:,:) = 0d0 Fc(:,:) = 0d0
do iEns=1,nEns do iEns=1,nEns
Fc(:,:) = Fc(:,:) + wEns(iEns)*FceLDA(:,:,iEns) Fc(:,:) = Fc(:,:) + wEns(iEns)*FceLDA(:,:,iEns)
enddo enddo
end subroutine RMFL20_lda_correlation_potential end subroutine RMFL20_lda_correlation_potential

View File

@ -39,15 +39,10 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
Ex = 0d0 Ex = 0d0
do iG=1,nGrid do iG=1,nGrid
r = max(0d0,rho(iG)) r = max(0d0,rho(iG))
if(r > threshold) then if(r > threshold) then
Ex = Ex + weight(iG)*Cxw*r**(4d0/3d0) Ex = Ex + weight(iG)*Cxw*r**(4d0/3d0)
endif endif
enddo enddo
end subroutine RMFL20_lda_exchange_energy end subroutine RMFL20_lda_exchange_energy

View File

@ -45,11 +45,9 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho
rI = max(0d0,rho(iG)) rI = max(0d0,rho(iG))
if(r > threshold .and. rI > threshold) then if(r > threshold .and. rI > threshold) then
e = Cxw*r**(1d0/3d0) e = Cxw*r**(1d0/3d0)
dedr = 1d0/3d0*Cxw*r**(-2d0/3d0) dedr = 1d0/3d0*Cxw*r**(-2d0/3d0)
Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r)
endif endif
enddo enddo

57
src/eDFT/huckel_guess.f90 Normal file
View File

@ -0,0 +1,57 @@
subroutine huckel_guess(nBas,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
! Hickel guess of the molecular orbitals for HF calculation
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
! Local variables
integer :: mu,nu
double precision :: a
! Output variables
double precision,intent(out) :: c(nBas,nBas)
a = 1.75d0
Fp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
call hartree_coulomb(nBas,P,ERI,J)
call fock_exchange_potential(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
do mu=1,nBas
do nu=mu+1,nBas
F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
F(nu,mu) = F(mu,nu)
enddo
enddo
Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
end subroutine

52
src/eDFT/mo_guess.f90 Normal file
View File

@ -0,0 +1,52 @@
subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
! Guess of the molecular orbitals for HF calculation
implicit none
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: guess_type
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
if(guess_type == 1) then
Fp = matmul(transpose(X),matmul(Hc,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
elseif(guess_type == 2) then
call huckel_guess(nBas,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
elseif(guess_type == 3) then
call random_number(c)
else
print*,'Wrong guess option'
stop
endif
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
end subroutine