working on structure eDFT code

src/eDFT/select_rung.f90

@@ -32,10 +32,6 @@ subroutine select_rung(rung,DFA)
     case(4)
       write(*,*) "* 4th rung of Jacob's ladder: hybrid functional *"
 
-! Hartree-Fock calculation
-    case(666)
-      write(*,*) "* rung 666: Hartree-Fock calculation                              *"
-
 ! Default    
     case default
       write(*,*) "!!! rung not available !!!"

src/eDFT/unrestricted_correlation_derivative_discontinuity.f90

@@ -45,6 +45,13 @@ subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,
       call print_warning('!!! derivative discontinuity NYI for GGAs !!!')
       stop
 
+!   MGGA functionals
+
+    case(3) 
+
+      call print_warning('!!! derivative discontinuity NYI for MGGAs !!!')
+      stop
+
 !   Hybrid functionals
 
     case(4) 
@@ -54,12 +61,6 @@ subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,
 
       aC = 0.81d0
 
-!   Hartree-Fock calculation
-
-    case(666) 
-
-      Ec(:,:) = 0d0
-
   end select
  
 end subroutine unrestricted_correlation_derivative_discontinuity

src/eDFT/unrestricted_correlation_energy.f90

@@ -46,6 +46,12 @@ subroutine unrestricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,d
 
       call unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
 
+!   MGGA functionals
+
+    case(3) 
+
+      call unrestricted_mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
+
 !   Hybrid functionals
 
     case(4) 
@@ -57,12 +63,6 @@ subroutine unrestricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,d
 
       Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) 
 
-!   Hartree-Fock calculation
-
-    case(666) 
-
-      Ec(:) = 0d0
-
   end select
  
 end subroutine unrestricted_correlation_energy

src/eDFT/unrestricted_correlation_individual_energy.f90

@@ -50,6 +50,13 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
       call print_warning('!!! Individual energies NYI for GGAs !!!')
       stop
 
+!   MGGA functionals
+
+    case(3) 
+
+      call print_warning('!!! Individual energies NYI for MGGAs !!!')
+      stop
+
 !   Hybrid functionals
 
     case(4) 
@@ -59,12 +66,6 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns
 
       aC = 0.81d0
 
-!   Hartree-Fock calculation
-
-    case(666) 
-
-      Ec(:) = 0d0
-
   end select
  
 end subroutine unrestricted_correlation_individual_energy

src/eDFT/unrestricted_correlation_potential.f90

@@ -51,6 +51,12 @@ subroutine unrestricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nB
 
       call unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
 
+!   MGGA functionals
+
+    case(3)
+
+      call unrestricted_mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
+
 !   Hybrid functionals
 
     case(4)
@@ -64,12 +70,6 @@ subroutine unrestricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nB
 
       Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:))
 
-!   Hartree-Fock calculation
-
-    case(666)
-
-      Fc(:,:,:) = 0d0
-
   end select
 
 end subroutine unrestricted_correlation_potential

src/eDFT/unrestricted_exchange_derivative_discontinuity.f90

@@ -43,13 +43,18 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
 
       call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),&
                                                               rhow(:),Cx_choice,doNcentered,kappa,ExDD(:))
-
 !   GGA functionals
 
     case(2) 
 
       call unrestricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:))
 
+!   MGGA functionals
+
+    case(3) 
+
+      call unrestricted_mgga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:))
+
 !   Hybrid functionals
 
     case(4) 
@@ -57,12 +62,6 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC
       call print_warning('!!! exchange part of derivative discontinuity NYI for hybrids !!!')
       stop
 
-!   Hartree-Fock calculation
-
-    case(666) 
-
-      ExDD(:) = 0d0
-
   end select
  
 end subroutine unrestricted_exchange_derivative_discontinuity

src/eDFT/unrestricted_exchange_energy.f90

@@ -82,14 +82,6 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,a
          + cX*(ExHF  - ExLDA) &
          + aX*(ExGGA - ExLDA) 
 
-!   Hartree-Fock calculation
-
-    case(666) 
-
-      call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF)
-
-      Ex = ExHF
-
   end select
  
 end subroutine unrestricted_exchange_energy

src/eDFT/unrestricted_exchange_individual_energy.f90

@@ -33,6 +33,7 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
 
   double precision              :: ExLDA
   double precision              :: ExGGA
+  double precision              :: ExMGGA
   double precision              :: ExHF
 
 ! Output variables
@@ -64,6 +65,14 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
 
       Ex = ExGGA
 
+!   MGGA functionals
+
+    case(3) 
+
+      call unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,ExMGGA)
+
+      Ex = ExMGGA
+
 !   Hybrid functionals
 
     case(4) 
@@ -71,14 +80,6 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE
       call print_warning('!!! Individual energies NYI for Hybrids !!!')
       stop
 
-!   Hartree-Fock calculation
-
-    case(666) 
-
-      call unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,ExHF)
-
-      Ex = ExHF
-
   end select
  
 end subroutine unrestricted_exchange_individual_energy

src/eDFT/unrestricted_exchange_potential.f90

@@ -57,6 +57,12 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
 
       call unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
 
+!   MGGA functionals
+
+    case(3) 
+
+      call unrestricted_mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+
 !   Hybrid functionals
 
     case(4) 
@@ -74,14 +80,6 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w
               + cX*(FxHF(:,:)  - FxLDA(:,:)) &
               + aX*(FxGGA(:,:) - FxLDA(:,:))  
 
-!   Hartree-Fock calculation
-
-    case(666) 
-
-      call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF)
-
-      Fx(:,:) = FxHF(:,:)
-
   end select
 
 end subroutine unrestricted_exchange_potential

src/eDFT/unrestricted_mgga_correlation_potential.f90

@@ -0,0 +1,38 @@
+subroutine unrestricted_mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
+
+! Compute unrestricted MGGA correlation potential
+
+  implicit none
+  include 'parameters.h'
+
+! 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
+  double precision,intent(in)   :: AO(nBas,nGrid)
+  double precision,intent(in)   :: dAO(3,nBas,nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+  double precision,intent(in)   :: drho(3,nGrid,nspin)
+
+! Local variables
+
+! Output variables
+
+  double precision,intent(out)  :: Fc(nBas,nBas,nspin)
+
+! Select MGGA exchange functional
+
+  select case (DFA)
+
+    case default
+
+      call print_warning('!!! MGGA correlation potential not available !!!')
+      stop
+
+  end select
+
+end subroutine unrestricted_mgga_correlation_potential

src/eDFT/unrestricted_mgga_exchange_derivative_discontinuity.f90

@@ -0,0 +1,36 @@
+subroutine unrestricted_mgga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD)
+
+! Compute the exchange MGGA part of the derivative discontinuity
+
+  implicit none
+  include 'parameters.h'
+
+! 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)   :: rhow(nGrid)
+  double precision,intent(in)   :: drhow(ncart,nGrid)
+
+! Local variables
+
+
+! Output variables
+
+  double precision,intent(out)  :: ExDD(nEns)
+
+! Select exchange functional
+
+  select case (DFA)
+
+    case default
+
+      call print_warning('!!! MGGA exchange derivative discontinuity not available !!!')
+      stop
+
+  end select
+ 
+end subroutine unrestricted_mgga_exchange_derivative_discontinuity

src/eDFT/unrestricted_mgga_exchange_individual_energy.f90

@@ -0,0 +1,35 @@
+subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex)
+
+! Compute MGGA exchange energy for individual states
+
+  implicit none
+  include 'parameters.h'
+
+! 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)   :: rhow(nGrid)
+  double precision,intent(in)   :: drhow(ncart,nGrid)
+  double precision,intent(in)   :: rho(nGrid)
+  double precision,intent(in)   :: drho(ncart,nGrid)
+
+! Output variables
+
+  double precision              :: Ex
+
+! Select correlation functional
+
+  select case (DFA)
+
+    case default
+
+      call print_warning('!!! MGGA exchange individual energy not available !!!')
+      stop
+
+  end select
+
+end subroutine unrestricted_mgga_exchange_individual_energy

src/eDFT/unrestricted_mgga_exchange_potential.f90

@@ -0,0 +1,36 @@
+subroutine unrestricted_mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+
+! Select MGGA exchange functional for potential calculation
+
+  implicit none
+  include 'parameters.h'
+
+! 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
+  double precision,intent(in)   :: AO(nBas,nGrid)
+  double precision,intent(in)   :: dAO(3,nBas,nGrid)
+  double precision,intent(in)   :: rho(nGrid)
+  double precision,intent(in)   :: drho(3,nGrid)
+
+! Output variables
+
+  double precision,intent(out)  :: Fx(nBas,nBas)
+
+! Select MGGA exchange functional
+
+  select case (DFA)
+
+    case default
+
+      call print_warning('!!! MGGA exchange potential not available !!!')
+      stop
+
+  end select
+
+end subroutine unrestricted_mgga_exchange_potential