diff --git a/examples/molecule.F2 b/examples/molecule.F2 index ac3ddcd..12bba91 100644 --- a/examples/molecule.F2 +++ b/examples/molecule.F2 @@ -2,4 +2,4 @@ 2 9 9 0 0 # Znuc x y z F 0. 0. 0. - F 0. 0. 2.1 + F 0. 0. 3.3 diff --git a/input/dft b/input/dft index a462f87..2b261ca 100644 --- a/input/dft +++ b/input/dft @@ -1,14 +1,14 @@ # Restricted or unrestricted KS calculation - RKS + GOK-RKS # exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666 - 1 S51 + 1 RS51 # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666 - 0 LF19 + 0 HF # quadrature grid SG-n 1 # Number of states in ensemble (nEns) - 3 + 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.00000 0.50000 + 0.50000 0.00000 # eKS: maxSCF thresh DIIS n_diis guess_type ortho_type 64 0.0000001 T 15 1 1 diff --git a/input/options b/input/options index 8489957..bc81042 100644 --- a/input/options +++ b/input/options @@ -7,10 +7,10 @@ # CIS/TDHF/BSE: singlet triplet T T # GF: maxSCF thresh DIIS n_diis renormalization - 256 0.00001 T 5 3 + 256 0.00001 T 5 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta - 1 0.00001 T 5 F F T F F F T 0.000 + 256 0.00001 T 5 F F T F F F F 0.000 # ACFDT: AC Kx XBS - T F T + T F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/scan_LiH.sh b/scan_LiH.sh index 4643c11..38759d5 100755 --- a/scan_LiH.sh +++ b/scan_LiH.sh @@ -2,8 +2,8 @@ MOL="LiH" BASIS="cc-pvqz" -R_START=3.5 -R_END=3.5 +R_START=2.5 +R_END=3.6 DR=0.1 for R in $(seq $R_START $DR $R_END) diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index c294628..3d3f331 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -182,7 +182,7 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute density matrix !------------------------------------------------------------------------ - call density_matrix(nBas,nEns,nO,c(:,:),P(:,:,:)) + call restricted_density_matrix(nBas,nEns,nO,c(:,:),P(:,:,:)) ! Weight-dependent density matrix diff --git a/src/eDFT/GOK_UKS.f90 b/src/eDFT/GOK_UKS.f90 index 0f1b4a9..41556ad 100644 --- a/src/eDFT/GOK_UKS.f90 +++ b/src/eDFT/GOK_UKS.f90 @@ -192,7 +192,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute density matrix !------------------------------------------------------------------------ - call density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:)) + call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:)) ! Weight-dependent density matrix @@ -224,7 +224,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres do ispin=1,nspin do iEns=1,nEns -! call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns)) + call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns)) end do end do diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 index fc5cd65..33adeb3 100644 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_energy.f90 @@ -16,8 +16,9 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Local variables integer :: iG + double precision :: Cx0 + double precision :: Cx1 double precision :: CxLDA - double precision :: Cx2 double precision :: Cxw double precision :: r @@ -27,10 +28,11 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Cx coefficient for Slater LDA exchange - CxLDA = -(4d0/3d0)*(1d0/pi)**(1d0/3d0) - Cx2 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0) + Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0) + Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0) + CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) - Cxw = (1d0 - wEns(2))*CxLDA + wEns(2)*(Cx2 - CxLDA) + Cxw = CxLDA + wEns(1)*(Cx1 - Cx0) ! Compute LDA exchange energy diff --git a/src/eDFT/RMFL20_lda_exchange_potential.f90 b/src/eDFT/RMFL20_lda_exchange_potential.f90 new file mode 100644 index 0000000..eebb5d4 --- /dev/null +++ b/src/eDFT/RMFL20_lda_exchange_potential.f90 @@ -0,0 +1,59 @@ +subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + +! Compute the restricted version of the weight-dependent MFL20 exchange potential + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: mu,nu,iG + double precision :: Cx0 + double precision :: Cx1 + double precision :: CxLDA + double precision :: Cxw + double precision :: r,vAO + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Cx coefficient for Slater LDA exchange + + Cx0 = -(4d0/3d0)*(1d0/pi)**(1d0/3d0) + Cx1 = -(176d0/105d0)*(1d0/pi)**(1d0/3d0) + CxLDA = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) + + Cxw = CxLDA + wEns(1)*(Cx1 - Cx0) + +! Compute LDA exchange matrix in the AO basis + + Fx(:,:) = 0d0 + do mu=1,nBas + do nu=1,nBas + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) + Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*Cxw*r**(1d0/3d0) + + endif + + enddo + enddo + enddo + +end subroutine RMFL20_lda_exchange_potential diff --git a/src/eDFT/RS51_lda_exchange_energy.f90 b/src/eDFT/RS51_lda_exchange_energy.f90 new file mode 100644 index 0000000..86cde1a --- /dev/null +++ b/src/eDFT/RS51_lda_exchange_energy.f90 @@ -0,0 +1,42 @@ +subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex) + +! Compute restriced version of Slater's LDA exchange energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: iG + double precision :: alpha,r + +! Output variables + + double precision :: Ex + +! Cx coefficient for Slater LDA exchange + + alpha = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) + +! Compute LDA exchange energy + + Ex = 0d0 + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) + + endif + + enddo + +end subroutine RS51_lda_exchange_energy diff --git a/src/eDFT/RS51_lda_exchange_potential.f90 b/src/eDFT/RS51_lda_exchange_potential.f90 new file mode 100644 index 0000000..27a9b97 --- /dev/null +++ b/src/eDFT/RS51_lda_exchange_potential.f90 @@ -0,0 +1,50 @@ +subroutine RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) + +! Compute the restricted version of Slater's LDA exchange potential + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + +! Local variables + + integer :: mu,nu,iG + double precision :: Cx + double precision :: r,vAO + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Cx coefficient for Slater LDA exchange + + Cx = -(3d0/4d0)*(3d0/pi)**(1d0/3d0) + +! Compute LDA exchange matrix in the AO basis + + Fx(:,:) = 0d0 + do mu=1,nBas + do nu=1,nBas + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) + Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*Cx*r**(1d0/3d0) + + endif + + enddo + enddo + enddo + +end subroutine RS51_lda_exchange_potential diff --git a/src/eDFT/density_matrix.f90 b/src/eDFT/density_matrix.f90 deleted file mode 100644 index a37c1ea..0000000 --- a/src/eDFT/density_matrix.f90 +++ /dev/null @@ -1,43 +0,0 @@ -subroutine density_matrix(nBas,nEns,nO,c,P) - -! Calculate density matrices - - implicit none - - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - integer,intent(in) :: nEns - integer,intent(in) :: nO(nspin) - double precision,intent(in) :: c(nBas,nBas,nspin) - -! Local variables - - integer :: ispin - -! Output variables - - double precision,intent(out) :: P(nBas,nBas,nspin,nEns) - -! Ground state density matrix - - do ispin=1,nspin - P(:,:,ispin,1) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) - end do - -! Singly-excited state density matrix - - P(:,:,1,2) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1))) & - + matmul(c(:,nO(1)+1:nO(1)+1,1),transpose(c(:,nO(1)+1:nO(1)+1,1))) - P(:,:,2,2) = matmul(c(:,1:nO(2),2),transpose(c(:,1:nO(2),2))) - -! Doubly-excited state density matrix - - do ispin=1,nspin - P(:,:,ispin,3) = matmul(c(:,1:nO(ispin)-1,ispin),transpose(c(:,1:nO(ispin)-1,ispin))) & - + matmul(c(:,nO(ispin)+1:nO(ispin)+1,ispin),transpose(c(:,nO(ispin)+1:nO(ispin)+1,ispin))) - end do - -end subroutine density_matrix diff --git a/src/eDFT/exchange_energy.f90 b/src/eDFT/exchange_energy.f90 index 34b2497..6c0b39e 100644 --- a/src/eDFT/exchange_energy.f90 +++ b/src/eDFT/exchange_energy.f90 @@ -69,6 +69,7 @@ subroutine exchange_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho, + aX*(ExGGA - ExLDA) ! Hartree-Fock calculation + case(666) call fock_exchange_energy(nBas,P,FxHF,ExHF) diff --git a/src/eDFT/lda_exchange_energy.f90 b/src/eDFT/lda_exchange_energy.f90 index 20c3754..023e15c 100644 --- a/src/eDFT/lda_exchange_energy.f90 +++ b/src/eDFT/lda_exchange_energy.f90 @@ -22,13 +22,21 @@ subroutine lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,Ex) select case (DFA) -! Slater's LDA correlation functional: Slater, Phys. Rev. 81 (1951) 385 +! Slater's LDA correlation functional + + case ('RS51') + + call RS51_lda_exchange_energy(nGrid,weight,rho,Ex) + +! Restricted version of Slater's LDA correlation functional case ('S51') call S51_lda_exchange_energy(nGrid,weight,rho,Ex) - case ('MFL20') +! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 exchange functional + + case ('RMFL20') call RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) diff --git a/src/eDFT/lda_exchange_potential.f90 b/src/eDFT/lda_exchange_potential.f90 index 9da6444..7bc8606 100644 --- a/src/eDFT/lda_exchange_potential.f90 +++ b/src/eDFT/lda_exchange_potential.f90 @@ -25,12 +25,24 @@ subroutine lda_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) select case (DFA) -! Slater's LDA correlation functional: Slater, Phys. Rev. 81 (1951) 385 +! Restricted version of Slater's LDA correlation functional + + case ('RS51') + + call RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) + +! Slater's LDA correlation functional case ('S51') call S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) +! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 functional + + case ('RMFL20') + + call RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + case default call print_warning('!!! LDA exchange functional not available !!!') diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90 index 2a7bc5e..91c5810 100644 --- a/src/eDFT/read_options.f90 +++ b/src/eDFT/read_options.f90 @@ -14,7 +14,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th character(len=7),intent(out) :: method integer,intent(out) :: x_rung,c_rung - character(len=3),intent(out) :: x_DFA, c_DFA + character(len=12),intent(out) :: x_DFA, c_DFA integer,intent(out) :: SGn integer,intent(out) :: nEns double precision,intent(out) :: wEns(maxEns) @@ -32,7 +32,7 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th ! Open file with method specification - open(unit=40,file='input/dft') + open(unit=1,file='input/dft') ! Default values @@ -108,6 +108,6 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,th ! Close file with options - close(unit=40) + close(unit=1) end subroutine read_options diff --git a/src/eDFT/select_rung.f90 b/src/eDFT/select_rung.f90 index 3d20c3a..8ae7343 100644 --- a/src/eDFT/select_rung.f90 +++ b/src/eDFT/select_rung.f90 @@ -14,7 +14,7 @@ subroutine select_rung(rung,DFA) ! Hartree calculation case(0) - write(*,*) "* 0th rung of Jacob's ladder: Hartree calculation *" + write(*,*) "* 0th rung of Jacob's ladder: Hartree calculation *" ! LDA functionals case(1) @@ -34,7 +34,7 @@ subroutine select_rung(rung,DFA) ! Hartree-Fock calculation case(666) - write(*,*) "* rung 666: Hartree-Fock calculation *" + write(*,*) "* rung 666: Hartree-Fock calculation *" ! Default case default