R and U density matrices

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

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

input/options

@@ -9,8 +9,8 @@
 # GF: maxSCF thresh  DIIS n_diis renormalization
       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

scan_LiH.sh

@@ -2,8 +2,8 @@
 
 MOL="LiH"
 BASIS="cc-pvqz"
-R_START=3.5
+R_START=2.5
-R_END=3.5
+R_END=3.6
 DR=0.1
 
 for R in $(seq $R_START $DR $R_END)

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
     

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
 

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)
+  Cx0   = -(4d0/3d0)*(1d0/pi)**(1d0/3d0)
-  Cx2   = -(176d0/105d0)*(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
 

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

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

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

src/eDFT/density_matrix.f90

@@ -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

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)

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)
 

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 !!!')

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