10
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00

shuffle SCF

This commit is contained in:
Pierre-Francois Loos 2020-03-30 17:50:07 +02:00
parent e4dcb5fadf
commit f7fa1de9b3
8 changed files with 68 additions and 77 deletions

View File

@ -1,26 +1,11 @@
1 8 1 1
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 3.1964631 -0.1126487 1 3.42525091 0.15432897
2 0.7478133 -0.2295064 2 0.62391373 0.53532814
3 0.2199663 1.1869167 3 0.16885540 0.44463454
P 3 2 1
1 3.1964631 0.0559802 S 3
2 0.7478133 0.2615506 1 3.42525091 0.15432897
3 0.2199663 0.7939723 2 0.62391373 0.53532814
S 1 3 0.16885540 0.44463454
1 0.0823099 1.0000000
P 1
1 0.0823099 1.0000000
S 1
1 0.0207000 1.0000000
P 1
1 0.0207000 1.0000000
D 1
1 0.4000000 1.0000000

View File

@ -1,5 +1,5 @@
# Restricted or unrestricted KS calculation # Restricted or unrestricted KS calculation
LIM-RKS GOK-RKS
# exchange rung: # exchange rung:
# Hartree = 0 # Hartree = 0
# LDA = 1: RS51,RMFL20 # LDA = 1: RS51,RMFL20
@ -19,6 +19,6 @@
# 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)
1.0 0.0 0.0000 0.00000
# 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 2 1 32 0.00001 T 5 1 1

View File

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

View File

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

View File

@ -1,26 +1,11 @@
1 8 1 1
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 3.1964631 -0.1126487 1 3.42525091 0.15432897
2 0.7478133 -0.2295064 2 0.62391373 0.53532814
3 0.2199663 1.1869167 3 0.16885540 0.44463454
P 3 2 1
1 3.1964631 0.0559802 S 3
2 0.7478133 0.2615506 1 3.42525091 0.15432897
3 0.2199663 0.7939723 2 0.62391373 0.53532814
S 1 3 0.16885540 0.44463454
1 0.0823099 1.0000000
P 1
1 0.0823099 1.0000000
S 1
1 0.0207000 1.0000000
P 1
1 0.0207000 1.0000000
D 1
1 0.4000000 1.0000000

View File

@ -21,7 +21,8 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS
double precision,intent(in) :: AO(nBas,nGrid) double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(ncart,nBas,nGrid) double precision,intent(in) :: dAO(ncart,nBas,nGrid)
integer,intent(in) :: nO,nV integer,intent(in) :: nO
integer,intent(in) :: nV
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: V(nBas,nBas)
@ -128,7 +129,21 @@ 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) call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,Fx,X,cp,F,Fp,eps,c,Pw) if(.not. restart) then
if(guess_type == 1) then
cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
call diagonalize_matrix(nBas,cp(:,:),eps(:))
c(:,:) = matmul(X(:,:),cp(:,:))
else
print*,'Wrong guess option'
stop
end if
end if
! Initialization ! Initialization
@ -163,19 +178,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS
nSCF = nSCF + 1 nSCF = nSCF + 1
! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp(:,:),eps(:))
! Back-transform eigenvectors in non-orthogonal basis
c(:,:) = matmul(X(:,:),cp(:,:))
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Compute density matrix ! Compute density matrix
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -254,6 +256,19 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS
if(abs(rcond) < 1d-15) n_diis = 0 if(abs(rcond) < 1d-15) n_diis = 0
! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp(:,:),eps(:))
! Back-transform eigenvectors in non-orthogonal basis
c(:,:) = matmul(X(:,:),cp(:,:))
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Compute KS energy ! Compute KS energy
!------------------------------------------------------------------------ !------------------------------------------------------------------------

View File

@ -166,7 +166,7 @@ program eDFT
call cpu_time(end_KS) call cpu_time(end_KS)
t_KS = end_KS - start_KS t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOC-RKS = ',t_KS,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOK-RKS = ',t_KS,' seconds'
write(*,*) write(*,*)
end if end if

View File

@ -49,6 +49,10 @@ subroutine print_RKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ', ENuc,' au' write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ', ENuc,' au'
write(*,'(A40,1X,F16.10,A3)') ' Kohn-Sham energy: ',Ew + ENuc,' au' write(*,'(A40,1X,F16.10,A3)') ' Kohn-Sham energy: ',Ew + ENuc,' au'
write(*,'(A60)') '-------------------------------------------------' 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 HOMO energy:',eps(HOMO)*HatoeV,' eV'
write(*,'(A40,F13.6,A3)') ' KS LUMO energy:',eps(LUMO)*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(*,'(A40,F13.6,A3)') ' KS HOMO-LUMO gap:',Gap *HatoeV,' eV'