From f7fa1de9b3efb8ac54c1ac5c0667008bb74dcbfe Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 30 Mar 2020 17:50:07 +0200 Subject: [PATCH] shuffle SCF --- input/basis | 35 ++++++++++---------------------- input/dft | 6 +++--- input/molecule | 5 +++-- input/molecule.xyz | 5 +++-- input/weight | 35 ++++++++++---------------------- src/eDFT/GOK_RKS.f90 | 45 ++++++++++++++++++++++++++++-------------- src/eDFT/eDFT.f90 | 2 +- src/eDFT/print_RKS.f90 | 12 +++++++---- 8 files changed, 68 insertions(+), 77 deletions(-) diff --git a/input/basis b/input/basis index d92a3df..a59a2e4 100644 --- a/input/basis +++ b/input/basis @@ -1,26 +1,11 @@ -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 +1 1 S 3 - 1 3.1964631 -0.1126487 - 2 0.7478133 -0.2295064 - 3 0.2199663 1.1869167 -P 3 - 1 3.1964631 0.0559802 - 2 0.7478133 0.2615506 - 3 0.2199663 0.7939723 -S 1 - 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 + 1 3.42525091 0.15432897 + 2 0.62391373 0.53532814 + 3 0.16885540 0.44463454 +2 1 +S 3 + 1 3.42525091 0.15432897 + 2 0.62391373 0.53532814 + 3 0.16885540 0.44463454 + diff --git a/input/dft b/input/dft index d705d62..35b19de 100644 --- a/input/dft +++ b/input/dft @@ -1,5 +1,5 @@ # Restricted or unrestricted KS calculation - LIM-RKS + GOK-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 @@ -19,6 +19,6 @@ # Number of states in ensemble (nEns) 2 # 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 - 32 0.00001 T 5 2 1 + 32 0.00001 T 5 1 1 diff --git a/input/molecule b/input/molecule index 6a6f6d1..8076140 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,5 @@ # nAt nEla nElb nCore nRyd - 1 2 2 0 0 + 2 1 1 0 0 # 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 diff --git a/input/molecule.xyz b/input/molecule.xyz index 8023e37..6edc99d 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -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 diff --git a/input/weight b/input/weight index d92a3df..a59a2e4 100644 --- a/input/weight +++ b/input/weight @@ -1,26 +1,11 @@ -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 +1 1 S 3 - 1 3.1964631 -0.1126487 - 2 0.7478133 -0.2295064 - 3 0.2199663 1.1869167 -P 3 - 1 3.1964631 0.0559802 - 2 0.7478133 0.2615506 - 3 0.2199663 0.7939723 -S 1 - 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 + 1 3.42525091 0.15432897 + 2 0.62391373 0.53532814 + 3 0.16885540 0.44463454 +2 1 +S 3 + 1 3.42525091 0.15432897 + 2 0.62391373 0.53532814 + 3 0.16885540 0.44463454 + diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index a7bd8bc..0c70040 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -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) :: 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) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) @@ -128,8 +129,22 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS ! 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 nSCF = 0 @@ -163,19 +178,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS 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 !------------------------------------------------------------------------ @@ -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 +! 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 !------------------------------------------------------------------------ diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index e8dd7f2..842b0cb 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -166,7 +166,7 @@ program eDFT call cpu_time(end_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(*,*) end if diff --git a/src/eDFT/print_RKS.f90 b/src/eDFT/print_RKS.f90 index 346d8b6..509559d 100644 --- a/src/eDFT/print_RKS.f90 +++ b/src/eDFT/print_RKS.f90 @@ -36,7 +36,7 @@ subroutine print_RKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) 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)') ' 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)') '-------------------------------------------------' @@ -45,13 +45,17 @@ subroutine print_RKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) 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)') ' 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(*,'(A40,F13.6,A3)') ' KS HOMO-LUMO gap:',Gap *HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' write(*,*)