From 3e72d87adfdb85fc097554a0c5eb01944250f0b7 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 28 Mar 2020 22:52:45 +0100 Subject: [PATCH] starting implementing B88 --- include/parameters.h | 1 + input/basis | 69 ++++++++++++++++++------ input/dft | 8 +-- input/weight | 69 ++++++++++++++++++------ src/IntPak/IntPak.f90 | 2 +- src/QuAcK/G0T0.f90 | 50 ++++++++++++++++- src/QuAcK/excitation_density_Tmatrix.f90 | 42 +++++++-------- src/QuAcK/self_energy_Tmatrix_diag.f90 | 4 -- src/QuAcK/soG0T0.f90 | 5 +- src/eDFT/B88_gga_exchange_energy.f90 | 9 ++-- src/eDFT/G96_gga_exchange_energy.f90 | 5 +- src/eDFT/G96_gga_exchange_potential.f90 | 5 +- src/eDFT/LIM_RKS.f90 | 18 ++++--- src/eDFT/RB88_gga_exchange_energy.f90 | 53 ++++++++++++++++++ src/eDFT/RB88_gga_exchange_potential.f90 | 66 +++++++++++++++++++++++ src/eDFT/gga_exchange_energy.f90 | 10 +++- src/eDFT/gga_exchange_potential.f90 | 10 +++- 17 files changed, 332 insertions(+), 94 deletions(-) create mode 100644 src/eDFT/RB88_gga_exchange_energy.f90 create mode 100644 src/eDFT/RB88_gga_exchange_potential.f90 diff --git a/include/parameters.h b/include/parameters.h index 279fe1e..e5bd88d 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -18,5 +18,6 @@ double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) +! double precision,parameter :: Cx1 = - 0.904d0*(4d0/3d0)*(1d0/pi)**(1d0/3d0) double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) diff --git a/input/basis b/input/basis index a1f7a8d..1ce59ba 100644 --- a/input/basis +++ b/input/basis @@ -1,27 +1,62 @@ -1 5 +1 14 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 1 0.8380000 1.0000000 +P 1 + 1 0.2920000 1.0000000 +P 1 + 1 0.0848000 1.0000000 +D 1 + 1 2.0620000 1.0000000 +D 1 + 1 0.6620000 1.0000000 +D 1 + 1 0.1900000 1.0000000 +F 1 + 1 1.3970000 1.0000000 +F 1 + 1 0.3600000 1.0000000 +2 14 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 - + 1 0.8380000 1.0000000 +P 1 + 1 0.2920000 1.0000000 +P 1 + 1 0.0848000 1.0000000 +D 1 + 1 2.0620000 1.0000000 +D 1 + 1 0.6620000 1.0000000 +D 1 + 1 0.1900000 1.0000000 +F 1 + 1 1.3970000 1.0000000 +F 1 + 1 0.3600000 1.0000000 diff --git a/input/dft b/input/dft index 1e4c470..b786990 100644 --- a/input/dft +++ b/input/dft @@ -3,10 +3,10 @@ # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 -# GGA = 2: +# GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 666 HF + 1 RMFL20 # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 @@ -19,6 +19,6 @@ # Number of states in ensemble (nEns) 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.00000 0.00000 + 1.00000 0.00000 # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type - 32 0.00001 T 5 1 1 + 32 0.00001 T 5 1 1 diff --git a/input/weight b/input/weight index a1f7a8d..1ce59ba 100644 --- a/input/weight +++ b/input/weight @@ -1,27 +1,62 @@ -1 5 +1 14 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 1 0.8380000 1.0000000 +P 1 + 1 0.2920000 1.0000000 +P 1 + 1 0.0848000 1.0000000 +D 1 + 1 2.0620000 1.0000000 +D 1 + 1 0.6620000 1.0000000 +D 1 + 1 0.1900000 1.0000000 +F 1 + 1 1.3970000 1.0000000 +F 1 + 1 0.3600000 1.0000000 +2 14 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 +S 1 + 1 0.0236300 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 - + 1 0.8380000 1.0000000 +P 1 + 1 0.2920000 1.0000000 +P 1 + 1 0.0848000 1.0000000 +D 1 + 1 2.0620000 1.0000000 +D 1 + 1 0.6620000 1.0000000 +D 1 + 1 0.1900000 1.0000000 +F 1 + 1 1.3970000 1.0000000 +F 1 + 1 0.3600000 1.0000000 diff --git a/src/IntPak/IntPak.f90 b/src/IntPak/IntPak.f90 index 6973651..9310168 100644 --- a/src/IntPak/IntPak.f90 +++ b/src/IntPak/IntPak.f90 @@ -311,7 +311,7 @@ program IntPak allocate(DG(1:KG),ExpG(1:KG)) DG = (/ 1d0 /) ExpG = (/ 1d0 /) - ExpS = ExpS*ExpS +! ExpS = ExpS*ExpS call cpu_time(start_2eInt(iType)) call Compute2eInt(debug,chemist_notation,iType,nShell, & diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 17eada4..c74aca7 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -102,7 +102,10 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute excitation densities for the T-matrix - call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & + rho1s(:,:,:) = 0d0 + rho2s(:,:,:) = 0d0 + + call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) !---------------------------------------------- @@ -124,7 +127,50 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute excitation densities for the T-matrix - call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & + rho1t(:,:,:) = 0d0 + rho2t(:,:,:) = 0d0 + + call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & + X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) + +!---------------------------------------------- +! Compute T-matrix version of the self-energy +!---------------------------------------------- + + SigT(:) = 0d0 + rho2s(:,:,:) = 0d0 + rho2t(:,:,:) = 0d0 + + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & + Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & + Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & + SigT(:)) +!---------------------------------------------- +! Singlet manifold +!---------------------------------------------- + + ispin = 1 + +! Compute excitation densities for the T-matrix + + rho1s(:,:,:) = 0d0 + rho2s(:,:,:) = 0d0 + + call excitation_density_Tmatrix(ispin,0d0,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & + X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) + +!---------------------------------------------- +! Triplet manifold +!---------------------------------------------- + + ispin = 2 + +! Compute excitation densities for the T-matrix + + rho1t(:,:,:) = 0d0 + rho2t(:,:,:) = 0d0 + + call excitation_density_Tmatrix(ispin,0d0,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) !---------------------------------------------- diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index f0a8aa3..fd900fc 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -1,4 +1,4 @@ -subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) +subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) ! Compute excitation densities for T-matrix self-energy @@ -7,6 +7,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r ! Input variables integer,intent(in) :: ispin + double precision,intent(in) :: db integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: X1(nVV,nVV) @@ -27,9 +28,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r double precision,intent(out) :: rho1(nBas,nO,nVV) double precision,intent(out) :: rho2(nBas,nV,nOO) - rho1(:,:,:) = 0d0 - rho2(:,:,:) = 0d0 - !---------------------------------------------- ! Singlet manifold !---------------------------------------------- @@ -45,10 +43,11 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do c=nO+1,nBas-nR do d=c,nBas-nR cd = cd + 1 -! rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,c,d)*X1(cd,ab) rho1(p,i,ab) = rho1(p,i,ab) & - + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) & - /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) +! + db*(ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) & +! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) & + + 0d0*db*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) & + + (1d0 - db)*ERI(p,i,c,d)*X1(cd,ab) end do end do @@ -56,10 +55,11 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do k=nC+1,nO do l=k,nO kl = kl + 1 -! rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,k,l)*Y1(kl,ab) rho1(p,i,ab) = rho1(p,i,ab) & - + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) & - /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) +! + db*(ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) & +! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) & + + 0d0*db*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) & + + (1d0 - db)*ERI(p,i,k,l)*Y1(kl,ab) end do end do @@ -73,10 +73,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do c=nO+1,nBas-nR do d=c,nBas-nR cd = cd + 1 -! rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,c,d)*X2(cd,ij) - rho2(p,a,ij) = rho2(p,a,ij) & - + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) & - /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d))) + rho2(p,a,ij) = rho2(p,a,ij) + db*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) end do end do @@ -84,10 +81,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do k=nC+1,nO do l=k,nO kl = kl + 1 -! rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,k,l)*Y2(kl,ij) - rho2(p,a,ij) = rho2(p,a,ij) & - + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) & - /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l))) + rho2(p,a,ij) = rho2(p,a,ij) + db*(ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) end do end do @@ -113,7 +107,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do c=nO+1,nBas-nR do d=c+1,nBas-nR cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) + 1.5d0*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) + rho1(p,i,ab) = rho1(p,i,ab) & + + 1.0d0*db*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) & + + (1d0-db)*0d0*(ERI(p,i,c,d))*X1(cd,ab) end do end do @@ -121,7 +117,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do k=nC+1,nO do l=k+1,nO kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) + 1.5d0*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) + rho1(p,i,ab) = rho1(p,i,ab) & + + 1.0d0*db*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) & + + (1d0-db)*0d0*(ERI(p,i,k,l))*Y1(kl,ab) end do end do @@ -135,7 +133,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do c=nO+1,nBas-nR do d=c+1,nBas-nR cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) + 0.5d0*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) + rho2(p,a,ij) = rho2(p,a,ij) + 1d0*(ERI(p,nO+a,c,d) - db*ERI(p,nO+a,d,c))*X2(cd,ij) end do end do @@ -143,7 +141,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do k=nC+1,nO do l=k+1,nO kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) + 0.5d0*(ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) + rho2(p,a,ij) = rho2(p,a,ij) + 1d0*(ERI(p,nO+a,k,l) - db*ERI(p,nO+a,l,k))*Y2(kl,ij) end do end do diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90 index 2ad76c8..3882dd1 100644 --- a/src/QuAcK/self_energy_Tmatrix_diag.f90 +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -32,10 +32,6 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, double precision,intent(out) :: SigT(nBas) -! Initialize - - SigT(:) = 0d0 - !---------------------------------------------- ! Singlet part of the T-matrix self-energy !---------------------------------------------- diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index 30c9f00..88793f9 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -88,7 +88,10 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) ! Compute excitation densities for the T-matrix - call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), & + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + call excitation_density_Tmatrix(ispin,1d0,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), & X1(:,:),Y1(:,:),rho1(:,:,:),X2(:,:),Y2(:,:),rho2(:,:,:)) !---------------------------------------------- diff --git a/src/eDFT/B88_gga_exchange_energy.f90 b/src/eDFT/B88_gga_exchange_energy.f90 index 439c220..d47fa7a 100644 --- a/src/eDFT/B88_gga_exchange_energy.f90 +++ b/src/eDFT/B88_gga_exchange_energy.f90 @@ -1,6 +1,6 @@ -subroutine B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) +subroutine B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) -! Compute Becke's 96 GGA exchange energy +! Compute Becke's 88 GGA exchange energy implicit none @@ -8,9 +8,6 @@ subroutine B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) ! Input variables - character(len=12),intent(in) :: DFA - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) @@ -19,7 +16,7 @@ subroutine B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) ! Local variables integer :: iG - double precision :: alpha,beta,gamma + double precision :: alpha,beta double precision :: r,g,x ! Output variables diff --git a/src/eDFT/G96_gga_exchange_energy.f90 b/src/eDFT/G96_gga_exchange_energy.f90 index b9cb6c4..90aa536 100644 --- a/src/eDFT/G96_gga_exchange_energy.f90 +++ b/src/eDFT/G96_gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine G96_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) +subroutine G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Compute Gill's 96 GGA exchange energy @@ -8,9 +8,6 @@ subroutine G96_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) ! Input variables - character(len=12),intent(in) :: DFA - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) diff --git a/src/eDFT/G96_gga_exchange_potential.f90 b/src/eDFT/G96_gga_exchange_potential.f90 index a164760..2243865 100644 --- a/src/eDFT/G96_gga_exchange_potential.f90 +++ b/src/eDFT/G96_gga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine G96_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Compute Gill's GGA exchange poential @@ -7,9 +7,6 @@ subroutine G96_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho ! Input variables - character(len=12),intent(in) :: DFA - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: nBas diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 index 6158264..25c4976 100644 --- a/src/eDFT/LIM_RKS.f90 +++ b/src/eDFT/LIM_RKS.f90 @@ -73,18 +73,16 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres !------------------------------------------------------------------------ write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' EQUI-WEIGHT CALCULATION ' + write(*,'(A40)') ' NON-ZERO-WEIGHT CALCULATION ' write(*,'(A40)') '*************************************************' - wLIM(1:nEns) = 1d0/dble(nEns) - do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns) + write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wEns(iEns) end do write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,nEns,wLIM,nGrid,weight,maxSCF,thresh, & + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,EwGICEW,F) !------------------------------------------------------------------------ @@ -92,10 +90,14 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres !------------------------------------------------------------------------ Om(1) = 0d0 - Om(2) = 2d0*(EwEW - EwZW) - OmGIC(1) = 0d0 - OmGIC(2) = 2d0*(EwGICEW - EwGICZW) + + Om(2) = 0d0 + OmGIC(2) = 0d0 + if(wEns(2) > 10d-3) then + Om(2) = (EwEW - EwZW)/wEns(2) + OmGIC(2) = (EwGICEW - EwGICZW)/wEns(2) + end if write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES ' diff --git a/src/eDFT/RB88_gga_exchange_energy.f90 b/src/eDFT/RB88_gga_exchange_energy.f90 new file mode 100644 index 0000000..7f7d48a --- /dev/null +++ b/src/eDFT/RB88_gga_exchange_energy.f90 @@ -0,0 +1,53 @@ +subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + +! Compute restricted version of Becke's 88 GGA 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) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + integer :: iG + double precision :: alpha + double precision :: beta + double precision :: r,g,x + +! Output variables + + double precision :: Ex + +! Coefficients for B88 GGA exchange functional + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + beta = 0.0042d0 + +! Compute GGA exchange energy + + Ex = 0d0 + + do iG=1,nGrid + + r = max(0d0,0.5d0*rho(iG)) + + if(r > threshold) then + g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 + x = sqrt(g)/r**(4d0/3d0) + + Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) & + - weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x)) + + end if + + end do + + Ex = 2d0*Ex + +end subroutine RB88_gga_exchange_energy diff --git a/src/eDFT/RB88_gga_exchange_potential.f90 b/src/eDFT/RB88_gga_exchange_potential.f90 new file mode 100644 index 0000000..48466f6 --- /dev/null +++ b/src/eDFT/RB88_gga_exchange_potential.f90 @@ -0,0 +1,66 @@ +subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + +! Compute restricted Becke's GGA 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) :: dAO(ncart,nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + integer :: mu,nu,iG + double precision :: alpha + double precision :: beta + double precision :: r,g,vAO,gAO + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Coefficients for B88 GGA exchange functional + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + beta = 0.0042d0 + +! Compute GGA exchange matrix in the AO basis + + Fx(:,:) = 0d0 + + do mu=1,nBas + do nu=1,nBas + do iG=1,nGrid + + r = max(0d0,0.5d0*rho(iG)) + + if(r > threshold) then + + g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2 + vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) + Fx(mu,nu) = Fx(mu,nu) & + + vAO*(4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) & + + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0)) + + gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & + + drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) & + + drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG)) + + gAO = weight(iG)*gAO + + Fx(mu,nu) = Fx(mu,nu) - 2d0*gAO*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0) + + end if + + end do + end do + end do + +end subroutine RB88_gga_exchange_potential diff --git a/src/eDFT/gga_exchange_energy.f90 b/src/eDFT/gga_exchange_energy.f90 index 06d72dc..b110a5b 100644 --- a/src/eDFT/gga_exchange_energy.f90 +++ b/src/eDFT/gga_exchange_energy.f90 @@ -26,13 +26,19 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) case ('G96') - call G96_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) + call G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + +! Becke's 88 exchange functional + + case ('RB88') + + call RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Becke's 88 exchange functional case ('B88') - call B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) + call B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) case default diff --git a/src/eDFT/gga_exchange_potential.f90 b/src/eDFT/gga_exchange_potential.f90 index 538b4a0..78f94a2 100644 --- a/src/eDFT/gga_exchange_potential.f90 +++ b/src/eDFT/gga_exchange_potential.f90 @@ -30,13 +30,19 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh case ('G96') - call G96_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +! Becke's 88 exchange functional + + case ('RB88') + + call RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + ! Becke's 88 exchange functional case ('B88') - call B88_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) case default