diff --git a/input/basis b/input/basis index a1f7a8d..d92a3df 100644 --- a/input/basis +++ b/input/basis @@ -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 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 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.1220000 1.0000000 + 1 0.0823099 1.0000000 +P 1 + 1 0.0823099 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.0207000 1.0000000 P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 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 - + 1 0.0207000 1.0000000 +D 1 + 1 0.4000000 1.0000000 diff --git a/input/dft b/input/dft index 21ea210..d705d62 100644 --- a/input/dft +++ b/input/dft @@ -13,12 +13,12 @@ # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 - 1 RVWN5 + 1 RMFL20 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) 2 # 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 - 32 0.00001 T 5 1 1 + 32 0.00001 T 5 2 1 diff --git a/input/molecule b/input/molecule index 8076140..6a6f6d1 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 1 2 2 0 0 # Znuc x y z - H 0.0 0.0 0.0 - H 0.0 0.0 1.4 + Be 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..8023e37 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7408481486 + Be 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index a1f7a8d..d92a3df 100644 --- a/input/weight +++ b/input/weight @@ -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 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 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.1220000 1.0000000 + 1 0.0823099 1.0000000 +P 1 + 1 0.0823099 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.0207000 1.0000000 P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 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 - + 1 0.0207000 1.0000000 +D 1 + 1 0.4000000 1.0000000 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index f2adb51..a7bd8bc 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -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 - if(.not. restart) then - if(guess_type == 1) then - - F(:,:) = Hc(:,:) - - else if(guess_type == 2) then - - call random_number(F(:,:)) - - end if - end if + if(.not. restart) call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,Fx,X,cp,F,Fp,eps,c,Pw) ! Initialization diff --git a/src/eDFT/MFL20_lda_correlation_individual_energy.f90 b/src/eDFT/MFL20_lda_correlation_individual_energy.f90 index 9521aee..29217c1 100644 --- a/src/eDFT/MFL20_lda_correlation_individual_energy.f90 +++ b/src/eDFT/MFL20_lda_correlation_individual_energy.f90 @@ -16,7 +16,7 @@ subroutine MFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,r ! Local variables - logical :: LDA_centered = .false. + logical :: LDA_centered = .true. integer :: iEns,isp double precision :: EcLDA(nsp) double precision,allocatable :: aMFL(:,:) diff --git a/src/eDFT/RMFL20_lda_correlation_energy.f90 b/src/eDFT/RMFL20_lda_correlation_energy.f90 index 7dd4b13..534ca8c 100644 --- a/src/eDFT/RMFL20_lda_correlation_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_energy.f90 @@ -43,9 +43,7 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) ! Compute correlation energy for ground and doubly-excited states do iEns=1,nEns - call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rho(:),EceLDA(iEns)) - end do ! LDA-centered functional @@ -61,11 +59,8 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) ! Weight-denpendent functional for ensembles Ec = 0d0 - do iEns=1,nEns - Ec = Ec + wEns(iEns)*EceLDA(iEns) - end do end subroutine RMFL20_lda_correlation_energy diff --git a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 index bbee4af..9b0d37d 100644 --- a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 @@ -48,8 +48,6 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow, ! LDA-centered functional - EcLDA = 0d0 - call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA) if(LDA_centered) then @@ -61,7 +59,6 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow, ! Weight-denpendent functional for ensembles Ec = 0d0 - do iEns=1,nEns Ec = Ec + wEns(iEns)*EceLDA(iEns) enddo diff --git a/src/eDFT/RMFL20_lda_correlation_potential.f90 b/src/eDFT/RMFL20_lda_correlation_potential.f90 index a747690..efc53b3 100644 --- a/src/eDFT/RMFL20_lda_correlation_potential.f90 +++ b/src/eDFT/RMFL20_lda_correlation_potential.f90 @@ -44,15 +44,11 @@ include 'parameters.h' ! Compute correlation energy for ground, singly-excited and doubly-excited states do iEns=1,nEns - call restricted_elda_correlation_potential(nEns,aMFL(:,iEns),nGrid,weight,nBas,AO,rho,FceLDA(:,:,iEns)) - end do ! LDA-centered functional - FcLDA(:,:) = 0d0 - call RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,FcLDA) if(LDA_centered) then @@ -64,11 +60,8 @@ include 'parameters.h' ! Weight-denpendent functional for ensembles Fc(:,:) = 0d0 - do iEns=1,nEns - Fc(:,:) = Fc(:,:) + wEns(iEns)*FceLDA(:,:,iEns) - enddo end subroutine RMFL20_lda_correlation_potential diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 index f38d77a..092f8a9 100644 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_energy.f90 @@ -39,15 +39,10 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) Ex = 0d0 do iG=1,nGrid - r = max(0d0,rho(iG)) - if(r > threshold) then - Ex = Ex + weight(iG)*Cxw*r**(4d0/3d0) - endif - enddo end subroutine RMFL20_lda_exchange_energy diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 index 52f76ab..f35a235 100644 --- a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 @@ -45,11 +45,9 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho rI = max(0d0,rho(iG)) if(r > threshold .and. rI > threshold) then - e = Cxw*r**(1d0/3d0) dedr = 1d0/3d0*Cxw*r**(-2d0/3d0) Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) - endif enddo diff --git a/src/eDFT/huckel_guess.f90 b/src/eDFT/huckel_guess.f90 new file mode 100644 index 0000000..c15b0fa --- /dev/null +++ b/src/eDFT/huckel_guess.f90 @@ -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 diff --git a/src/eDFT/mo_guess.f90 b/src/eDFT/mo_guess.f90 new file mode 100644 index 0000000..0a6e600 --- /dev/null +++ b/src/eDFT/mo_guess.f90 @@ -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