merging quack with eDFT

input/geminal

@@ -1 +0,0 @@
-1.0

src/IntPak/ReadGeminal.f90

@@ -1,25 +0,0 @@
-subroutine ReadGeminal(ExpS)
-
-! Read the geminal information
-
-  implicit none
-
-! Input variables
-  double precision,intent(out)  :: ExpS
-
-! Open file with geometry specification
-  open(unit=4,file='input/geminal')
-
-! Read exponent of Slater geminal
-  read(4,*) ExpS
-
-
-  write(*,'(A28)') '------------------'
-  write(*,'(A28,1X,F16.10)') 'Slater geminal exponent',ExpS
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Close file with geminal information
-  close(unit=4)
-
-end subroutine ReadGeminal

src/IntPak/ReadGeometry.f90

@@ -1,40 +0,0 @@
-subroutine ReadGeometry(NAtoms,ZNuc,XYZAtoms)
-
-! Read molecular geometry
-
-  implicit none
-
-! Input variables
-  integer,intent(in)            :: NAtoms
-  double precision,intent(out)  :: ZNuc(NAtoms),XYZAtoms(NAtoms,3)
-
-! Local variables
-  integer                       :: i
-
-! Open file with geometry specification
-  open(unit=1,file='input/molecule')
-
-! Read number of atoms
-  read(1,*)
-  read(1,*)
-  read(1,*)
-
-  do i=1,NAtoms
-    read(1,*) ZNuc(i),XYZAtoms(i,1),XYZAtoms(i,2),XYZAtoms(i,3)
-  enddo
-
-! Print geometry
-  write(*,'(A28)') 'Molecular geometry'
-  write(*,'(A28)') '------------------'
-  do i=1,NAtoms
-    write(*,'(A28,1X,I16)') 'Atom n. ',i
-    write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i)
-    write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',XYZAtoms(i,1),XYZAtoms(i,2),XYZAtoms(i,3)
-  enddo
-  write(*,'(A28)') '------------------'
-  write(*,*)
-
-! Close file with geometry specification
-  close(unit=1)
-
-end subroutine ReadGeometry

src/IntPak/read_geometry.f90

@@ -0,0 +1,68 @@
+subroutine read_geometry(nAt,ZNuc,rA,ENuc)
+
+! Read molecular geometry
+
+  implicit none
+
+  include 'parameters.h'
+
+! Ouput variables
+
+  integer,intent(in)            :: nAt
+
+! Local variables
+
+  integer                       :: i,j
+  double precision              :: RAB
+  character(len=2)              :: El
+  integer,external              :: element_number
+
+! Ouput variables
+
+  double precision,intent(out)  :: ZNuc(NAt),rA(nAt,ncart),ENuc
+
+! Open file with geometry specification
+
+  open(unit=1,file='input/molecule')
+
+! Read geometry
+
+  read(1,*) 
+  read(1,*) 
+  read(1,*) 
+
+  do i=1,nAt
+    read(1,*) El,rA(i,1),rA(i,2),rA(i,3)
+    ZNuc(i) = element_number(El)
+  enddo
+
+! Compute nuclear repulsion energy
+
+  ENuc = 0
+
+  do i=1,nAt-1
+    do j=i+1,nAt
+      RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2
+      ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB)
+    enddo
+  enddo
+
+! Close file with geometry specification
+  close(unit=1)
+
+! Print geometry
+  write(*,'(A28)') '------------------'
+  write(*,'(A28)') 'Molecular geometry'
+  write(*,'(A28)') '------------------'
+  do i=1,NAt
+    write(*,'(A28,1X,I16)') 'Atom n. ',i
+    write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i)
+    write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rA(i,j),j=1,ncart)
+  enddo
+  write(*,*)
+  write(*,'(A28)') '------------------'
+  write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc
+  write(*,'(A28)') '------------------'
+  write(*,*)
+
+end subroutine read_geometry

src/xcDFT/AO_values_grid.f90

@@ -0,0 +1,101 @@
+subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & 
+                          nGrid,root,AO,dAO)
+
+! Compute values of the AOs and their derivatives with respect to the cartesian coordinates
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nShell
+  double precision,intent(in)   :: CenterShell(maxShell,3)
+  integer,intent(in)            :: TotAngMomShell(maxShell)
+  integer,intent(in)            :: KShell(maxShell)
+  double precision,intent(in)   :: DShell(maxShell,maxK)
+  double precision,intent(in)   :: ExpShell(maxShell,maxK)
+  double precision,intent(in)   :: root(3,nGrid)
+  integer,intent(in)            :: nGrid
+
+! Local variables
+
+  integer                       :: atot,nShellFunction,a(3)
+  integer,allocatable           :: ShellFunction(:,:)
+  double precision              :: rASq,xA,yA,zA,NormCoeff,prim
+
+  integer                       :: iSh,iShF,iK,iG,iBas
+
+! Output variables
+
+  double precision,intent(out)  :: AO(nBas,nGrid)
+  double precision,intent(out)  :: dAO(3,nBas,nGrid)
+
+! Initialization
+
+  iBas = 0
+  AO(:,:)    = 0d0
+  dAO(:,:,:) = 0d0
+
+!------------------------------------------------------------------------
+! Loops over shells
+!------------------------------------------------------------------------
+  do iSh=1,nShell
+
+    atot = TotAngMomShell(iSh)
+    nShellFunction = (atot*atot + 3*atot + 2)/2
+    allocate(ShellFunction(1:nShellFunction,1:3))
+    call generate_shell(atot,nShellFunction,ShellFunction)
+
+    do iShF=1,nShellFunction
+
+      iBas = iBas + 1
+      a(:) = ShellFunction(iShF,:)
+
+      do iG=1,nGrid
+
+        xA = root(1,iG) - CenterShell(iSh,1)
+        yA = root(2,iG) - CenterShell(iSh,2)
+        zA = root(3,iG) - CenterShell(iSh,3)
+
+!       Calculate distance for exponential
+
+        rASq = xA**2 + yA**2 + zA**2
+
+!------------------------------------------------------------------------
+!         Loops over contraction degrees
+!-------------------------------------------------------------------------
+        do iK=1,KShell(iSh)
+          
+!         Calculate the exponential part
+
+          prim = DShell(iSh,iK)*NormCoeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq)        
+          AO(iBas,iG) = AO(iBas,iG) + prim
+
+          prim = -2d0*ExpShell(iSh,iK)*prim
+          dAO(:,iBas,iG) = dAO(:,iBas,iG) + prim
+        
+        enddo
+
+        dAO(1,iBas,iG) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(1,iBas,iG)
+        if(a(1) > 0) dAO(1,iBas,iG) = dAO(1,iBas,iG) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
+        
+        dAO(2,iBas,iG) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(2,iBas,iG)
+        if(a(2) > 0) dAO(2,iBas,iG) = dAO(2,iBas,iG) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iBas,iG)
+        
+        dAO(3,iBas,iG) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(3,iBas,iG)
+        if(a(3) > 0) dAO(3,iBas,iG) = dAO(3,iBas,iG) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iBas,iG)
+
+!       Calculate polynmial part
+
+        AO(iBas,iG) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
+        
+      enddo
+      
+    enddo
+    deallocate(ShellFunction)
+  enddo
+!------------------------------------------------------------------------
+! End loops over shells 
+!------------------------------------------------------------------------
+
+end subroutine AO_values_grid

src/xcDFT/B88_gga_exchange_energy.f90

@@ -0,0 +1,53 @@
+subroutine B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
+
+! Compute Becke's 96 GGA exchange energy
+
+  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)   :: rho(nGrid)
+  double precision,intent(in)   :: drho(3,nGrid)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: alpha,beta,gamma
+  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,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
+
+end subroutine B88_gga_exchange_energy

src/xcDFT/B88_gga_exchange_potential.f90

@@ -0,0 +1,68 @@
+subroutine B88_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+
+! Compute Becke's GGA exchange 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)
+  double precision,intent(in)   :: drho(3,nGrid)
+
+! Local variables
+
+  integer                       :: mu,nu,iG
+  double precision              :: alpha,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,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 B88_gga_exchange_potential

src/xcDFT/C16_lda_correlation_energy.f90

@@ -0,0 +1,93 @@
+subroutine C16_lda_correlation_energy(nGrid,weight,rho,Ec)
+
+! Compute Chachiyo's LDA correlation energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,rs
+  double precision              :: a_p,b_p,ec_p
+  double precision              :: a_f,b_f,ec_f
+  double precision              :: z,fz,ec_z
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Coefficients for Chachiyo's LDA correlation
+
+  a_p = (log(2d0) - 1d0)/(2d0*pi**2)
+  b_p = 20.4562557d0
+
+  a_f = (log(2d0) - 1d0)/(4d0*pi**2)
+  b_f = 27.4203609d0
+
+! Compute LDA correlation energy
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+!   Spin-up and spin-down densities
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+
+!   Total density
+
+    r = ra + rb
+
+!   Spin-up part contribution
+
+    if(ra > threshold) then
+
+      rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
+      ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
+
+      Ec(1) = Ec(1) + weight(iG)*ec_f*ra
+
+    endif
+
+!   Opposite-spin contribution
+
+    if(r > threshold) then
+
+      rs  = (4d0*pi*r/3d0)**(-1d0/3d0)
+      ec_p = a_p*log(1d0 + b_p/rs + b_p/rs**2)
+      ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
+
+      z = (ra-rb)/r
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+
+      ec_z = ec_p + (ec_f - ec_p)*fz
+
+      Ec(2) = Ec(2) + weight(iG)*ec_z*r
+
+      endif
+
+!   Spin-down contribution
+
+    if(rb > threshold) then
+
+      rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
+      ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
+
+      Ec(3) = Ec(3) + weight(iG)*ec_f*rb
+
+    endif
+
+  enddo
+
+  Ec(2) = Ec(2) - Ec(1) - Ec(3)
+
+end subroutine C16_lda_correlation_energy

src/xcDFT/C16_lda_correlation_potential.f90

@@ -0,0 +1,131 @@
+subroutine C16_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
+
+! Compute LDA correlation 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,nspin)
+
+! Local variables
+
+  integer                       :: mu,nu,iG
+  double precision              :: ra,rb,r,rs
+  double precision              :: a_p,b_p,ec_p,decdrs_p,decdra_p,decdrb_p
+  double precision              :: a_f,b_f,ec_f,decdrs_f,decdra_f,decdrb_f
+  double precision              :: ec_z,decdra_z,decdrb_z
+  double precision              :: z,dzdra,dzdrb,fz,dfzdz,dfzdra,dfzdrb
+  double precision              :: drsdra,drsdrb,dFcdra,dFcdrb
+
+! Output variables
+
+  double precision,intent(out)  :: Fc(nBas,nBas,nspin)
+
+! Coefficients for Chachiyo's LDA correlation
+
+  a_p = (log(2d0) - 1d0)/(2d0*pi**2)
+  b_p = 20.4562557d0
+
+  a_f = (log(2d0) - 1d0)/(4d0*pi**2)
+  b_f = 27.4203609d0
+
+! Compute LDA correlation matrix in the AO basis
+
+  Fc(:,:,:) = 0d0
+
+  do mu=1,nBas
+    do nu=1,nBas
+      do iG=1,nGrid
+
+!       Spin-up and spin-down densities
+
+        ra = max(0d0,rho(iG,1))
+        rb = max(0d0,rho(iG,2))
+      
+!       Total density
+
+        r = ra + rb
+
+!       Spin-up part contribution
+
+        if(ra > threshold) then
+
+          rs  = (4d0*pi*r/3d0)**(-1d0/3d0)
+
+          ec_p = a_p*log(1d0 + b_p/rs + b_p/rs**2)
+          ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
+
+          z = (ra-rb)/r
+
+          fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+          fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+
+          ec_z = ec_p + (ec_f - ec_p)*fz
+
+          dzdra  = (1d0 - z)/r
+          dfzdz  = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+          dfzdra = dzdra*dfzdz
+
+          drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+
+          decdrs_p = - a_p/rs**2*(b_p + 2d0*b_p/rs)/(1d0 + b_p/rs + b_p/rs**2)
+          decdrs_f = - a_f/rs**2*(b_f + 2d0*b_f/rs)/(1d0 + b_f/rs + b_f/rs**2)
+
+          decdra_p = drsdra*decdrs_p
+          decdra_f = drsdra*decdrs_f
+
+          decdra_z = decdra_p + (decdra_f - decdra_p)*fz + (ec_f - ec_p)*dfzdra
+
+          dFcdra = decdra_z*r + ec_z
+
+          Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdra
+       
+        endif
+
+!       Spin-down part contribution
+
+        if(rb > threshold) then
+
+          rs  = (4d0*pi*r/3d0)**(-1d0/3d0)
+
+          ec_p = a_p*log(1d0 + b_p/rs + b_p/rs**2)
+          ec_f = a_f*log(1d0 + b_f/rs + b_f/rs**2)
+
+          z = (ra-rb)/r
+
+          fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+          fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+
+          ec_z = ec_p + (ec_f - ec_p)*fz
+
+          dzdrb  = - (1d0 + z)/r
+          dfzdz  = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+          dfzdrb = dzdrb*dfzdz
+
+          drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+
+          decdrs_p = - a_p/rs**2*(b_p + 2d0*b_p/rs)/(1d0 + b_p/rs + b_p/rs**2)
+          decdrs_f = - a_f/rs**2*(b_f + 2d0*b_f/rs)/(1d0 + b_f/rs + b_f/rs**2)
+
+          decdrb_p = drsdrb*decdrs_p
+          decdrb_f = drsdrb*decdrs_f
+
+          decdrb_z = decdrb_p + (decdrb_f - decdrb_p)*fz + (ec_f - ec_p)*dfzdrb
+
+          dFcdrb = decdrb_z*r + ec_z
+
+          Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdrb
+
+        endif
+
+      enddo
+    enddo
+  enddo
+
+end subroutine C16_lda_correlation_potential

src/xcDFT/DIIS_extrapolation.f90

@@ -0,0 +1,54 @@
+subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout)
+
+! Perform DIIS extrapolation
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: n_err,n_e
+  double precision,intent(in)   :: error_in(n_err),error(n_err,n_diis),e(n_e,n_diis)
+
+! Local variables
+
+  double precision,allocatable  :: A(:,:),b(:),w(:)
+
+! Output variables
+
+  double precision,intent(out)  :: rcond
+  integer,intent(inout)         :: n_diis
+  double precision,intent(inout):: e_inout(n_e)
+
+! Memory allocaiton
+
+  allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1))
+
+! Update DIIS "history"
+
+  call prepend(n_err,n_diis,error,error_in)
+  call prepend(n_e,n_diis,e,e_inout)
+
+!  Build A matrix
+
+  A(1:n_diis,1:n_diis) = matmul(transpose(error),error)
+
+  A(1:n_diis,n_diis+1) = -1d0
+  A(n_diis+1,1:n_diis) = -1d0
+  A(n_diis+1,n_diis+1) = +0d0
+
+! Build x matrix
+
+  b(1:n_diis) = +0d0
+  b(n_diis+1) = -1d0
+
+! Solve linear system
+
+  call linear_solve(n_diis+1,A,b,w,rcond)
+
+! Extrapolate
+
+    e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))
+
+end subroutine DIIS_extrapolation

src/xcDFT/G96_gga_exchange_energy.f90

@@ -0,0 +1,52 @@
+subroutine G96_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
+
+! Compute Gill's 96 GGA exchange energy
+
+  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)   :: rho(nGrid)
+  double precision,intent(in)   :: drho(3,nGrid)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: alpha,beta
+  double precision              :: r,g
+
+! Output variables
+
+  double precision              :: Ex
+
+! Coefficients for G96 GGA exchange functional
+
+  alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
+  beta  = 1d0/137d0
+
+! Compute GGA exchange energy
+
+  Ex = 0d0
+
+  do iG=1,nGrid
+
+    r = max(0d0,rho(iG))
+
+    if(r > threshold) then 
+ 
+      g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
+
+      Ex = Ex + weight(iG)*r**(4d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2)
+
+    end if
+
+  end do
+
+end subroutine G96_gga_exchange_energy

src/xcDFT/G96_gga_exchange_potential.f90

@@ -0,0 +1,69 @@
+subroutine G96_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+
+! Compute Gill's GGA exchange poential
+
+  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)
+
+! Local variables
+
+  integer                       :: mu,nu,iG
+  double precision              :: alpha,beta
+  double precision              :: r,g,vAO,gAO
+
+! Output variables
+
+  double precision,intent(out)  :: Fx(nBas,nBas)
+
+! Coefficients for G96 GGA exchange functional
+
+  alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
+  beta  = +1d0/137d0
+  beta  = 0d0
+
+! 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,rho(iG))
+        g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
+ 
+        if(r > threshold) then
+
+          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)
+
+        endif
+
+      enddo
+    enddo
+  enddo
+
+end subroutine G96_gga_exchange_potential

src/xcDFT/Kohn_Sham.f90

@@ -0,0 +1,362 @@
+subroutine Kohn_Sham(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,Ew)
+
+! Perform unrestricted Kohn-Sham calculation for ensembles
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: x_rung,c_rung
+  character(len=12),intent(in)  :: x_DFA,c_DFA
+  integer,intent(in)            :: nEns
+  double precision,intent(in)   :: wEns(nEns)
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  integer,intent(in)            :: maxSCF,max_diis,guess_type
+  double precision,intent(in)   :: thresh
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: AO(nBas,nGrid)
+  double precision,intent(in)   :: dAO(ncart,nBas,nGrid)
+
+  integer,intent(in)            :: nO(nspin),nV(nspin)
+  double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: T(nBas,nBas)
+  double precision,intent(in)   :: V(nBas,nBas)
+  double precision,intent(in)   :: Hc(nBas,nBas) 
+  double precision,intent(in)   :: X(nBas,nBas) 
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ENuc
+
+! Local variables
+
+  integer                       :: xc_rung
+  integer                       :: nSCF,nBasSq
+  integer                       :: n_diis
+  double precision              :: conv
+  double precision              :: rcond(nspin)
+  double precision              :: ET(nspin)
+  double precision              :: EV(nspin)
+  double precision              :: EJ(nsp)
+  double precision              :: Ex(nspin)
+  double precision              :: Ec(nsp)
+  double precision              :: Ew
+
+  double precision,allocatable  :: eps(:,:)
+  double precision,allocatable  :: c(:,:,:)
+  double precision,allocatable  :: cp(:,:,:)
+  double precision,allocatable  :: J(:,:,:)
+  double precision,allocatable  :: F(:,:,:)
+  double precision,allocatable  :: Fp(:,:,:)
+  double precision,allocatable  :: Fx(:,:,:)
+  double precision,allocatable  :: FxHF(:,:,:)
+  double precision,allocatable  :: Fc(:,:,:)
+  double precision,allocatable  :: err(:,:,:)
+  double precision,allocatable  :: err_diis(:,:,:)
+  double precision,allocatable  :: F_diis(:,:,:)
+  double precision,external     :: trace_matrix
+  double precision,external     :: electron_number
+
+  double precision,allocatable  :: Pw(:,:,:)
+  double precision,allocatable  :: rhow(:,:)
+  double precision,allocatable  :: drhow(:,:,:)
+  double precision              :: nEl(nspin)
+
+  double precision,allocatable  :: P(:,:,:,:)
+  double precision,allocatable  :: rho(:,:,:)
+  double precision,allocatable  :: drho(:,:,:,:)
+
+  double precision              :: E(nEns)
+  double precision              :: Om(nEns)
+
+  integer                       :: ispin,iEns
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'************************************************'
+  write(*,*)'*    Unrestricted Kohn-Sham calculation        *'
+  write(*,*)'*           *** for ensembles ***              *'
+  write(*,*)'************************************************'
+  write(*,*)
+
+! Useful stuff
+
+  nBasSq = nBas*nBas
+
+!------------------------------------------------------------------------
+! Rung of Jacob's ladder
+!------------------------------------------------------------------------
+
+! Select rung for exchange 
+
+  write(*,*)
+  write(*,*) '*******************************************************************'
+  write(*,*) '*                        Exchange rung                            *'
+  write(*,*) '*******************************************************************'
+
+  call select_rung(x_rung,x_DFA)
+
+! Select rung for correlation
+
+  write(*,*)
+  write(*,*) '*******************************************************************'
+  write(*,*) '*                       Correlation rung                          *'
+  write(*,*) '*******************************************************************'
+
+  call select_rung(c_rung,c_DFA)
+
+! Overall rung
+
+  xc_rung = max(x_rung,c_rung)
+
+! Memory allocation
+
+  allocate(eps(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),        &
+           J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin),     & 
+           Fx(nBas,nBas,nspin),FxHF(nBas,nBas,nspin),Fc(nBas,nBas,nspin),err(nBas,nBas,nspin), &
+           Pw(nBas,nBas,nspin),rhow(nGrid,nspin),drhow(ncart,nGrid,nspin), &
+           err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin),  &
+           P(nBas,nBas,nspin,nEns),rho(nGrid,nspin,nEns),drho(ncart,nGrid,nspin,nEns))
+
+! Guess coefficients and eigenvalues
+
+  if(guess_type == 1) then
+
+    do ispin=1,nspin
+      F(:,:,ispin) = Hc(:,:)
+    end do
+
+  else if(guess_type == 2) then
+
+    do ispin=1,nspin
+      call random_number(F(:,:,ispin))
+    end do
+
+  end if
+
+! Initialization
+
+  nSCF = 0
+  conv = 1d0
+
+  nEl(:) = 0d0
+  Ex(:)  = 0d0
+  Ec(:)  = 0d0
+
+  Fx(:,:,:)   = 0d0
+  FxHF(:,:,:) = 0d0
+  Fc(:,:,:)   = 0d0
+
+  n_diis          = 0
+  F_diis(:,:,:)   = 0d0
+  err_diis(:,:,:) = 0d0
+
+!------------------------------------------------------------------------
+! Main SCF loop
+!------------------------------------------------------------------------
+
+  write(*,*)
+  write(*,*)'------------------------------------------------------------------------------------------'
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & 
+            '|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|'
+  write(*,*)'------------------------------------------------------------------------------------------'
+  
+  do while(conv > thresh .and. nSCF < maxSCF)
+
+!   Increment 
+
+    nSCF = nSCF + 1
+
+!  Transform Fock matrix in orthogonal basis
+
+    do ispin=1,nspin
+      Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
+    end do
+
+!  Diagonalize Fock matrix to get eigenvectors and eigenvalues
+
+    cp(:,:,:) = Fp(:,:,:)
+    do ispin=1,nspin
+      call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin))
+    end do
+    
+!   Back-transform eigenvectors in non-orthogonal basis
+
+    do ispin=1,nspin
+      c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
+    end do
+
+!------------------------------------------------------------------------
+!   Compute density matrix 
+!------------------------------------------------------------------------
+
+    call density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:))
+ 
+!   Weight-dependent density matrix
+    
+    Pw(:,:,:) = 0d0
+    do iEns=1,nEns
+      Pw(:,:,:) = Pw(:,:,:) + wEns(iEns)*P(:,:,:,iEns)
+    end do
+
+!------------------------------------------------------------------------
+!   Compute one-electron density and its gradient if necessary
+!------------------------------------------------------------------------
+
+    do ispin=1,nspin
+      do iEns=1,nEns
+        call density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),rho(:,ispin,iEns))
+      end do
+    end do
+
+!   Weight-dependent one-electron density 
+
+    rhow(:,:) = 0d0
+    do iEns=1,nEns
+      rhow(:,:) = rhow(:,:) + wEns(iEns)*rho(:,:,iEns) 
+    end do
+
+    if(xc_rung > 1 .and. xc_rung /= 666) then 
+
+!     Ground state density 
+
+      do ispin=1,nspin
+        do iEns=1,nEns
+          call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns))
+        end do
+      end do
+
+!     Weight-dependent one-electron density 
+
+      drhow(:,:,:) = 0d0
+      do iEns=1,nEns
+        drhow(:,:,:) = drhow(:,:,:) + wEns(iEns)*drho(:,:,:,iEns)
+      end do
+
+    end if
+
+!   Build Coulomb repulsion
+
+    do ispin=1,nspin
+      call hartree_coulomb(nBas,Pw(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin))
+    end do
+
+!   Compute exchange potential
+
+    do ispin=1,nspin
+      call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,Pw(:,:,ispin),ERI(:,:,:,:), &
+                              AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin),Fx(:,:,ispin),FxHF(:,:,ispin))
+    end do
+
+!   Compute correlation potential
+
+    call correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rhow(:,:),drhow(:,:,:),Fc(:,:,:))
+
+!   Build Fock operator
+    do ispin=1,nspin
+      F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + Fx(:,:,ispin) + Fc(:,:,ispin)
+    end do
+
+!   Check convergence 
+
+    do ispin=1,nspin
+      err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin))
+    end do
+
+    conv = maxval(abs(err(:,:,:)))
+    
+!   DIIS extrapolation
+
+    n_diis = min(n_diis+1,max_diis)
+    do ispin=1,nspin
+      call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))
+    end do
+
+!   Reset DIIS if required
+
+    if(minval(rcond(:)) < 1d-15) n_diis = 0
+
+!------------------------------------------------------------------------
+!   Compute KS energy
+!------------------------------------------------------------------------
+
+!  Kinetic energy
+
+    do ispin=1,nspin
+      ET(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),T(:,:)))
+    end do
+
+!  Potential energy
+
+    do ispin=1,nspin
+      EV(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),V(:,:)))
+    end do
+
+!  Coulomb energy
+
+    EJ(1) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1)))
+    EJ(2) = trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2)))
+    EJ(3) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2)))
+
+!   Exchange energy
+
+    do ispin=1,nspin
+      call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas, &
+                           Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin))
+    end do
+
+!   Correlation energy
+
+    call correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),Ec)
+
+!   Total energy
+
+    Ew = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:))
+
+!   Check the grid accuracy by computing the number of electrons 
+
+    do ispin=1,nspin
+      nEl(ispin) = electron_number(nGrid,weight(:),rhow(:,ispin))
+    end do
+
+!   Dump results
+
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & 
+      '|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',conv,'|',sum(nEl(:)),'|'
+ 
+  end do
+  write(*,*)'------------------------------------------------------------------------------------------'
+!------------------------------------------------------------------------
+! End of SCF loop
+!------------------------------------------------------------------------
+
+! Did it actually converge?
+
+  if(nSCF == maxSCF) then
+
+    write(*,*)
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)'                 Convergence failed                 '
+    write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+    write(*,*)
+
+    stop
+
+  end if
+
+! Compute final KS energy
+
+  call print_KS(nBas,nO(:),eps(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Ew)
+
+!------------------------------------------------------------------------
+! Compute individual energies from ensemble energy
+!------------------------------------------------------------------------
+
+  call individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas,     &
+                         AO(:,:),dAO(:,:,:),nO(:),nV(:),T(:,:),V(:,:),ERI(:,:,:,:),ENuc,  & 
+                         Pw(:,:,:),rhow(:,:),drhow(:,:,:),J(:,:,:),Fx(:,:,:),FxHF(:,:,:), & 
+                         Fc(:,:,:),P(:,:,:,:),rho(:,:,:),drho(:,:,:,:),E(:),Om(:))
+
+end subroutine Kohn_Sham 

src/xcDFT/LF19_lda_correlation_Levy_Zahariev_shift.f90

@@ -0,0 +1,82 @@
+subroutine LF19_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight,rho,EcLZ)
+
+! Compute Loos-Fromager's LDA contribution to Levy-Zahariev shift
+
+  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)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  logical                       :: LDA_centered = .false.
+  integer                       :: iEns
+  double precision              :: EcLZLDA(nsp)
+  double precision,allocatable  :: aLF(:,:)
+  double precision,allocatable  :: EcLZeLDA(:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: EcLZ(nsp)
+
+
+! Allocation
+
+  allocate(aLF(3,nEns),EcLZeLDA(nsp,nEns))
+
+! Parameters for weight-dependent LDA correlation functional
+
+  aLF(1,1) = -0.0238184d0
+  aLF(2,1) = +0.00575719d0
+  aLF(3,1) = +0.0830576d0
+
+  aLF(1,2) = -0.0282814d0
+  aLF(2,2) = +0.00340758d0
+  aLF(3,2) = +0.0663967d0
+
+  aLF(1,3) = -0.0144633d0
+  aLF(2,3) = -0.0504501d0
+  aLF(3,3) = +0.0331287d0
+
+! Compute correlation energy for ground, singly-excited and doubly-excited states
+
+  do iEns=1,nEns
+
+    call elda_correlation_Levy_Zahariev_shift(nEns,aLF(:,iEns),nGrid,weight(:),rho(:,:),EcLZeLDA(:,iEns))
+
+  end do
+
+! LDA-centered functional
+
+  EcLZLDA(:) = 0d0
+
+  if(LDA_centered) then
+
+    call VWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:,:),EcLZLDA(:))
+
+    do iEns=1,nEns
+ 
+      EcLZeLDA(:,iEns) = EcLZeLDA(:,iEns) + EcLZLDA(:)- EcLZeLDA(:,1)
+ 
+    end do
+
+  end if
+
+! Weight-denpendent functional for ensembles
+
+  EcLZ(:) = 0d0
+
+  do iEns=1,nEns
+
+    EcLZ(:) = EcLZ(:) + wEns(iEns)*EcLZeLDA(:,iEns)
+
+  enddo
+
+end subroutine LF19_lda_correlation_Levy_Zahariev_shift

src/xcDFT/LF19_lda_correlation_derivative_discontinuity.f90

@@ -0,0 +1,63 @@
+subroutine LF19_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,Ec)
+
+! Compute eLDA correlation part of the derivative discontinuity
+
+  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)
+  double precision,intent(in)   :: rhow(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iEns,jEns
+  double precision,allocatable  :: aLF(:,:)
+  double precision              :: dEc(nsp,nEns)
+  double precision,external     :: Kronecker_delta
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp,nEns)
+
+! Allocation
+
+  allocate(aLF(3,nEns))
+
+! Parameters for weight-dependent LDA correlation functional
+
+  aLF(1,1) = -0.0238184d0
+  aLF(2,1) = +0.00575719d0
+  aLF(3,1) = +0.0830576d0
+
+  aLF(1,2) = -0.0282814d0
+  aLF(2,2) = +0.00340758d0
+  aLF(3,2) = +0.0663967d0
+
+  aLF(1,3) = -0.0144633d0
+  aLF(2,3) = -0.0504501d0
+  aLF(3,3) = +0.0331287d0
+
+! Compute correlation energy for ground, singly-excited and doubly-excited states
+
+  do iEns=1,nEns
+
+    call elda_correlation_energy(nEns,aLF(:,iEns),nGrid,weight(:),rhow(:,:),dEc(:,iEns))
+
+  end do
+
+  Ec(:,:) = 0d0
+
+  do iEns=1,nEns
+    do jEns=1,nEns
+
+      Ec(:,iEns) = Ec(:,iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*(dEc(:,jEns) - dEc(:,1))
+
+    end do
+  end do
+
+end subroutine LF19_lda_correlation_derivative_discontinuity

src/xcDFT/LF19_lda_correlation_energy.f90

@@ -0,0 +1,84 @@
+subroutine LF19_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec)
+
+! Compute eLDA correlation energy 
+
+  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)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  logical                       :: LDA_centered = .false.
+  integer                       :: iEns,isp
+  double precision              :: EcLDA(nsp)
+  double precision,allocatable  :: aLF(:,:)
+  double precision,allocatable  :: EceLDA(:,:)
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Allocation
+
+  allocate(aLF(3,nEns),EceLDA(nsp,nEns))
+
+! Parameters for weight-dependent LDA correlation functional
+
+  aLF(1,1) = -0.0238184d0
+  aLF(2,1) = +0.00575719d0
+  aLF(3,1) = +0.0830576d0
+
+  aLF(1,2) = -0.0282814d0
+  aLF(2,2) = +0.00340758d0
+  aLF(3,2) = +0.0663967d0
+
+  aLF(1,3) = -0.0144633d0
+  aLF(2,3) = -0.0504501d0
+  aLF(3,3) = +0.0331287d0
+
+! Compute correlation energy for ground, singly-excited and doubly-excited states
+
+  do iEns=1,nEns
+
+    call elda_correlation_energy(nEns,aLF(:,iEns),nGrid,weight(:),rho(:,:),EceLDA(:,iEns))
+
+  end do
+
+! LDA-centered functional
+
+  EcLDA(:) = 0d0
+
+  if(LDA_centered) then 
+
+    call VWN5_lda_correlation_energy(nGrid,weight(:),rho(:,:),EcLDA(:))
+
+    do iEns=1,nEns
+      do isp=1,nsp
+
+        EceLDA(isp,iEns) = EceLDA(isp,iEns) + EcLDA(isp) - EceLDA(isp,1)
+
+      end do
+    end do
+
+  end if
+
+! Weight-denpendent functional for ensembles
+
+  Ec(:) = 0d0
+
+  do iEns=1,nEns
+    do isp=1,nsp
+
+      Ec(isp) = Ec(isp) + wEns(iEns)*EceLDA(isp,iEns)
+
+    end do
+  end do
+
+end subroutine LF19_lda_correlation_energy

src/xcDFT/LF19_lda_correlation_individual_energy.f90

@@ -0,0 +1,81 @@
+subroutine LF19_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ec)
+
+! Compute eLDA correlation energy 
+
+  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)
+  double precision,intent(in)   :: rhow(nGrid,nspin)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  logical                       :: LDA_centered = .false.
+  integer                       :: iEns,isp
+  double precision              :: EcLDA(nsp)
+  double precision,allocatable  :: aLF(:,:)
+  double precision,allocatable  :: EceLDA(:,:)
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Allocation
+
+  allocate(aLF(3,nEns),EceLDA(nsp,nEns))
+
+! Parameters for weight-dependent LDA correlation functional
+
+  aLF(1,1) = -0.0238184d0
+  aLF(2,1) = +0.00575719d0
+  aLF(3,1) = +0.0830576d0
+
+  aLF(1,2) = -0.0282814d0
+  aLF(2,2) = +0.00340758d0
+  aLF(3,2) = +0.0663967d0
+
+  aLF(1,3) = -0.0144633d0
+  aLF(2,3) = -0.0504501d0
+  aLF(3,3) = +0.0331287d0
+
+! Compute correlation energy for ground, singly-excited and doubly-excited states
+
+  do iEns=1,nEns
+
+    call elda_correlation_individual_energy(nEns,aLF(:,iEns),nGrid,weight(:),rhow(:,:),rho(:,:),EceLDA(:,iEns))
+
+  end do
+
+! LDA-centered functional
+
+  EcLDA(:) = 0d0
+
+  if(LDA_centered) then
+
+    call W38_lda_correlation_individual_energy(nGrid,weight(:),rhow(:,:),rho(:,:),EcLDA(:))
+
+    do iEns=1,nEns
+      do isp=1,nsp
+        EceLDA(isp,iEns) = EceLDA(isp,iEns) + EcLDA(isp) - EceLDA(isp,1)
+      end do
+    end do
+
+  end if
+
+! Weight-denpendent functional for ensembles
+
+  Ec(:) = 0d0
+
+  do iEns=1,nEns
+    do isp=1,nsp
+      Ec(isp) = Ec(isp) + wEns(iEns)*EceLDA(isp,iEns)
+    enddo
+  enddo
+
+end subroutine LF19_lda_correlation_individual_energy

src/xcDFT/LF19_lda_correlation_potential.f90

@@ -0,0 +1,82 @@
+subroutine LF19_lda_correlation_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fc)
+
+! Compute Loos-Fromager weight-dependent LDA correlation 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,nspin)
+
+! Local variables
+
+  logical                       :: LDA_centered = .false.
+  integer                       :: iEns
+  double precision,allocatable  :: aLF(:,:)
+  double precision,allocatable  :: FcLDA(:,:,:)
+  double precision,allocatable  :: FceLDA(:,:,:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: Fc(nBas,nBas,nspin)
+
+! Allocation
+
+  allocate(aLF(3,nEns),FcLDA(nBas,nBas,nspin),FceLDA(nBas,nBas,nspin,nEns))
+
+! Parameters for weight-dependent LDA correlation functional
+
+  aLF(1,1) = -0.0238184d0
+  aLF(2,1) = +0.00575719d0
+  aLF(3,1) = +0.0830576d0
+
+  aLF(1,2) = -0.0282814d0
+  aLF(2,2) = +0.00340758d0
+  aLF(3,2) = +0.0663967d0
+
+  aLF(1,3) = -0.0144633d0
+  aLF(2,3) = -0.0504501d0
+  aLF(3,3) = +0.0331287d0
+
+! Compute correlation energy for ground, singly-excited and doubly-excited states
+
+  do iEns=1,nEns
+
+    call elda_correlation_potential(nEns,aLF(:,iEns),nGrid,weight,nBas,AO,rho,FceLDA(:,:,:,iEns))
+
+  end do
+
+! LDA-centered functional
+
+  FcLDA(:,:,:) = 0d0
+
+  if(LDA_centered) then 
+
+    call VWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,FcLDA)
+
+    do iEns=1,nEns
+
+      FceLDA(:,:,:,iEns) = FceLDA(:,:,:,iEns) + FcLDA(:,:,:) - FceLDA(:,:,:,1)
+
+    end do
+
+  end if
+
+! Weight-denpendent functional for ensembles
+
+  Fc(:,:,:) = 0d0
+
+  do iEns=1,nEns
+
+    Fc(:,:,:) = Fc(:,:,:) + wEns(iEns)*FceLDA(:,:,:,iEns)
+
+  enddo
+
+end subroutine LF19_lda_correlation_potential

src/xcDFT/Makefile

@@ -0,0 +1,35 @@
+IDIR =../../include
+BDIR =../../bin
+ODIR = obj
+SDIR =.
+FC = gfortran -I$(IDIR) 
+ifeq ($(DEBUG),1)
+FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
+else
+FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O2 
+endif
+
+LIBS = ~/Dropbox/quack/lib/*.a
+#LIBS = -lblas -llapack
+
+SRCF90 = $(wildcard *.f90)
+
+SRC = $(wildcard *.f)
+
+OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC))
+
+$(ODIR)/%.o: %.f90
+	$(FC) -c -o $@ $< $(FFLAGS) 
+
+$(ODIR)/%.o: %.f
+	$(FC) -c -o $@ $< $(FFLAGS) 
+
+$(BDIR)/xcDFT: $(OBJ)
+	$(FC) -o $@ $^ $(FFLAGS) $(LIBS)
+
+debug: 
+	DEBUG=1 make $(BDIR)/xcDFT
+#DEBUG=1 make clean $(BDIR)/xcDFT
+
+clean:
+	rm -f $(ODIR)/*.o $(BDIR)/xcDFT $(BDIR)/debug

src/xcDFT/S51_lda_exchange_energy.f90

@@ -0,0 +1,42 @@
+subroutine S51_lda_exchange_energy(nGrid,weight,rho,Ex)
+
+! Compute 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/2d0)*(3d0/(4d0*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 S51_lda_exchange_energy

src/xcDFT/S51_lda_exchange_potential.f90

@@ -0,0 +1,50 @@
+subroutine S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
+
+! Compute 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              :: alpha 
+  double precision              :: r,vAO
+
+! Output variables
+
+  double precision,intent(out)  :: Fx(nBas,nBas)
+
+! Cx coefficient for Slater LDA exchange
+
+  alpha = -(3d0/2d0)*(3d0/(4d0*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*alpha*r**(1d0/3d0)
+
+        endif
+
+      enddo
+    enddo
+  enddo
+
+end subroutine S51_lda_exchange_potential

src/xcDFT/VWN5_lda_correlation_Levy_Zahariev_shift.f90

@@ -0,0 +1,199 @@
+subroutine VWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight,rho,EcLZ)
+
+! Compute Wigner's LDA contribution to Levy-Zahariev shift
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,rs,x,z
+  double precision              :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
+  double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
+  double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
+  double precision              :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
+  double precision              :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
+  double precision              :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
+
+  double precision              :: ec_z,ec_p,ec_f,ec_a
+  double precision              :: fz,d2fz
+
+! Output variables
+
+  double precision,intent(out)  :: EcLZ(nsp)
+
+! Parameters of the functional
+
+  a_p  = +0.0621814D0/2D0
+  x0_p = -0.10498d0
+  b_p  = +3.72744d0
+  c_p  = +12.9352d0
+
+  a_f  = +0.0621814D0/4D0
+  x0_f = -0.325d0
+  b_f  = +7.06042d0
+  c_f  = +18.0578d0
+
+  a_a  = -1d0/(6d0*pi**2)
+  x0_a = -0.0047584D0
+  b_a  = +1.13107d0
+  c_a  = +13.0045d0
+
+! Initialization
+
+  EcLZ(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+
+!   spin-up contribution
+   
+    if(ra > threshold) then
+   
+      r = ra + rb
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      z = (ra - rb)/r
+      x = sqrt(rs)
+   
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+   
+      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+   
+      x_p = x*x + b_p*x + c_p
+      x_f = x*x + b_f*x + c_f
+      x_a = x*x + b_a*x + c_a
+   
+      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+   
+      q_p = sqrt(4d0*c_p - b_p*b_p)
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+      q_a = sqrt(4d0*c_a - b_a*b_a)
+   
+      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+   
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+   
+      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+   
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+
+      dzdra = (1d0 - z)/r
+      dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+      dfzdra = dzdra*dfzdz
+
+      drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_p = 2d0*x + b_p
+      dxdx_f = 2d0*x + b_f
+      dxdx_a = 2d0*x + b_a
+
+      decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                      - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                      - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+      decdra_p = drsdra*dxdrs*decdx_p
+      decdra_f = drsdra*dxdrs*decdx_f
+      decdra_a = drsdra*dxdrs*decdx_a
+
+      decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+             + (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
+
+      EcLZ(2) = EcLZ(2) + weight(iG)*decdra*r*r
+   
+    end if
+
+!   spin-down contribution
+   
+    if(rb > threshold) then
+   
+      r = ra + rb
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      z = (ra - rb)/r
+      x = sqrt(rs)
+   
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+   
+      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+   
+      x_p = x*x + b_p*x + c_p
+      x_f = x*x + b_f*x + c_f
+      x_a = x*x + b_a*x + c_a
+
+      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+   
+      q_p = sqrt(4d0*c_p - b_p*b_p)
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+      q_a = sqrt(4d0*c_a - b_a*b_a)
+   
+      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+   
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+   
+      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+   
+      dzdrb = -(1d0 + z)/r
+      dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+      dfzdrb = dzdrb*dfzdz
+
+      drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_p = 2d0*x + b_p
+      dxdx_f = 2d0*x + b_f
+      dxdx_a = 2d0*x + b_a
+
+      decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                      - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                      - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+      decdrb_p = drsdrb*dxdrs*decdx_p
+      decdrb_f = drsdrb*dxdrs*decdx_f
+      decdrb_a = drsdrb*dxdrs*decdx_a
+
+      decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+             + (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
+
+      EcLZ(2) = EcLZ(2) + weight(iG)*decdrb*r*r
+   
+    end if
+
+  end do
+
+end subroutine VWN5_lda_correlation_Levy_Zahariev_shift

src/xcDFT/VWN5_lda_correlation_energy.f90

@@ -0,0 +1,137 @@
+subroutine VWN5_lda_correlation_energy(nGrid,weight,rho,Ec)
+
+! Compute VWN5 LDA correlation energy
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,rs,x,z
+  double precision              :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
+  double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
+  double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
+
+  double precision              :: ec_z,ec_p,ec_f,ec_a
+  double precision              :: fz,d2fz
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Parameters of the functional
+
+  a_p  = +0.0621814D0/2D0
+  x0_p = -0.10498d0
+  b_p  = +3.72744d0
+  c_p  = +12.9352d0
+
+  a_f  = +0.0621814D0/4D0
+  x0_f = -0.325d0
+  b_f  = +7.06042d0
+  c_f  = +18.0578d0
+
+  a_a  = -1d0/(6d0*pi**2)
+  x0_a = -0.0047584D0
+  b_a  = 1.13107d0
+  c_a  = 13.0045d0
+
+! Initialization
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+
+!   alpha-alpha contribution
+
+    if(ra > threshold) then
+
+      rs = (4d0*pi*ra/3d0)**(-1d0/3d0)
+      x = sqrt(rs)
+
+      x_f   = x*x + b_f*x + c_f
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      q_f   = sqrt(4d0*c_f - b_f*b_f)
+
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+
+      Ec(1) = Ec(1) + weight(iG)*ec_f*ra
+
+    end if
+
+!   alpha-beta contribution
+
+    if(ra > threshold .and. rb > threshold) then
+
+      r = ra + rb
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      z = (ra - rb)/r
+      x = sqrt(rs)
+
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+
+      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+
+      x_p = x*x + b_p*x + c_p
+      x_f = x*x + b_f*x + c_f
+      x_a = x*x + b_a*x + c_a
+
+      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+
+      q_p = sqrt(4d0*c_p - b_p*b_p)
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+      q_a = sqrt(4d0*c_a - b_a*b_a)
+
+      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+
+      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+
+      Ec(2) = Ec(2) + weight(iG)*ec_z*r
+
+    end if
+
+!   beta-beta contribution
+
+    if(rb > threshold) then
+
+      rs = (4d0*pi*rb/3d0)**(-1d0/3d0)
+      x = sqrt(rs)
+
+      x_f   = x*x + b_f*x + c_f
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      q_f   = sqrt(4d0*c_f - b_f*b_f)
+
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) &
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+
+      Ec(3) = Ec(3) + weight(iG)*ec_f*rb
+
+    end if
+
+  end do
+
+  Ec(2) = Ec(2) - Ec(1) - Ec(3)
+
+end subroutine VWN5_lda_correlation_energy

src/xcDFT/VWN5_lda_correlation_individual_energy.f90

@@ -0,0 +1,204 @@
+subroutine VWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
+
+! Compute VWN5 LDA correlation potential
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rhow(nGrid,nspin)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,raI,rbI,rI,rs,x,z
+  double precision              :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
+  double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
+  double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
+  double precision              :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
+  double precision              :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
+  double precision              :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
+  double precision              :: ec_z,ec_p,ec_f,ec_a
+  double precision              :: fz,d2fz
+
+! Output variables
+
+  double precision              :: Ec(nspin)
+
+! Parameters of the functional
+
+  a_p  = +0.0621814D0/2D0
+  x0_p = -0.10498d0
+  b_p  = +3.72744d0
+  c_p  = +12.9352d0
+
+  a_f  = +0.0621814D0/4D0
+  x0_f = -0.325d0
+  b_f  = +7.06042d0
+  c_f  = +18.0578d0
+
+  a_a  = -1d0/(6d0*pi**2)
+  x0_a = -0.0047584D0
+  b_a  = +1.13107d0
+  c_a  = +13.0045d0
+
+! Initialization
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rhow(iG,1))
+    rb = max(0d0,rhow(iG,2))
+
+    raI = max(0d0,rho(iG,1))
+    rbI = max(0d0,rho(iG,2))
+   
+!   spin-up contribution
+   
+    if(ra > threshold .and. raI > threshold) then
+   
+      r = ra + rb
+      rI = raI + rbI
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      z = (ra - rb)/r
+      x = sqrt(rs)
+   
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+   
+      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+   
+      x_p = x*x + b_p*x + c_p
+      x_f = x*x + b_f*x + c_f
+      x_a = x*x + b_a*x + c_a
+   
+      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+   
+      q_p = sqrt(4d0*c_p - b_p*b_p)
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+      q_a = sqrt(4d0*c_a - b_a*b_a)
+   
+      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+   
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+   
+      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+   
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+
+      dzdra = (1d0 - z)/r
+      dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+      dfzdra = dzdra*dfzdz
+
+      drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_p = 2d0*x + b_p
+      dxdx_f = 2d0*x + b_f
+      dxdx_a = 2d0*x + b_a
+
+      decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                      - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                      - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+      decdra_p = drsdra*dxdrs*decdx_p
+      decdra_f = drsdra*dxdrs*decdx_f
+      decdra_a = drsdra*dxdrs*decdx_a
+
+      decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+             + (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
+
+      Ec(2) = Ec(2) + weight(iG)*(ec_z + decdra*r)*rI
+   
+    end if
+
+!   spin-down contribution
+   
+    if(rb > threshold .and. rbI > threshold) then
+   
+      r = ra + rb
+      rI = raI + rbI
+
+      rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+      z = (ra - rb)/r
+      x = sqrt(rs)
+   
+      fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+      fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+   
+      d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+   
+      x_p = x*x + b_p*x + c_p
+      x_f = x*x + b_f*x + c_f
+      x_a = x*x + b_a*x + c_a
+
+      xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+      xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+      xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+   
+      q_p = sqrt(4d0*c_p - b_p*b_p)
+      q_f = sqrt(4d0*c_f - b_f*b_f)
+      q_a = sqrt(4d0*c_a - b_a*b_a)
+   
+      ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                 - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+   
+      ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                 - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+   
+      ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                 - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+
+      ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+   
+      dzdrb = -(1d0 + z)/r
+      dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+      dfzdrb = dzdrb*dfzdz
+
+      drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+      dxdrs  = 0.5d0/sqrt(rs)
+
+      dxdx_p = 2d0*x + b_p
+      dxdx_f = 2d0*x + b_f
+      dxdx_a = 2d0*x + b_a
+
+      decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                      - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+      decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                      - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+      decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                      - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+      decdrb_p = drsdrb*dxdrs*decdx_p
+      decdrb_f = drsdrb*dxdrs*decdx_f
+      decdrb_a = drsdrb*dxdrs*decdx_a
+
+      decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+             + (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
+
+      Ec(2) = Ec(2) + weight(iG)*(ec_z + decdrb*r)*rI
+   
+    end if
+
+  end do
+
+end subroutine VWN5_lda_correlation_individual_energy

src/xcDFT/VWN5_lda_correlation_potential.f90

@@ -0,0 +1,202 @@
+subroutine VWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
+
+! Compute VWN5 LDA correlation 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,nspin)
+
+! Local variables
+
+  integer                       :: mu,nu,iG
+  double precision              :: ra,rb,r,rs,x,z
+  double precision              :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p
+  double precision              :: a_f,x0_f,xx0_f,b_f,c_f,x_f,q_f
+  double precision              :: a_a,x0_a,xx0_a,b_a,c_a,x_a,q_a
+  double precision              :: dfzdz,dxdrs,dxdx_p,dxdx_f,dxdx_a,decdx_p,decdx_f,decdx_a
+  double precision              :: dzdra,dfzdra,drsdra,decdra_p,decdra_f,decdra_a,decdra
+  double precision              :: dzdrb,dfzdrb,drsdrb,decdrb_p,decdrb_f,decdrb_a,decdrb
+
+  double precision              :: ec_z,ec_p,ec_f,ec_a
+  double precision              :: fz,d2fz
+
+! Output variables
+
+  double precision              :: Fc(nBas,nBas,nspin)
+
+! Parameters of the functional
+
+  a_p  = +0.0621814D0/2D0
+  x0_p = -0.10498d0
+  b_p  = +3.72744d0
+  c_p  = +12.9352d0
+
+  a_f  = +0.0621814D0/4D0
+  x0_f = -0.325d0
+  b_f  = +7.06042d0
+  c_f  = +18.0578d0
+
+  a_a  = -1d0/(6d0*pi**2)
+  x0_a = -0.0047584D0
+  b_a  = +1.13107d0
+  c_a  = +13.0045d0
+
+! Initialization
+
+  Fc(:,:,:) = 0d0
+
+  do mu=1,nBas
+    do nu=1,nBas
+      do iG=1,nGrid
+
+        ra = max(0d0,rho(iG,1))
+        rb = max(0d0,rho(iG,2))
+       
+!       spin-up contribution
+       
+        if(ra > threshold) then
+       
+          r = ra + rb
+          rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+          z = (ra - rb)/r
+          x = sqrt(rs)
+       
+          fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+          fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+       
+          d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+       
+          x_p = x*x + b_p*x + c_p
+          x_f = x*x + b_f*x + c_f
+          x_a = x*x + b_a*x + c_a
+       
+          xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+          xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+          xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+       
+          q_p = sqrt(4d0*c_p - b_p*b_p)
+          q_f = sqrt(4d0*c_f - b_f*b_f)
+          q_a = sqrt(4d0*c_a - b_a*b_a)
+       
+          ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                     - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+       
+          ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                     - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+       
+          ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                     - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+       
+          ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+
+          dzdra = (1d0 - z)/r
+          dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+          dfzdra = dzdra*dfzdz
+
+          drsdra = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+          dxdrs  = 0.5d0/sqrt(rs)
+
+          dxdx_p = 2d0*x + b_p
+          dxdx_f = 2d0*x + b_f
+          dxdx_a = 2d0*x + b_a
+
+          decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                          - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+          decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                          - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+          decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                          - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+          decdra_p = drsdra*dxdrs*decdx_p
+          decdra_f = drsdra*dxdrs*decdx_f
+          decdra_a = drsdra*dxdrs*decdx_a
+
+          decdra = decdra_p + decdra_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdra/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdra*z**3 &
+                 + (decdra_f - decdra_p)*fz*z**4 + (ec_f - ec_p)*dfzdra*z**4 + 4d0*(ec_f - ec_p)*fz*dzdra*z**3
+
+          Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdra*r)
+       
+        end if
+
+!       spin-down contribution
+       
+        if(rb > threshold) then
+       
+          r = ra + rb
+          rs = (4d0*pi*r/3d0)**(-1d0/3d0)
+          z = (ra - rb)/r
+          x = sqrt(rs)
+       
+          fz = (1d0 + z)**(4d0/3d0) + (1d0 - z)**(4d0/3d0) - 2d0
+          fz = fz/(2d0*(2d0**(1d0/3d0) - 1d0))
+       
+          d2fz = 4d0/(9d0*(2**(1d0/3d0) - 1d0))
+       
+          x_p = x*x + b_p*x + c_p
+          x_f = x*x + b_f*x + c_f
+          x_a = x*x + b_a*x + c_a
+
+          xx0_p = x0_p*x0_p + b_p*x0_p + c_p
+          xx0_f = x0_f*x0_f + b_f*x0_f + c_f
+          xx0_a = x0_a*x0_a + b_a*x0_a + c_a
+       
+          q_p = sqrt(4d0*c_p - b_p*b_p)
+          q_f = sqrt(4d0*c_f - b_f*b_f)
+          q_a = sqrt(4d0*c_a - b_a*b_a)
+       
+          ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & 
+                     - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) )
+       
+          ec_f = a_f*( log(x**2/x_f) + 2d0*b_f/q_f*atan(q_f/(2d0*x + b_f)) & 
+                     - b_f*x0_f/xx0_f*( log((x - x0_f)**2/x_f) + 2d0*(b_f + 2d0*x0_f)/q_f*atan(q_f/(2d0*x + b_f)) ) )
+       
+          ec_a = a_a*( log(x**2/x_a) + 2d0*b_a/q_a*atan(q_a/(2d0*x + b_a)) & 
+                     - b_a*x0_a/xx0_a*( log((x - x0_a)**2/x_a) + 2d0*(b_a + 2d0*x0_a)/q_a*atan(q_a/(2d0*x + b_a)) ) )
+
+          ec_z = ec_p + ec_a*fz/d2fz*(1d0-z**4) + (ec_f - ec_p)*fz*z**4
+       
+          dzdrb = - (1d0 + z)/r
+          dfzdz = (4d0/3d0)*((1d0 + z)**(1d0/3d0) - (1d0 - z)**(1d0/3d0))/(2d0*(2d0**(1d0/3d0) - 1d0))
+          dfzdrb = dzdrb*dfzdz
+
+          drsdrb = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0)
+          dxdrs  = 0.5d0/sqrt(rs)
+
+          dxdx_p = 2d0*x + b_p
+          dxdx_f = 2d0*x + b_f
+          dxdx_a = 2d0*x + b_a
+
+          decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p &
+                          - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) )
+
+          decdx_f = a_f*( 2d0/x - 4d0*b_f/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f &
+                          - b_f*x0_f/xx0_f*( 2/(x-x0_f) - 4d0*(b_f+2d0*x0_f)/( (b_f+2d0*x)**2 + q_f**2) - dxdx_f/x_f ) )
+
+          decdx_a = a_a*( 2d0/x - 4d0*b_a/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a &
+                          - b_a*x0_a/xx0_a*( 2/(x-x0_a) - 4d0*(b_a+2d0*x0_a)/( (b_a+2d0*x)**2 + q_a**2) - dxdx_a/x_a ) )
+
+          decdrb_p = drsdrb*dxdrs*decdx_p
+          decdrb_f = drsdrb*dxdrs*decdx_f
+          decdrb_a = drsdrb*dxdrs*decdx_a
+
+          decdrb = decdrb_p + decdrb_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 &
+                 + (decdrb_f - decdrb_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3
+          Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_z + decdrb*r)
+       
+        end if
+
+      end do
+    end do
+  end do
+
+end subroutine VWN5_lda_correlation_potential

src/xcDFT/W38_lda_correlation_Levy_Zahariev_shift.f90

@@ -0,0 +1,56 @@
+subroutine W38_lda_correlation_Levy_Zahariev_shift(nGrid,weight,rho,EcLZ)
+
+! Compute Wigner's LDA contribution to Levy-Zahariev shift
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r
+  double precision              :: a,d,epsc
+  double precision              :: dFcdra,dFcdrb
+
+! Output variables
+
+  double precision,intent(out)  :: EcLZ(nsp)
+
+! Coefficients for Wigner's LDA correlation
+
+  a = 0.04918d0
+  d = 0.349d0
+
+! Compute LDA correlation matrix in the AO basis
+
+  EcLZ(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+
+    r  = ra + rb
+
+    if(r > threshold) then
+
+      epsc  = ra*rb/(r + d*r**(2d0/3d0))
+      dFcdra = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 2d0/r + 1d0/ra)
+      dFcdrb = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 2d0/r + 1d0/rb)
+  
+      EcLZ(2) = EcLZ(2) - weight(iG)*r*r*(dFcdra + dFcdrb)
+
+    endif
+
+  enddo
+
+  EcLZ(:) = -4d0*a*EcLZ(:)
+
+end subroutine W38_lda_correlation_Levy_Zahariev_shift

src/xcDFT/W38_lda_correlation_energy.f90

@@ -0,0 +1,52 @@
+subroutine W38_lda_correlation_energy(nGrid,weight,rho,Ec)
+
+! Compute Wigner's LDA correlation energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r
+  double precision              :: a,d,epsc
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Coefficients for Wigner's LDA correlation
+
+  a = 0.04918d0
+  d = 0.349d0
+
+! Compute LDA correlation energy
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+
+    if(ra > threshold .and. rb > threshold) then
+
+      r = ra + rb
+
+      epsc = ra*rb/(r + d*r**(2d0/3d0))
+
+      Ec(2) = Ec(2) + weight(iG)*epsc
+
+    endif
+
+  enddo
+
+  Ec(2) = -4d0*a*Ec(2)
+
+end subroutine W38_lda_correlation_energy

src/xcDFT/W38_lda_correlation_individual_energy.f90

@@ -0,0 +1,62 @@
+subroutine W38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec)
+
+! Compute Wigner's LDA individual energy
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rhow(nGrid,nspin)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r
+  double precision              :: raI,rbI,rI
+  double precision              :: a,d,epsc
+  double precision              :: dFcdra,dFcdrb
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp)
+
+! Coefficients for Wigner's LDA correlation
+
+  a = 0.04918d0
+  d = 0.349d0
+
+! Compute LDA correlation individual energy
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rhow(iG,1))
+    rb = max(0d0,rhow(iG,2))
+
+    raI = max(0d0,rho(iG,1))
+    rbI = max(0d0,rho(iG,2))
+    
+    r  = ra + rb
+    rI = raI + rbI
+
+    if(r > threshold .and. rI > threshold) then
+
+      epsc  = ra*rb/(r + d*r**(2d0/3d0))
+      dFcdra = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/ra)
+      dFcdrb = epsc*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/rb)
+  
+      Ec(2) = Ec(2) + weight(iG)*rI*0.5d0*(dFcdra + dFcdrb)
+
+    endif
+
+  enddo
+
+  Ec(2) = -4d0*a*Ec(2)
+
+end subroutine W38_lda_correlation_individual_energy

src/xcDFT/W38_lda_correlation_potential.f90

@@ -0,0 +1,76 @@
+subroutine W38_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
+
+! Compute Wigner's LDA correlation 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,nspin)
+
+! Local variables
+
+  integer                       :: mu,nu,iG
+  double precision              :: ra,rb,r
+  double precision              :: a,d,ec
+  double precision              :: dFcdr
+
+! Output variables
+
+  double precision,intent(out)  :: Fc(nBas,nBas,nspin)
+
+! Coefficients for Wigner's LDA correlation
+
+  a = 0.04918d0
+  d = 0.349d0
+
+! Compute LDA correlation matrix in the AO basis
+
+  Fc(:,:,:) = 0d0
+
+  do mu=1,nBas
+    do nu=1,nBas
+      do iG=1,nGrid
+
+        ra = max(0d0,rho(iG,1))
+        rb = max(0d0,rho(iG,2))
+      
+!       Spin-up part contribution
+
+        if(ra > threshold) then
+
+          r = ra + rb
+
+          ec    = ra*rb/(r + d*r**(2d0/3d0))
+          dFcdr = ec*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/ra)
+  
+          Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdr
+
+        endif
+
+!       Spin-down part contribution
+
+        if(rb > threshold) then
+
+          r = ra + rb
+
+          ec    = ra*rb/(r + d*r**(2d0/3d0))
+          dFcdr = ec*(d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) - 1d0/r + 1d0/rb)
+  
+          Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdr
+
+        endif
+
+      enddo
+
+    enddo
+  enddo
+
+  Fc(:,:,:) = -4d0*a*Fc(:,:,:)
+
+end subroutine W38_lda_correlation_potential

src/xcDFT/correlation_Levy_Zahariev_shift.f90

@@ -0,0 +1,69 @@
+subroutine correlation_Levy_Zahariev_shift(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,EcLZ)
+
+! Compute the correlation part of the Levy-Zahariev shift
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  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,nspin)
+  double precision,intent(in)   :: drho(ncart,nGrid,nspin)
+
+! Local variables
+
+  double precision              :: EcLZLDA(nsp),EcLZGGA(nsp)
+  double precision              :: aC
+
+! Output variables
+
+  double precision,intent(out)  :: EcLZ(nsp)
+
+  select case (rung)
+
+!   Hartree calculation
+
+    case(0) 
+
+      EcLZ(:) = 0d0
+
+!   LDA functionals
+
+    case(1) 
+
+      call lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),EcLZ(:))
+
+!   GGA functionals
+
+    case(2) 
+
+      call print_warning('!!! Individual energies NYI for GGAs !!!')
+      stop
+
+!   Hybrid functionals
+
+    case(4) 
+
+      call print_warning('!!! Individual energies NYI for hybrids !!!')
+      stop
+
+      aC = 0.81d0
+
+      EcLZ(:) = EcLZLDA(:) + aC*(EcLZGGA(:) - EcLZLDA(:)) 
+
+!   Hartree-Fock calculation
+
+    case(666) 
+
+      EcLZ(:) = 0d0
+
+  end select
+ 
+end subroutine correlation_Levy_Zahariev_shift

src/xcDFT/correlation_derivative_discontinuity.f90

@@ -0,0 +1,65 @@
+subroutine correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
+
+! Compute the correlation part of the derivative discontinuity
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  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,nspin)
+  double precision,intent(in)   :: drhow(ncart,nGrid,nspin)
+
+! Local variables
+
+  double precision              :: aC
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp,nEns)
+
+  select case (rung)
+
+!   Hartree calculation
+
+    case(0) 
+
+      Ec(:,:) = 0d0
+
+!   LDA functionals
+
+    case(1) 
+
+      call lda_correlation_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),Ec(:,:))
+
+!   GGA functionals
+
+    case(2) 
+
+      call print_warning('!!! derivative discontinuity NYI for GGAs !!!')
+      stop
+
+!   Hybrid functionals
+
+    case(4) 
+
+      call print_warning('!!! derivative discontinuity NYI for hybrids !!!')
+      stop
+
+      aC = 0.81d0
+
+!   Hartree-Fock calculation
+
+    case(666) 
+
+      Ec(:,:) = 0d0
+
+  end select
+ 
+end subroutine correlation_derivative_discontinuity

src/xcDFT/correlation_energy.f90

@@ -0,0 +1,68 @@
+subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
+
+! Compute the correlation energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  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,nspin)
+  double precision,intent(in)   :: drho(3,nGrid,nspin)
+
+! Local variables
+
+  double precision              :: EcLDA(nsp)
+  double precision              :: EcGGA(nsp)
+  double precision              :: aC
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp)
+
+  select case (rung)
+
+!   Hartree calculation
+
+    case(0) 
+
+      Ec(:) = 0d0
+
+!   LDA functionals
+
+    case(1) 
+
+      call lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),Ec(:))
+
+!   GGA functionals
+
+    case(2) 
+
+      call gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),Ec(:))
+
+!   Hybrid functionals
+
+    case(4) 
+
+      aC = 0.81d0
+
+      call lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),EcLDA(:))
+      call gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),EcGGA(:))
+
+      Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) 
+
+!   Hartree-Fock calculation
+
+    case(666) 
+
+      Ec(:) = 0d0
+
+  end select
+ 
+end subroutine correlation_energy

src/xcDFT/correlation_individual_energy.f90

@@ -0,0 +1,69 @@
+subroutine correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec)
+
+! Compute the correlation energy of individual states
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  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,nspin)
+  double precision,intent(in)   :: drhow(ncart,nGrid,nspin)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+  double precision,intent(in)   :: drho(ncart,nGrid,nspin)
+
+! Local variables
+
+  double precision              :: EcLDA(nsp)
+  double precision              :: EcGGA(nsp)
+  double precision              :: aC
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp)
+
+  select case (rung)
+
+!   Hartree calculation
+
+    case(0) 
+
+      Ec(:) = 0d0
+
+!   LDA functionals
+
+    case(1) 
+
+      call lda_correlation_individual_energy(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:))
+
+!   GGA functionals
+
+    case(2) 
+
+      call print_warning('!!! Individual energies NYI for GGAs !!!')
+      stop
+
+!   Hybrid functionals
+
+    case(4) 
+
+      call print_warning('!!! Individual energies NYI for hybrids !!!')
+      stop
+
+      aC = 0.81d0
+
+!   Hartree-Fock calculation
+
+    case(666) 
+
+      Ec(:) = 0d0
+
+  end select
+ 
+end subroutine correlation_individual_energy

src/xcDFT/correlation_potential.f90

@@ -0,0 +1,75 @@
+subroutine correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
+
+! Compute the correlation potential 
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  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(ncart,nBas,nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+  double precision,intent(in)   :: drho(ncart,nGrid,nspin)
+
+! Local variables
+
+  double precision,allocatable  :: FcLDA(:,:,:)
+  double precision,allocatable  :: FcGGA(:,:,:)
+  double precision              :: aC
+
+! Output variables
+
+  double precision,intent(out)  :: Fc(nBas,nBas,nspin)
+
+! Memory allocation
+
+  select case (rung)
+
+!   Hartree calculation
+
+    case(0) 
+
+      Fc(:,:,:) = 0d0
+
+!   LDA functionals
+
+    case(1) 
+
+      call lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
+
+!   GGA functionals
+
+    case(2) 
+
+      call gga_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rho(:,:),drho(:,:,:),Fc(:,:,:))
+
+!   Hybrid functionals
+
+    case(4) 
+
+      allocate(FcLDA(nBas,nBas,nspin),FcGGA(nBas,nBas,nspin))
+
+      aC = 0.81d0
+
+      call lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:,:),FcLDA(:,:,:))
+      call gga_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rho(:,:),drho(:,:,:),FcGGA(:,:,:))
+
+      Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:))  
+
+!   Hartree-Fock calculation
+
+    case(666) 
+
+      Fc(:,:,:) = 0d0
+
+  end select
+
+end subroutine correlation_potential

src/xcDFT/density.f90

@@ -0,0 +1,38 @@
+subroutine density(nGrid,nBas,P,AO,rho)
+
+! Calculate one-electron density
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,parameter    :: thresh = 1d-15
+
+  integer,intent(in)            :: nGrid
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: AO(nBas,nGrid)
+
+! Local variables
+
+  integer                       :: iG,mu,nu
+
+! Output variables
+
+  double precision,intent(out)  :: rho(nGrid) 
+
+  rho(:) = 0d0
+  do iG=1,nGrid
+    do mu=1,nBas
+      do nu=1,nBas
+        rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
+      enddo
+    enddo
+  enddo
+
+! do iG=1,nGrid
+!   rho(iG) = max(rho(iG),thresh)
+! enddo
+
+end subroutine density

src/xcDFT/density_matrix.f90

@@ -0,0 +1,43 @@
+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/xcDFT/eDFT.f90

@@ -0,0 +1,133 @@
+program eDFT
+
+! exchange-correlation density-functional theory calculations
+
+  include 'parameters.h'
+
+  integer                       :: nAt,nBas,nEl(nspin),nO(nspin),nV(nspin)
+  double precision              :: ENuc,EKS
+
+  double precision,allocatable  :: ZNuc(:),rAt(:,:)
+
+  integer                       :: nShell
+  integer,allocatable           :: TotAngMomShell(:)
+  integer,allocatable           :: KShell(:)
+  double precision,allocatable  :: CenterShell(:,:)
+  double precision,allocatable  :: DShell(:,:)
+  double precision,allocatable  :: ExpShell(:,:)
+
+  double precision,allocatable  :: S(:,:),T(:,:),V(:,:),Hc(:,:),X(:,:)
+  double precision,allocatable  :: ERI(:,:,:,:)
+
+  integer                       :: x_rung,c_rung
+  character(len=12)             :: x_DFA ,c_DFA
+  integer                       :: SGn
+  integer                       :: nRad,nAng,nGrid
+  double precision,allocatable  :: root(:,:)
+  double precision,allocatable  :: weight(:)
+  double precision,allocatable  :: AO(:,:)
+  double precision,allocatable  :: dAO(:,:,:)
+
+  double precision              :: start_KS,end_KS,t_KS
+
+  integer                       :: nEns
+  double precision,allocatable  :: wEns(:)
+
+  integer                       :: maxSCF,max_diis
+  double precision              :: thresh
+  logical                       :: DIIS,guess_type,ortho_type
+
+! Hello World
+
+  write(*,*)
+  write(*,*) '******************************************'
+  write(*,*) '* eDFT: density-functional for ensembles *'
+  write(*,*) '******************************************'
+  write(*,*)
+
+!------------------------------------------------------------------------
+! Read input information
+!------------------------------------------------------------------------
+
+! Read number of atoms, number of electrons of the system
+! nO   = number of occupied orbitals
+! nV   = number of virtual orbitals (see below)
+! nBas = number of basis functions (see below)
+!      = nO + nV
+
+  call read_molecule(nAt,nEl,nO)
+  allocate(ZNuc(nAt),rAt(nAt,ncart))
+
+! Read geometry
+
+  call read_geometry(nAt,ZNuc,rAt,ENuc)
+
+  allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), &
+           DShell(maxShell,maxK),ExpShell(maxShell,maxK))
+
+!------------------------------------------------------------------------
+! Read basis set information
+!------------------------------------------------------------------------
+
+  call read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
+
+!------------------------------------------------------------------------
+! Read one- and two-electron integrals
+!------------------------------------------------------------------------
+
+! Memory allocation for one- and two-electron integrals
+
+  allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), &
+           ERI(nBas,nBas,nBas,nBas))
+
+!   Read integrals
+
+  call read_integrals(nBas,S,T,V,Hc,ERI)  
+
+! Orthogonalization X = S^(-1/2)
+
+  call orthogonalization_matrix(nBas,S,X)
+
+!------------------------------------------------------------------------
+! DFT options
+!------------------------------------------------------------------------
+
+! Allocate ensemble weights
+
+  allocate(wEns(maxEns))
+
+  call read_options(x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
+
+!------------------------------------------------------------------------
+! Construct quadrature grid
+!------------------------------------------------------------------------
+  call read_grid(SGn,nRad,nAng,nGrid)
+
+  allocate(root(ncart,nGrid),weight(nGrid))
+  call quadrature_grid(nRad,nAng,nGrid,root,weight)
+
+!------------------------------------------------------------------------
+! Calculate AO values at grid points
+!------------------------------------------------------------------------
+
+  allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid))
+  call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
+                      nGrid,root,AO,dAO)
+
+!------------------------------------------------------------------------
+! Compute KS energy
+!------------------------------------------------------------------------
+
+    call cpu_time(start_KS)
+    call Kohn_Sham(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(1:nEns),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & 
+                   nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS)
+    call cpu_time(end_KS)
+
+    t_KS = end_KS - start_KS
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for KS = ',t_KS,' seconds'
+    write(*,*)
+
+!------------------------------------------------------------------------
+! End of eDFT
+!------------------------------------------------------------------------
+end program eDFT

src/xcDFT/elda_correlation_Levy_Zahariev_shift.f90

@@ -0,0 +1,62 @@
+subroutine elda_correlation_Levy_Zahariev_shift(nEns,aLF,nGrid,weight,rho,EcLZ)
+
+! Compute Levy-Zahariev LDA correlation shift of 2-glomium for various states
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nEns
+  double precision,intent(in)   :: aLF(nEns)
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,ec_p
+  double precision              :: dFcdra,dFcdrb
+
+! Output variables
+
+  double precision,intent(out)  :: EcLZ
+
+! Compute Levy-Zahariev eLDA correlation shift
+
+  EcLZ = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+    r = ra + rb
+
+    if(ra > threshold) then
+
+      ec_p   = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+      dFcdra = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
+      dFcdra = dFcdra/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+      dFcdra = ec_p*dFcdra/6d0
+      dFcdra = ec_p + dFcdra
+
+      EcLZ = EcLZ - weight(iG)*r*r*dFcdra
+
+     end if
+
+    if(rb > threshold) then
+
+      ec_p   = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+      dFcdrb = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
+      dFcdrb = dFcdrb/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+      dFcdrb = ec_p*dFcdrb/6d0
+      dFcdrb = ec_p + dFcdrb
+
+      EcLZ = EcLZ - weight(iG)*r*r*dFcdrb
+
+     end if
+
+  end do
+
+end subroutine elda_correlation_Levy_Zahariev_shift

src/xcDFT/elda_correlation_energy.f90

@@ -0,0 +1,70 @@
+subroutine elda_correlation_energy(nEns,aLF,nGrid,weight,rho,Ec)
+
+! Compute LDA correlation energy of 2-glomium for various states
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nEns
+  double precision,intent(in)   :: aLF(nEns)
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,r,ec_p
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp)
+
+
+! Compute eLDA correlation energy
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+    r  = ra + rb
+
+!   Spin-up contribution
+
+    if(ra > threshold) then
+
+      ec_p = aLF(1)/(1d0 + aLF(2)*ra**(-1d0/6d0) + aLF(3)*ra**(-1d0/3d0))
+
+      Ec(1) = Ec(1) + weight(iG)*ec_p*ra
+
+    end if
+
+!   Opposite-spin contribution
+
+    if(r > threshold) then
+
+      ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+
+      Ec(2) = Ec(2) + weight(iG)*ec_p*r
+
+    end if
+
+!   Spin-down contribution
+
+    if(rb > threshold) then
+
+      ec_p = aLF(1)/(1d0 + aLF(2)*rb**(-1d0/6d0) + aLF(3)*rb**(-1d0/3d0))
+
+      Ec(3) = Ec(3) + weight(iG)*ec_p*rb
+
+    end if
+
+  end do
+
+  Ec(2) = Ec(2) - Ec(1) - Ec(3)
+
+end subroutine elda_correlation_energy

src/xcDFT/elda_correlation_individual_energy.f90

@@ -0,0 +1,58 @@
+subroutine elda_correlation_individual_energy(nEns,aLF,nGrid,weight,rhow,rho,Ec)
+
+! Compute LDA correlation individual energy of 2-glomium for various states
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nEns
+  double precision,intent(in)   :: aLF(nEns)
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: weight(nGrid)
+  double precision,intent(in)   :: rhow(nGrid,nspin)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra, rb, r
+  double precision              :: raI,rbI,rI
+  double precision              :: ec_p,dFcdr
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp)
+
+! Compute eLDA correlation potential
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = max(0d0,rho(iG,1))
+    rb = max(0d0,rho(iG,2))
+
+    raI = max(0d0,rho(iG,1))
+    rbI = max(0d0,rho(iG,2))
+
+    r  = ra  + rb
+    rI = raI + rbI
+
+    if(r > threshold .and. rI > threshold) then
+
+      ec_p  = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+
+      dFcdr = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
+      dFcdr = dFcdr/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+      dFcdr = ec_p*dFcdr/(6d0*r)
+      dFcdr = ec_p + dFcdr*r
+
+      Ec(2) = Ec(2) + weight(iG)*rI*dFcdr
+
+     end if
+
+  end do
+
+end subroutine elda_correlation_individual_energy

src/xcDFT/elda_correlation_potential.f90

@@ -0,0 +1,71 @@
+subroutine elda_correlation_potential(nEns,aLF,nGrid,weight,nBas,AO,rho,Fc)
+
+! Compute LDA correlation energy of 2-glomium for various states
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nEns
+  double precision,intent(in)   :: aLF(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,nspin)
+
+! Local variables
+
+  integer                       :: mu,nu,iG
+  double precision              :: ra,rb,r,ec_p
+  double precision              :: dFcdra,dFcdrb
+
+! Output variables
+
+  double precision,intent(out)  :: Fc(nBas,nBas,nspin)
+
+! Compute eLDA correlation potential
+
+  Fc(:,:,:) = 0d0
+
+  do mu=1,nBas
+    do nu=1,nBas
+      do iG=1,nGrid
+
+        ra = max(0d0,rho(iG,1))
+        rb = max(0d0,rho(iG,2))
+
+        if(ra > threshold) then
+
+          r = ra + rb
+
+          ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+          dFcdra = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
+          dFcdra = dFcdra/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+          dFcdra = ec_p*dFcdra/(6d0*r)
+          dFcdra = ec_p + dFcdra*r
+
+          Fc(mu,nu,1) = Fc(mu,nu,1) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdra
+
+         endif
+
+        if(rb > threshold) then
+
+          r = ra + rb
+
+          ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+          dFcdrb = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
+          dFcdrb = dFcdrb/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
+          dFcdrb = ec_p*dFcdrb/(6d0*r)
+          dFcdrb = ec_p + dFcdrb*r
+
+          Fc(mu,nu,2) = Fc(mu,nu,2) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdrb
+
+         endif
+
+      end do
+    end do
+  end do
+
+end subroutine elda_correlation_potential

src/xcDFT/electron_number.f90

@@ -0,0 +1,20 @@
+function electron_number(nGrid,w,rho) result(nEl)
+
+! Compute the number of electrons via integration of the one-electron density
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  double precision,intent(in)   :: w(nGrid)
+  double precision,intent(in)   :: rho(nGrid)
+
+! Output variables
+
+  double precision              :: nEl
+
+  nEl = dot_product(w,rho)
+
+end function electron_number

src/xcDFT/elements.f90

@@ -0,0 +1,170 @@
+function element_number(element_name)
+
+ implicit none
+
+ integer,parameter                  :: nelement_max = 103
+ character(len=2),intent(in)        :: element_name
+ integer                            :: element_number
+ character(len=2),parameter         :: element_list(nelement_max) =                            &
+  (/' H',                                                                                'He', &  !   2
+    'Li','Be',                                                  ' B',' C',' N',' O',' F','Ne', &  !  10
+    'Na','Mg',                                                  'Al','Si',' P',' S','Cl','Ar', &  !  18
+    ' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', &  !  36
+    'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe', &  !  54
+    'Cs','Ba',                                                                                 &  !  56
+    'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb',                     &  !  70
+              'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn', &  !  86
+    'Fr','Ra',                                                                                 &  !  88
+              'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No',           &  ! 102
+    'Lr'                                                                                       &  ! 103
+  /)
+
+!=====
+ integer :: ielement
+!=====
+
+ ielement=1
+ do while( ADJUSTL(element_name) /= ADJUSTL(element_list(ielement)) )
+   if( ielement == nelement_max ) then
+     write(*,'(a,a)')    ' Input symbol ',element_name
+     write(*,'(a,i3,a)') ' Element symbol is not one of first ',nelement_max,' elements'
+     write(*,*) '!!! element symbol not understood !!!'
+     stop
+   endif
+   ielement = ielement + 1
+ enddo
+
+ element_number = ielement
+
+end function element_number
+
+function element_core(zval,zatom)
+ implicit none
+ double precision,intent(in) :: zval
+ double precision,intent(in) :: zatom
+ integer                     :: element_core
+!=====
+
+ !
+ ! If zval /= zatom, this is certainly an effective core potential
+ ! and no core states should be frozen.
+ if( ABS(zval - zatom) > 1d0-3 ) then
+   element_core = 0
+ else
+
+   if( zval <= 4.00001d0 ) then  ! up to Be
+     element_core = 0
+   else if( zval <= 12.00001d0 ) then  ! up to Mg
+     element_core = 1
+   else if( zval <= 30.00001d0 ) then  ! up to Ca
+     element_core = 5
+   else if( zval <= 48.00001d0 ) then  ! up to Sr
+     element_core = 9
+   else
+     write(*,*) '!!! not imlemented in element_core !!!'
+     stop
+   endif
+
+ endif
+
+
+end function element_core
+
+function element_covalent_radius(zatom)
+
+! Return covalent radius of an atom
+
+ implicit none
+ include 'parameters.h'
+
+ integer,intent(in) :: zatom
+ double precision   :: element_covalent_radius
+
+ !
+ ! Data from Cambridge Structural Database
+ ! http://en.wikipedia.org/wiki/Covalent_radius
+ !
+ ! Values are first given in picometer
+ ! They will be converted in bohr later on
+ select case(zatom)
+ case( 1)
+   element_covalent_radius =  31.
+ case( 2)
+   element_covalent_radius =  28.
+ case( 3)
+   element_covalent_radius = 128.
+ case( 4)
+   element_covalent_radius =  96.
+ case( 5)
+   element_covalent_radius =  84.
+ case( 6)
+   element_covalent_radius =  73.
+ case( 7)
+   element_covalent_radius =  71.
+ case( 8)
+   element_covalent_radius =  66.
+ case( 9)
+   element_covalent_radius =  57.
+ case(10) ! Ne.
+   element_covalent_radius =  58.
+ case(11)
+   element_covalent_radius = 166.
+ case(12)
+   element_covalent_radius = 141.
+ case(13)
+   element_covalent_radius = 121.
+ case(14)
+   element_covalent_radius = 111.
+ case(15)
+   element_covalent_radius = 107.
+ case(16)
+   element_covalent_radius = 105.
+ case(17)
+   element_covalent_radius = 102.
+ case(18) ! Ar.
+   element_covalent_radius = 106.
+ case(19)
+   element_covalent_radius = 203.
+ case(20)
+   element_covalent_radius = 176.
+ case(21)
+   element_covalent_radius = 170.
+ case(22)
+   element_covalent_radius = 160.
+ case(23)
+   element_covalent_radius = 153.
+ case(24)
+   element_covalent_radius = 139.
+ case(25)
+   element_covalent_radius = 145.
+ case(26)
+   element_covalent_radius = 145.
+ case(27)
+   element_covalent_radius = 140.
+ case(28)
+   element_covalent_radius = 124.
+ case(29)
+   element_covalent_radius = 132.
+ case(30)
+   element_covalent_radius = 122.
+ case(31)
+   element_covalent_radius = 120.
+ case(32)
+   element_covalent_radius = 119.
+ case(34)
+   element_covalent_radius = 120.
+ case(35)
+   element_covalent_radius = 120.
+ case(36) ! Kr.
+   element_covalent_radius = 116.
+ case default
+   write(*,*) '!!! covalent radius not available !!!'
+   stop
+ end select
+
+ ! pm to bohr conversion
+ element_covalent_radius = element_covalent_radius*pmtoau
+
+
+end function element_covalent_radius
+

src/xcDFT/exchange_energy.f90

@@ -0,0 +1,80 @@
+subroutine exchange_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,FxHF,rho,drho,Ex)
+
+! Compute the exchange energy
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  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)   :: P(nBas,nBas)
+  double precision,intent(in)   :: FxHF(nBas,nBas)
+  double precision,intent(in)   :: rho(nGrid)
+  double precision,intent(in)   :: drho(3,nGrid)
+
+! Local variables
+
+  double precision              :: ExLDA,ExGGA,ExHF
+  double precision              :: cX,aX,aC
+
+! Output variables
+
+  double precision,intent(out)  :: Ex
+
+  select case (rung)
+
+!   Hartree calculation
+
+    case(0) 
+
+      Ex = 0d0
+
+!   LDA functionals
+
+    case(1) 
+
+      call lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA)
+
+      Ex = ExLDA
+
+!   GGA functionals
+
+    case(2) 
+
+      call gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
+
+      Ex = ExGGA
+
+!   Hybrid functionals
+
+    case(4) 
+
+      cX = 0.20d0
+      aX = 0.72d0
+      aC = 0.81d0
+
+      call lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA)
+      call gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
+      call fock_exchange_energy(nBas,P,FxHF,ExHF)
+
+      Ex = ExLDA              & 
+         + cX*(ExHF  - ExLDA) &
+         + aX*(ExGGA - ExLDA) 
+
+!   Hartree-Fock calculation
+    case(666) 
+
+      call fock_exchange_energy(nBas,P,FxHF,ExHF)
+
+      Ex = ExHF
+
+  end select
+ 
+end subroutine exchange_energy

src/xcDFT/exchange_potential.f90

@@ -0,0 +1,81 @@
+subroutine exchange_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,P,ERI,AO,dAO,rho,drho,Fx,FxHF)
+
+! Compute the exchange potential 
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  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)   :: P(nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,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
+
+  double precision,allocatable  :: FxLDA(:,:),FxGGA(:,:)
+  double precision              :: cX,aX
+
+! Output variables
+
+  double precision,intent(out)  :: Fx(nBas,nBas),FxHF(nBas,nBas)
+
+! Memory allocation
+
+  select case (rung)
+
+!   Hartree calculation
+    case(0) 
+
+      Fx(:,:) = 0d0
+
+!   LDA functionals
+
+    case(1) 
+
+      call lda_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
+
+!   GGA functionals
+
+    case(2) 
+
+      call gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+
+!   Hybrid functionals
+
+    case(4) 
+
+      allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas))
+
+      cX = 0.20d0
+      aX = 0.72d0
+
+      call lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA)
+      call gga_exchange_potential(DFA,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
+      call fock_exchange_potential(nBas,P,ERI,FxHF)
+
+      Fx(:,:) = FxLDA(:,:)                   &
+              + cX*(FxHF(:,:)  - FxLDA(:,:)) &
+              + aX*(FxGGA(:,:) - FxLDA(:,:))  
+
+!   Hartree-Fock calculation
+
+    case(666) 
+
+      call fock_exchange_potential(nBas,P,ERI,FxHF)
+
+      Fx(:,:) = FxHF(:,:)
+
+  end select
+
+end subroutine exchange_potential

src/xcDFT/fock_exchange_energy.f90

@@ -0,0 +1,25 @@
+subroutine fock_exchange_energy(nBas,P,Fx,Ex)
+
+! Compute the (exact) Fock exchange energy
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: Fx(nBas,nBas)
+
+! Local variables
+
+  double precision,external     :: trace_matrix
+
+! Output variables
+
+  double precision,intent(out)  :: Ex
+
+! Compute HF exchange energy
+
+  Ex = trace_matrix(nBas,matmul(P,Fx))
+
+end subroutine fock_exchange_energy

src/xcDFT/fock_exchange_potential.f90

@@ -0,0 +1,34 @@
+subroutine fock_exchange_potential(nBas,P,ERI,Fx)
+
+! Compute the Fock exchange potential
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: mu,nu,la,si
+
+! Output variables
+
+  double precision,intent(out)  :: Fx(nBas,nBas)
+
+! Compute HF exchange matrix
+
+  Fx(:,:) = 0d0
+  do nu=1,nBas
+    do si=1,nBas
+      do la=1,nBas
+        do mu=1,nBas
+          Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,si,la,nu)
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine fock_exchange_potential

src/xcDFT/generate_shell.f90

@@ -0,0 +1,32 @@
+subroutine generate_shell(atot,nShellFunction,ShellFunction)
+
+! Generate shells for a given total angular momemtum
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)   :: atot,nShellFunction
+
+! Local variables
+
+  integer              :: ax,ay,az,ia
+
+! Output variables
+
+  integer,intent(out)  :: ShellFunction(nShellFunction,3)
+
+  ia = 0
+  do ax=atot,0,-1
+    do az=0,atot
+      ay = atot - ax - az
+      if(ay >= 0) then
+        ia = ia + 1
+        ShellFunction(ia,1) = ax
+        ShellFunction(ia,2) = ay
+        ShellFunction(ia,3) = az
+      endif
+    enddo
+  enddo
+
+end subroutine generate_shell

src/xcDFT/gga_correlation_energy.f90

@@ -0,0 +1,42 @@
+subroutine gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec)
+
+! Compute GGA correlation energy
+
+  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)   :: rho(nGrid,nspin)
+  double precision,intent(in)   :: drho(ncart,nGrid,nspin)
+
+! Local variables
+
+  integer                       :: iG
+  double precision              :: ra,rb,ga,gb
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Coefficients for ??? GGA exchange functional
+
+! Compute GGA exchange energy
+
+  Ec(:) = 0d0
+
+  do iG=1,nGrid
+
+    ra = rho(iG,1)
+    rb = rho(iG,2)
+    ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**2
+    gb = drho(1,iG,2)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2
+
+  enddo
+
+end subroutine gga_correlation_energy

src/xcDFT/gga_correlation_potential.f90

@@ -0,0 +1,31 @@
+subroutine gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
+
+! Compute GGA 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)
+
+! Coefficients for GGA correlation functional
+
+! Compute GGA correlation matrix in the AO basis
+
+end subroutine gga_correlation_potential

src/xcDFT/gga_exchange_energy.f90

@@ -0,0 +1,44 @@
+subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
+
+! Select GGA exchange functional for energy 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)
+  double precision,intent(in)   :: rho(nGrid)
+  double precision,intent(in)   :: drho(3,nGrid)
+
+! Output variables
+
+  double precision              :: Ex
+
+  select case (DFA)
+
+!   Gill's 96 exchange functional 
+
+    case ('G96')
+
+      call G96_gga_exchange_energy(DFA,nEns,wEns,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)
+
+    case default
+
+      call print_warning('!!! GGA exchange functional not available !!!')
+      stop
+
+  end select
+
+end subroutine gga_exchange_energy

src/xcDFT/gga_exchange_potential.f90

@@ -0,0 +1,48 @@
+subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+
+! Select GGA 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 GGA exchange functional
+
+  select case (DFA)
+
+!   Gill's 96 exchange functional 
+
+    case ('G96')
+
+      call G96_gga_exchange_potential(DFA,nEns,wEns,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)
+
+    case default
+
+      call print_warning('!!! GGA exchange functional not available !!!')
+      stop
+
+  end select
+
+end subroutine gga_exchange_potential

src/xcDFT/gradient_density.f90

@@ -0,0 +1,45 @@
+subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho)
+
+! Calculate gradient of the one-electron density
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,parameter    :: thresh = 1d-15
+
+  integer,intent(in)            :: nGrid
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: AO(nBas,nGrid)
+  double precision,intent(in)   :: dAO(3,nBas,nGrid)
+
+! Local variables
+
+  integer                       :: ixyz,iG,mu,nu
+  double precision,external     :: trace_matrix
+
+! Output variables
+
+  double precision,intent(out)  :: drho(3,nGrid)
+
+  drho(:,:) = 0d0
+  do iG=1,nGrid
+    do mu=1,nBas
+      do nu=1,nBas
+        do ixyz=1,3
+          drho(ixyz,iG) = drho(ixyz,iG) & 
+                        + P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG))
+        enddo
+      enddo
+    enddo
+  enddo
+
+  do iG=1,nGrid
+    do ixyz=1,3
+      if(abs(drho(ixyz,iG)) < thresh) drho(ixyz,iG) = thresh
+    enddo
+  enddo
+
+end subroutine gradient_density

src/xcDFT/hartree_coulomb.f90

@@ -0,0 +1,33 @@
+subroutine hartree_coulomb(nBas,P,ERI,J)
+
+! Compute Coulomb matrix
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: mu,nu,la,si
+
+! Output variables
+
+  double precision,intent(out)  :: J(nBas,nBas)
+
+  J = 0d0
+  do mu=1,nBas
+    do nu=1,nBas
+      do la=1,nBas
+        do si=1,nBas
+          J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,nu,la,si)
+        enddo
+      enddo
+    enddo
+  enddo
+
+
+end subroutine hartree_coulomb

src/xcDFT/individual_energy.f90

@@ -0,0 +1,156 @@
+subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO, & 
+                             nO,nV,T,V,ERI,ENuc,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,E,Om)
+
+! Compute individual energies as well as excitation energies
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: x_rung,c_rung
+  character(len=12),intent(in)  :: x_DFA,c_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(ncart,nBas,nGrid)
+
+  integer,intent(in)            :: nO(nspin),nV(nspin)
+  double precision,intent(in)   :: T(nBas,nBas)
+  double precision,intent(in)   :: V(nBas,nBas)
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ENuc
+
+  double precision,intent(in)   :: Pw(nBas,nBas)
+  double precision,intent(in)   :: rhow(nGrid,nspin)
+  double precision,intent(in)   :: drhow(ncart,nGrid,nspin)
+
+  double precision,intent(in)   :: P(nBas,nBas,nspin,nEns)
+  double precision,intent(in)   :: rho(nGrid,nspin,nEns)
+  double precision,intent(in)   :: drho(ncart,nGrid,nspin,nEns)
+
+  double precision,intent(in)   :: J(nBas,nBas,nspin)
+  double precision,intent(in)   :: Fx(nBas,nBas,nspin)
+  double precision,intent(in)   :: FxHF(nBas,nBas,nspin)
+  double precision,intent(in)   :: Fc(nBas,nBas,nspin)
+
+! Local variables
+
+  double precision              :: ET(nspin,nEns)
+  double precision              :: EV(nspin,nEns)
+  double precision              :: EJ(nsp,nEns)
+  double precision              :: Ex(nspin,nEns)
+  double precision              :: Ec(nsp,nEns)
+  double precision              :: EcLZ(nsp)
+  double precision              :: EcDD(nsp,nEns) 
+
+  double precision,external     :: trace_matrix
+
+  integer                       :: ispin,iEns
+
+! Output variables
+
+  double precision,intent(out)  :: E(nEns)
+  double precision,intent(out)  :: Om(nEns)
+
+!------------------------------------------------------------------------
+! Kinetic energy
+!------------------------------------------------------------------------
+
+  do ispin=1,nspin
+    do iEns=1,nEns
+      ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:)))
+    end do
+  end do
+
+!------------------------------------------------------------------------
+! Potential energy
+!------------------------------------------------------------------------
+
+  do iEns=1,nEns
+    do ispin=1,nspin
+      EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:)))
+    end do
+  end do
+
+!------------------------------------------------------------------------
+! Hartree energy
+!------------------------------------------------------------------------
+
+  do iEns=1,nEns
+
+    do ispin=1,nspin
+      call hartree_coulomb(nBas,P(:,:,ispin,iEns),ERI,J(:,:,ispin))
+    end do
+
+    EJ(1,iEns) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1)))
+    EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2)))
+    EJ(3,iEns) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2)))
+
+  end do
+
+!------------------------------------------------------------------------
+! Exchange energy
+!------------------------------------------------------------------------
+
+  do iEns=1,nEns
+    do ispin=1,nspin
+
+      call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),ERI(:,:,:,:), & 
+                              AO(:,:),dAO(:,:,:),rho(:,ispin,iEns),drho(:,:,ispin,iEns),Fx(:,:,ispin),FxHF(:,:,ispin))
+      call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),FxHF(:,:,ispin), &
+                           rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns))
+
+    end do
+  end do
+
+!------------------------------------------------------------------------
+! Correlation energy
+!------------------------------------------------------------------------
+
+  do iEns=1,nEns
+
+    call correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:), & 
+                                       rho(:,:,iEns),drho(:,:,:,iEns),Ec(:,iEns))
+
+  end do
+
+!------------------------------------------------------------------------
+! Compute Levy-Zahariev shift
+!------------------------------------------------------------------------
+
+  call correlation_Levy_Zahariev_shift(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:,:),drho(:,:,:,:),EcLZ(:))
+
+!------------------------------------------------------------------------
+! Compute derivative discontinuities
+!------------------------------------------------------------------------
+
+  call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),EcDD(:,:))
+
+!------------------------------------------------------------------------
+! Total energy
+!------------------------------------------------------------------------
+
+  do iEns=1,nEns
+    E(iEns) = ENuc + sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) &
+                   + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(EcLZ(:)) + sum(EcDD(:,iEns))
+  end do
+
+!------------------------------------------------------------------------
+! Excitation energies
+!------------------------------------------------------------------------
+
+  do iEns=1,nEns
+    Om(iEns) = E(iEns) - E(1)
+  end do
+
+!------------------------------------------------------------------------
+! Dump results
+!------------------------------------------------------------------------
+
+  call print_individual_energy(nEns,EJ,Ex,Ec,EcLZ,EcDD,E,Om)
+
+end subroutine individual_energy

src/xcDFT/lda_correlation_Levy_Zahariev_shift.f90

@@ -0,0 +1,51 @@
+subroutine lda_correlation_Levy_Zahariev_shift(DFA,nEns,wEns,nGrid,weight,rho,EcLZ)
+
+! Compute the lda correlation part of the Levy-Zahariev shift
+
+  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)   :: rho(nGrid,nspin)
+
+! Output variables
+
+  double precision,intent(out)  :: EcLZ(nsp)
+
+! Select correlation functional
+
+  select case (DFA)
+
+!   Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
+
+    case ('W38')
+
+      call W38_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:,:),EcLZ(:))
+
+!   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
+
+    case ('VWN5')
+
+      call VWN5_lda_correlation_Levy_Zahariev_shift(nGrid,weight(:),rho(:,:),EcLZ(:))
+
+!   Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation)
+
+    case ('LF19')
+
+      call LF19_lda_correlation_Levy_Zahariev_shift(nEns,wEns,nGrid,weight(:),rho(:,:),EcLZ(:))
+
+    case default
+
+      call print_warning('!!! LDA correlation functional not available !!!')
+      stop
+
+  end select
+ 
+end subroutine lda_correlation_Levy_Zahariev_shift

src/xcDFT/lda_correlation_derivative_discontinuity.f90

@@ -0,0 +1,54 @@
+subroutine lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
+
+! Compute the correlation LDA 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,nspin)
+
+! Local variables
+
+  double precision              :: aC
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp,nEns)
+
+! Select correlation functional
+
+  select case (DFA)
+
+!   Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
+
+    case ('W38')
+
+      Ec(:,:) = 0d0
+
+!   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
+
+    case ('VWN5')
+
+      Ec(:,:) = 0d0
+
+!   Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation)
+
+    case ('LF19')
+
+      call LF19_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:,:),Ec(:,:))
+
+    case default
+
+      call print_warning('!!! LDA correlation functional not available !!!')
+      stop
+
+  end select
+ 
+end subroutine lda_correlation_derivative_discontinuity

src/xcDFT/lda_correlation_energy.f90

@@ -0,0 +1,60 @@
+subroutine lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec)
+
+! Select LDA correlation functional
+
+  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)   :: rho(nGrid,nspin)
+
+! Output variables
+
+  double precision,intent(out)  :: Ec(nsp)
+
+! Select correlation functional
+
+  select case (DFA)
+
+!   Hartree-Fock
+
+    case ('HF')
+
+      Ec(:) = 0d0
+
+    case ('W38')
+
+      call W38_lda_correlation_energy(nGrid,weight(:),rho(:,:),Ec(:))
+
+!   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
+
+    case ('VWN5')
+
+      call VWN5_lda_correlation_energy(nGrid,weight(:),rho(:,:),Ec(:))
+
+!   Chachiyo's LDA correlation functional: Chachiyo, JCP 145 (2016) 021101
+
+    case ('C16')
+
+      call C16_lda_correlation_energy(nGrid,weight(:),rho(:,:),Ec(:))
+
+!   Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation)
+
+    case ('LF19')
+
+      call LF19_lda_correlation_energy(nEns,wEns(:),nGrid,weight(:),rho(:,:),Ec(:))
+
+    case default
+
+      call print_warning('!!! LDA correlation functional not available !!!')
+      stop
+
+  end select
+
+end subroutine lda_correlation_energy

src/xcDFT/lda_correlation_individual_energy.f90

@@ -0,0 +1,51 @@
+subroutine lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ec)
+
+! Compute LDA correlation 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,nspin)
+  double precision,intent(in)   :: rho(nGrid,nspin)
+
+! Output variables
+
+  double precision              :: Ec(nsp)
+
+! Select correlation functional
+
+  select case (DFA)
+
+!   Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
+
+    case ('W38')
+
+      call W38_lda_correlation_individual_energy(nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:))
+
+!   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
+
+    case ('VWN5')
+
+      call VWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:))
+
+!   Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation)
+
+    case ('LF19')
+
+      call LF19_lda_correlation_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:))
+
+    case default
+
+      call print_warning('!!! LDA correlation functional not available !!!')
+      stop
+
+  end select
+
+end subroutine lda_correlation_individual_energy

src/xcDFT/lda_correlation_potential.f90

@@ -0,0 +1,64 @@
+subroutine lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc)
+
+! Select LDA 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)   :: rho(nGrid,nspin)
+
+! Output variables
+
+  double precision,intent(out)  :: Fc(nBas,nBas,nspin)
+
+! Select correlation functional
+
+  select case (DFA)
+
+!   Hartree-Fock 
+
+    case ('HF')
+
+      Fc(:,:,:) = 0d0
+
+!   Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678
+
+    case ('W38')
+
+      call W38_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
+
+!   Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
+
+    case ('VWN5')
+
+      call VWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
+
+!   Chachiyo's LDA correlation functional: Chachiyo, JCP 145 (2016) 021101
+
+    case ('C16')
+
+      call C16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
+
+!   Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation)
+
+    case ('LF19')
+
+      call LF19_lda_correlation_potential(nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:))
+
+    case default
+
+      call print_warning('!!! LDA correlation functional not available !!!')
+      stop
+
+  end select
+
+end subroutine lda_correlation_potential

src/xcDFT/lda_exchange_energy.f90

@@ -0,0 +1,38 @@
+subroutine lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,Ex)
+
+! Select LDA exchange functional
+
+  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)   :: rho(nGrid)
+
+! Output variables
+
+  double precision,intent(out)  :: Ex
+
+! Select correlation functional
+
+  select case (DFA)
+
+!   Slater's LDA correlation functional: Slater, Phys. Rev. 81 (1951) 385
+
+    case ('S51')
+
+      call S51_lda_exchange_energy(nGrid,weight,rho,Ex)
+
+    case default
+
+      call print_warning('!!! LDA exchange functional not available !!!')
+      stop
+
+  end select
+
+end subroutine lda_exchange_energy

src/xcDFT/lda_exchange_potential.f90

@@ -0,0 +1,41 @@
+subroutine lda_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
+
+! Select LDA 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)   :: rho(nGrid)
+
+! Output variables
+
+  double precision,intent(out)  :: Fx(nBas,nBas)
+
+! Select exchange functional
+
+  select case (DFA)
+
+!   Slater's LDA correlation functional: Slater, Phys. Rev. 81 (1951) 385
+
+    case ('S51')
+
+      call S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
+
+    case default
+
+      call print_warning('!!! LDA exchange functional not available !!!')
+      stop
+
+  end select
+
+end subroutine lda_exchange_potential

src/xcDFT/one_electron_density.f90

@@ -0,0 +1,47 @@
+subroutine one_electron_density(nGrid,nBas,P,AO,dAO,rho,drho)
+
+! Calculate one-electron density
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nGrid
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: P(nBas,nBas)
+  double precision,intent(in)   :: AO(nBas,nGrid)
+  double precision,intent(in)   :: dAO(3,nBas,nGrid)
+
+! Local variables
+
+  integer                       :: ixyz,iG,mu,nu
+  double precision,external     :: trace_matrix
+
+! Output variables
+
+  double precision,intent(out)  :: rho(nGrid) 
+  double precision,intent(out)  :: drho(3,nGrid)
+
+  rho(:) = 0d0
+  do iG=1,nGrid
+    do mu=1,nBas
+      do nu=1,nBas
+        rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
+      enddo
+    enddo
+  enddo
+ 
+  drho(:,:) = 0d0
+  do ixyz=1,3
+    do iG=1,nGrid
+      do mu=1,nBas
+        do nu=1,nBas
+          drho(ixyz,iG) = drho(ixyz,iG) & 
+                        + P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG))
+        enddo
+      enddo
+    enddo
+  enddo
+
+end subroutine one_electron_density

src/xcDFT/orthogonalization_matrix.f90

@@ -0,0 +1,63 @@
+subroutine orthogonalization_matrix(nBas,S,X)
+
+! Compute the orthogonalization matrix X = S^(-1/2)
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  double precision,intent(in)   :: S(nBas,nBas)
+
+! Local variables
+
+  logical                       :: debug
+  double precision,allocatable  :: UVec(:,:),Uval(:)
+  double precision,parameter    :: thresh = 1d-6
+
+  integer                       :: i
+
+! Output variables
+
+  double precision,intent(out)  :: X(nBas,nBas)
+
+  debug = .false.
+
+  allocate(Uvec(nBas,nBas),Uval(nBas))
+
+  write(*,*)
+  write(*,*) ' *** Lowdin orthogonalization X = S^(-1/2)  *** '
+  write(*,*)
+
+  Uvec = S
+  call diagonalize_matrix(nBas,Uvec,Uval)
+
+  do i=1,nBas
+
+    if(Uval(i) > thresh) then 
+
+      Uval(i) = 1d0/sqrt(Uval(i))
+
+    else
+
+      write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization'
+
+    endif
+
+  enddo
+  
+  call ADAt(nBas,Uvec,Uval,X)
+
+! Print results
+
+  if(debug) then
+
+    write(*,'(A28)') '----------------------'
+    write(*,'(A28)') 'Orthogonalization matrix'
+    write(*,'(A28)') '----------------------'
+    call matout(nBas,nBas,X)
+    write(*,*)
+
+  endif
+
+end subroutine orthogonalization_matrix

src/xcDFT/print_KS.f90

@@ -0,0 +1,102 @@
+subroutine print_KS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew)
+
+! Print one- and two-electron energies and other stuff for KS calculation
+
+  implicit none
+  include 'parameters.h'
+
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO(nspin)
+  double precision,intent(in)        :: eps(nBas,nspin)
+  double precision,intent(in)        :: c(nBas,nBas,nspin)
+  double precision,intent(in)        :: ENuc
+  double precision,intent(in)        :: ET(nspin)
+  double precision,intent(in)        :: EV(nspin)
+  double precision,intent(in)        :: EJ(nsp)
+  double precision,intent(in)        :: Ex(nspin)
+  double precision,intent(in)        :: Ec(nsp)
+  double precision,intent(in)        :: Ew
+
+  integer                            :: HOMO(nspin)
+  integer                            :: LUMO(nspin)
+  double precision                   :: Gap(nspin)
+
+! HOMO and LUMO
+
+  HOMO(:) = nO(:)
+
+  LUMO(:) = HOMO(:) + 1
+
+  Gap(1) = eps(LUMO(1),1) - eps(HOMO(1),1)
+  Gap(2) = eps(LUMO(2),2) - eps(HOMO(2),2)
+
+! Dump results
+
+
+  write(*,*)
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,'(A40)')              ' Summary              '
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,'(A40,1X,F16.10,A3)') ' One-electron    energy: ',sum(ET(:))  + sum(EV(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' One-electron a  energy: ',ET(1) + EV(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' One-electron b  energy: ',ET(2) + EV(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kinetic         energy: ',sum(ET(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kinetic      a  energy: ',ET(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kinetic      b  energy: ',ET(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Potential       energy: ',sum(EV(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Potential    a  energy: ',EV(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Potential    b  energy: ',EV(2),' au'
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron a  energy: ',sum(EJ(:)) + sum(Ex(:))  + sum(Ec(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1) + Ec(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2) + Ec(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2) + Ec(3),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Coulomb         energy: ',sum(EJ(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Coulomb      aa energy: ',EJ(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Coulomb      ab energy: ',EJ(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Coulomb      bb energy: ',EJ(3),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Exchange        energy: ',sum(Ex(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Exchange     a  energy: ',Ex(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Exchange     b  energy: ',Ex(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Correlation     energy: ',sum(Ec(:)),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Correlation  aa energy: ',Ec(1),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Correlation  ab energy: ',Ec(2),' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Correlation  bb energy: ',Ec(3),' au'
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,'(A40,1X,F16.10,A3)') ' Electronic      energy: ',Ew,' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Nuclear      repulsion: ',ENuc,' au'
+  write(*,'(A40,1X,F16.10,A3)') ' Kohn-Sham       energy: ',Ew + ENuc,' au'
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,'(A40,F13.6,A3)')     ' KS HOMO a    energy:',eps(HOMO(1),1)*HatoeV,' eV'
+  write(*,'(A40,F13.6,A3)')     ' KS LUMO a    energy:',eps(LUMO(1),1)*HatoeV,' eV'
+  write(*,'(A40,F13.6,A3)')     ' KS HOMOa-LUMOa  gap:',Gap(1)*HatoeV,' eV'
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,'(A40,F13.6,A3)')     ' KS HOMO b    energy:',eps(HOMO(2),2)*HatoeV,' eV'
+  write(*,'(A40,F13.6,A3)')     ' KS LUMO b    energy:',eps(LUMO(2),2)*HatoeV,' eV'
+  write(*,'(A40,F13.6,A3)')     ' KS HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV'
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,*)
+
+! Print results
+
+  write(*,'(A50)') '-----------------------------------------'
+  write(*,'(A50)') 'Kohn-Sham spin-up   orbital coefficients '
+  write(*,'(A50)') '-----------------------------------------'
+  call matout(nBas,nBas,c(:,:,1))
+  write(*,'(A50)') '-----------------------------------------'
+  write(*,'(A50)') 'Kohn-Sham spin-down orbital coefficients '
+  write(*,'(A50)') '-----------------------------------------'
+  call matout(nBas,nBas,c(:,:,2))
+  write(*,*)
+  write(*,'(A50)') '---------------------------------------'
+  write(*,'(A50)') ' Kohn-Sham spin-up   orbital energies  '
+  write(*,'(A50)') '---------------------------------------'
+  call matout(nBas,1,eps(:,1))
+  write(*,*)
+  write(*,'(A50)') '---------------------------------------'
+  write(*,'(A50)') ' Kohn-Sham spin-down orbital energies  '
+  write(*,'(A50)') '---------------------------------------'
+  call matout(nBas,1,eps(:,2))
+  write(*,*)
+
+end subroutine print_KS

src/xcDFT/print_individual_energy.f90

@@ -0,0 +1,105 @@
+subroutine print_individual_energy(nEns,EJ,Ex,Ec,EcLZ,EcDD,E,Om)
+
+! Print individual energies for eDFT calculation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)                 :: nEns
+  double precision,intent(in)        :: EJ(nsp,nEns)
+  double precision,intent(in)        :: Ex(nspin,nEns)
+  double precision,intent(in)        :: Ec(nsp,nEns)
+  double precision,intent(in)        :: EcLZ(nsp)
+  double precision,intent(in)        :: EcDD(nsp,nEns)
+  double precision,intent(in)        :: E(nEns)
+  double precision,intent(in)        :: Om(nEns)
+
+! Local variables
+
+  integer                            :: iEns
+
+!------------------------------------------------------------------------
+! Hartree energy
+!------------------------------------------------------------------------
+
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,'(A50)')           ' Individual Hartree     energies'
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=1,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree     energy state ',iEns,': ',sum(EJ(:,iEns)),' au'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
+
+!------------------------------------------------------------------------
+! Exchange energy
+!------------------------------------------------------------------------
+
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,'(A50)')           ' Individual exchange    energies'
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=1,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange    energy state ',iEns,': ',sum(Ex(:,iEns)),' au'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
+
+!------------------------------------------------------------------------
+! Correlation energy
+!------------------------------------------------------------------------
+
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,'(A50)')           ' Individual correlation energies'
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=1,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
+
+!------------------------------------------------------------------------
+! Compute Levy-Zahariev shift
+!------------------------------------------------------------------------
+
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,'(A40,2X,2X,F16.10,A3)') ' Levy-Zahariev shifts: ',sum(EcLZ(:)),' au'
+  write(*,'(A60)')              '-------------------------------------------------'
+  write(*,*)
+
+!------------------------------------------------------------------------
+! Compute derivative discontinuities
+!------------------------------------------------------------------------
+
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,'(A50)')           ' Derivative discontinuities (DD)    '
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=1,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation part of DD ',iEns,': ',sum(EcDD(:,iEns)),' au'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
+
+!------------------------------------------------------------------------
+! Total and Excitation energies
+!------------------------------------------------------------------------
+
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,'(A50)')           ' Individual and excitation energies '
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=1,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns),' au'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=2,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy  1 ->',iEns,': ',Om(iEns),' au'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+  do iEns=2,nEns
+    write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy  1 ->',iEns,': ',Om(iEns)*HaToeV,' eV'
+  end do
+  write(*,'(A60)')           '-------------------------------------------------'
+  write(*,*)
+
+end subroutine print_individual_energy

src/xcDFT/quadrature_grid.f90

@@ -0,0 +1,77 @@
+subroutine quadrature_grid(nRad,nAng,nGrid,root,weight)
+
+! Build roots and weights of quadrature grid
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)                 :: nRad,nAng,nGrid
+
+! Local variables
+
+  integer                            :: i,j,k
+  double precision                   :: scale
+  double precision,allocatable       :: Radius(:)
+  double precision,allocatable       :: RadWeight(:)
+  double precision,allocatable       :: XYZ(:,:)
+  double precision,allocatable       :: XYZWeight(:)
+
+! Output variables
+
+  double precision,intent(out)       :: root(3,nGrid)
+  double precision,intent(out)       :: weight(nGrid)
+
+! Memory allocation
+
+  allocate(Radius(nRad),RadWeight(nRad),XYZ(3,nAng),XYZWeight(nAng))
+
+! Findthe radial grid
+
+  scale = 1d0
+  call EulMac(Radius,RadWeight,nRad,scale)
+
+  write(*,20)
+  write(*,30)
+  write(*,20)
+  do i = 1,nRad
+     write(*,40) i,Radius(i),RadWeight(i)
+  end do
+  write(*,20)
+  write(*,*)
+
+! Find the angular grid
+
+  call Lebdev(XYZ,XYZWeight,nAng)
+
+  write(*,20)
+  write(*,50)
+  write(*,20)
+  do j = 1,nAng
+     write(*,60) j,(XYZ(k,j),k=1,3), XYZWeight(j)
+  end do
+  write(*,20)
+
+! Form the roots and weights
+
+  k = 0
+  do i=1,nRad
+    do j=1,nAng
+       k = k + 1
+       root(:,k) = Radius(i)*XYZ(:,j)
+       weight(k) = RadWeight(i)*XYZWeight(j)
+    enddo
+  enddo
+ 
+! Compute values of the basis functions (and the its gradient if required) at each grid point
+
+20 format(T2,58('-'))
+30 format(T20,'Radial Quadrature',/, &
+          T6,'I',T26,'Radius',T50,'Weight')
+40 format(T3,I4,T18,F17.10,T35,F25.10)
+50 format(T20,'Angular Quadrature',/, &
+          T6,'I',T19,'X',T29,'Y',T39,'Z',T54,'Weight')
+60 format(T3,I4,T13,3F10.5,T50,F10.5)
+
+end subroutine quadrature_grid

src/xcDFT/read_basis.f90

@@ -0,0 +1,119 @@
+subroutine read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
+
+! Read basis set information
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nAt,nO(nspin)
+  double precision,intent(in)   :: rAt(nAt,ncart)
+
+! Local variables
+
+  integer                       :: nShAt,iAt,iShell
+  integer                       :: i,j,k
+  character                     :: shelltype
+
+! Output variables
+
+  integer,intent(out)           :: nShell,nBas
+  double precision,intent(out)  :: CenterShell(maxShell,ncart)
+  integer,intent(out)           :: TotAngMomShell(maxShell),KShell(maxShell)
+  double precision,intent(out)  :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
+  integer,intent(out)           :: nV(nspin)
+
+!------------------------------------------------------------------------
+! Primary basis set information
+!------------------------------------------------------------------------
+
+! Open file with basis set specification
+
+  open(unit=2,file='input/basis')
+
+! Read basis information
+
+  write(*,'(A28)') 'Gaussian basis set'
+  write(*,'(A28)') '------------------'
+
+  nShell = 0
+  do i=1,nAt
+    read(2,*) iAt,nShAt
+    write(*,'(A28,1X,I16)') 'Atom n. ',iAt
+    write(*,'(A28,1X,I16)') 'number of shells ',nShAt
+    write(*,'(A28)') '------------------'
+
+!   Basis function centers
+
+    do j=1,nShAt
+      nShell = nShell + 1
+      do k=1,ncart
+        CenterShell(nShell,k) = rAt(iAt,k)
+      enddo
+
+! Shell type and contraction degree
+
+      read(2,*) shelltype,KShell(nShell)
+        if(shelltype == "S") then
+          TotAngMomShell(nShell) = 0
+          write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
+        elseif(shelltype == "P") then
+          TotAngMomShell(nShell) = 1
+          write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
+        elseif(shelltype == "D") then
+          TotAngMomShell(nShell) = 2
+          write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
+        elseif(shelltype == "F") then
+          TotAngMomShell(nShell) = 3
+          write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
+        elseif(shelltype == "G") then
+          TotAngMomShell(nShell) = 4
+          write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
+        elseif(shelltype == "H") then
+          TotAngMomShell(nShell) = 5
+          write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
+        elseif(shelltype == "I") then
+          TotAngMomShell(nShell) = 6
+          write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
+        endif
+
+! Read exponents and contraction coefficients
+
+        write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
+        do k=1,Kshell(nShell)
+          read(2,*) ExpShell(nShell,k),DShell(nShell,k)
+          write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
+        enddo
+    enddo
+    write(*,'(A28)') '------------------'
+  enddo
+
+! Total number of shells
+
+  write(*,'(A28,1X,I16)') 'Number of shells',nShell
+  write(*,'(A28)') '------------------'
+  write(*,*)
+
+! Close file with basis set specification
+
+  close(unit=2)
+
+! Calculate number of basis functions
+
+  nBas = 0
+  do iShell=1,nShell
+    nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2
+  enddo
+
+  write(*,'(A28)') '------------------'
+  write(*,'(A28,1X,I16)') 'Number of basis functions',NBas
+  write(*,'(A28)') '------------------'
+  write(*,*)
+
+! Number of virtual orbitals
+
+  nV(1) = nBas - nO(1)
+  nV(2) = nBas - nO(2)
+
+end subroutine read_basis

src/xcDFT/read_geometry.f90

@@ -0,0 +1,68 @@
+subroutine read_geometry(nAt,ZNuc,rA,ENuc)
+
+! Read molecular geometry
+
+  implicit none
+
+  include 'parameters.h'
+
+! Ouput variables
+
+  integer,intent(in)            :: nAt
+
+! Local variables
+
+  integer                       :: i,j
+  double precision              :: RAB
+  character(len=2)              :: El
+  integer,external              :: element_number
+
+! Ouput variables
+
+  double precision,intent(out)  :: ZNuc(NAt),rA(nAt,ncart),ENuc
+
+! Open file with geometry specification
+
+  open(unit=1,file='input/molecule')
+
+! Read geometry
+
+  read(1,*) 
+  read(1,*) 
+  read(1,*) 
+
+  do i=1,nAt
+    read(1,*) El,rA(i,1),rA(i,2),rA(i,3)
+    ZNuc(i) = element_number(El)
+  enddo
+
+! Compute nuclear repulsion energy
+
+  ENuc = 0
+
+  do i=1,nAt-1
+    do j=i+1,nAt
+      RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2
+      ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB)
+    enddo
+  enddo
+
+! Close file with geometry specification
+  close(unit=1)
+
+! Print geometry
+  write(*,'(A28)') '------------------'
+  write(*,'(A28)') 'Molecular geometry'
+  write(*,'(A28)') '------------------'
+  do i=1,NAt
+    write(*,'(A28,1X,I16)') 'Atom n. ',i
+    write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i)
+    write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rA(i,j),j=1,ncart)
+  enddo
+  write(*,*)
+  write(*,'(A28)') '------------------'
+  write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc
+  write(*,'(A28)') '------------------'
+  write(*,*)
+
+end subroutine read_geometry

src/xcDFT/read_grid.f90

@@ -0,0 +1,47 @@
+subroutine read_grid(SGn,nRad,nAng,nGrid)
+
+! Read grid type
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: SGn
+
+! Output variables
+
+  integer,intent(out)           :: nRad
+  integer,intent(out)           :: nAng
+  integer,intent(out)           :: nGrid
+
+  write(*,*)'----------------------------------------------------------'
+  write(*,'(A22,I1)')'  Quadrature grid: SG-',SGn
+  write(*,*)'----------------------------------------------------------'
+
+  select case (SGn)
+
+    case(0)
+     nRad = 23
+     nAng = 170
+
+    case(1)
+      nRad = 50
+      nAng = 194
+
+    case(2)
+      nRad = 75
+      nAng = 302
+
+    case(3)
+      nRad = 99
+      nAng = 590
+
+    case default
+      call print_warning('!!! Quadrature grid not available !!!')
+      stop
+
+  end select
+
+  nGrid = nRad*nAng
+
+end subroutine read_grid

src/xcDFT/read_integrals.f90

@@ -0,0 +1,114 @@
+subroutine read_integrals(nBas,S,T,V,Hc,G)
+
+! Read one- and two-electron integrals from files
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+
+! Local variables
+
+  logical                       :: debug
+  integer                       :: mu,nu,la,si
+  double precision              :: Ov,Kin,Nuc,ERI
+
+! Output variables
+
+  double precision,intent(out)  :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas)
+
+! Open file with integrals
+
+  debug = .false.
+
+  open(unit=8 ,file='int/Ov.dat')
+  open(unit=9 ,file='int/Kin.dat')
+  open(unit=10,file='int/Nuc.dat')
+  open(unit=11,file='int/ERI.dat')
+
+! Read overlap integrals
+
+  S = 0d0
+  do 
+    read(8,*,end=8) mu,nu,Ov
+    S(mu,nu) = Ov
+  enddo
+  8 close(unit=8)
+
+! Read kinetic integrals
+
+  T = 0d0
+  do 
+    read(9,*,end=9) mu,nu,Kin
+    T(mu,nu) = Kin
+  enddo
+  9 close(unit=9)
+
+! Read nuclear integrals
+
+  V = 0d0
+  do 
+    read(10,*,end=10) mu,nu,Nuc
+    V(mu,nu) = Nuc
+  enddo
+  10 close(unit=10)
+
+! Define core Hamiltonian
+
+  Hc = T + V
+
+! Read nuclear integrals
+
+  G = 0d0
+  do 
+    read(11,*,end=11) mu,nu,la,si,ERI
+!   (12|34)
+    G(mu,nu,la,si) = ERI
+!   (21|34)
+    G(nu,mu,la,si) = ERI
+!   (12|43)
+    G(mu,nu,si,la) = ERI
+!   (21|43)
+    G(nu,mu,si,la) = ERI
+!   (34|12)
+    G(la,si,mu,nu) = ERI
+!   (43|12)
+    G(si,la,mu,nu) = ERI
+!   (34|21)
+    G(la,si,nu,mu) = ERI
+!   (43|21)
+    G(si,la,nu,mu) = ERI
+  enddo
+  11 close(unit=11)
+
+
+! Print results
+  if(debug) then
+    write(*,'(A28)') '----------------------'
+    write(*,'(A28)') 'Overlap integrals'
+    write(*,'(A28)') '----------------------'
+    call matout(nBas,nBas,S)
+    write(*,*)
+    write(*,'(A28)') '----------------------'
+    write(*,'(A28)') 'Kinetic integrals'
+    write(*,'(A28)') '----------------------'
+    call matout(nBas,nBas,T)
+    write(*,*)
+    write(*,'(A28)') '----------------------'
+    write(*,'(A28)') 'Nuclear integrals'
+    write(*,'(A28)') '----------------------'
+    call matout(nBas,nBas,V)
+    write(*,*)
+    write(*,'(A28)') '----------------------'
+    write(*,'(A28)') 'Electron repulsion integrals'
+    write(*,'(A28)') '----------------------'
+    do la=1,nBas
+      do si=1,nBas
+        call matout(nBas,nBas,G(1,1,la,si))
+      enddo
+    enddo
+    write(*,*)
+  endif
+
+end subroutine read_integrals

src/xcDFT/read_molecule.f90

@@ -0,0 +1,52 @@
+subroutine read_molecule(nAt,nEl,nO)
+
+! Read number of atoms nAt and number of electrons nEl
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(out)           :: nAt,nEl(nspin),nO(nspin)
+
+! Local variables
+
+  integer                       :: n
+
+! Open file with geometry specification
+
+  open(unit=1,file='input/molecule')
+
+! Read number of atoms and number of electrons
+
+  read(1,*) 
+  read(1,*) nAt,n,nEl(1),nEl(2)
+
+! Check imput consistency
+
+  if(n /= sum(nEl(:))) then
+    write(*,*) 'number of electrons inconsistent'
+    stop
+  endif
+
+  nO(:) = nEl(:)
+
+! Print results
+
+  write(*,'(A28)') '----------------------'
+  write(*,'(A28,1X,I16)') 'Number of atoms',nAt
+  write(*,'(A28)') '----------------------'
+  write(*,*)
+  write(*,'(A28)') '----------------------'
+  write(*,'(A28,1X,I16)') 'Number of spin-up   electrons',nEl(1)
+  write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nEl(2)
+  write(*,'(A28,1X,I16)') 'Total number of     electrons',sum(nEl(:))
+  write(*,'(A28)') '----------------------'
+  write(*,*)
+
+! Close file with geometry specification
+
+  close(unit=1)
+
+end subroutine read_molecule

src/xcDFT/read_options.f90

@@ -0,0 +1,106 @@
+subroutine read_options(x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type)
+
+! Read DFT options
+
+  implicit none
+
+  include 'parameters.h'
+
+! Local variables
+
+  integer                       :: I
+
+! Output variables
+
+  integer,intent(out)           :: x_rung,c_rung
+  character(len=12),intent(out) :: x_DFA, c_DFA
+  integer,intent(out)           :: SGn
+  integer,intent(out)           :: nEns
+  double precision,intent(out)  :: wEns(maxEns)
+
+  integer,intent(out)           :: maxSCF
+  double precision,intent(out)  :: thresh
+  logical,intent(out)           :: DIIS
+  integer,intent(out)           :: max_diis
+  integer,intent(out)           :: guess_type
+  integer,intent(out)           :: ortho_type
+
+! Local variables
+
+  character(len=1)              :: answer
+
+! Open file with method specification
+
+  open(unit=1,file='input/options')
+
+! Default values
+
+  x_rung  = 1
+  c_rung  = 1
+  x_DFA   = 'S51'
+  c_DFA   = 'W38'
+  SGn     = 0
+  wEns(:) = 0d0
+
+! EXCHANGE: read rung of Jacob's ladder
+
+  read(1,*)
+  read(1,*) x_rung,x_DFA
+
+! CORRELATION: read rung of Jacob's ladder
+
+  read(1,*)
+  read(1,*) c_rung,c_DFA
+
+! Read SG-n grid
+
+  read(1,*)
+  read(1,*) SGn
+
+! Read number of states in ensemble
+
+  read(1,*)
+  read(1,*) nEns
+
+  if(nEns.gt.maxEns) then
+    write(*,*) ' Number of states in ensemble too big!! ' 
+    stop
+  endif
+
+  write(*,*)'----------------------------------------------------------'
+  write(*,'(A33,I3)')'  Number of states in ensemble = ',nEns
+  write(*,*)'----------------------------------------------------------'
+  write(*,*) 
+  
+! Read ensemble weights
+  read(1,*)
+  read(1,*) (wEns(I),I=2,nEns)
+  wEns(1) = 1d0 - sum(wEns)
+
+  write(*,*)'----------------------------------------------------------'
+  write(*,*)' Ensemble weights '
+  write(*,*)'----------------------------------------------------------'
+  call matout(nEns,1,wEns)
+  write(*,*) 
+
+! Read KS options
+
+  maxSCF     = 64
+  thresh     = 1d-6
+  DIIS       = .false.
+  max_diis   = 5
+  guess_type = 1
+  ortho_type = 1
+
+  read(1,*)
+  read(1,*) maxSCF,thresh,answer,max_diis,guess_type,ortho_type
+
+  if(answer == 'T') DIIS = .true.
+
+  if(.not.DIIS) max_diis = 1
+
+! Close file with options
+
+  close(unit=1)
+
+end subroutine read_options

src/xcDFT/select_rung.f90

@@ -0,0 +1,53 @@
+subroutine select_rung(rung,DFA)
+
+! Select rung of Jacob's ladder
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: rung
+  character(len=12),intent(in)  :: DFA
+
+  select case (rung)
+
+! Hartree calculation
+    case(0)
+      write(*,*) "* 0th rung of Jacob's ladder: Hartree calculation *"
+
+! LDA functionals
+    case(1)
+      write(*,*) "* 1st rung of Jacob's ladder: local-density approximation (LDA)   *"
+
+! GGA functionals
+    case(2)
+      write(*,*) "* 2nd rung of Jacob's ladder: generalized gradient approximation (GGA) *"
+
+! meta-GGA functionals
+    case(3)
+      write(*,*) "* 3rd rung of Jacob's ladder: meta-GGA functional (MGGA) *"
+
+! Hybrid functionals
+    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 !!!"
+      stop
+
+  end select
+
+  ! Print selected functional
+
+  write(*,*) '* You have selected the following functional: ',DFA,'        *'
+  write(*,*) '*******************************************************************'
+  write(*,*)
+
+
+end subroutine select_rung

src/xcDFT/utils.f90

@@ -0,0 +1,400 @@
+!------------------------------------------------------------------------
+function Kronecker_delta(i,j) result(delta)
+
+! Kronecker Delta
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: i,j
+
+! Output variables
+
+  double precision              :: delta
+
+  if(i == j) then
+    delta = 1d0
+  else
+    delta = 0d0
+  endif
+
+end function Kronecker_delta
+
+!------------------------------------------------------------------------
+subroutine matout(m,n,A)
+
+! Print the MxN array A
+
+  implicit none
+
+  integer,parameter             :: ncol = 5
+  double precision,parameter    :: small = 1d-10
+  integer,intent(in)            :: m,n
+  double precision,intent(in)   :: A(m,n)
+  double precision              :: B(ncol)
+  integer                       :: ilower,iupper,num,i,j
+  
+  do ilower=1,n,ncol
+    iupper = min(ilower + ncol - 1,n)
+    num = iupper - ilower + 1
+    write(*,'(3X,10(9X,I6))') (j,j=ilower,iupper)
+    do i=1,m
+      do j=ilower,iupper
+        B(j-ilower+1) = A(i,j)
+      enddo
+      do j=1,num
+        if(abs(B(j)) < small) B(j) = 0d0
+      enddo
+      write(*,'(I7,10F15.8)') i,(B(j),j=1,num)
+    enddo
+  enddo
+
+end subroutine matout
+
+!------------------------------------------------------------------------
+subroutine trace_vector(n,v,Tr)
+
+! Calculate the trace of the vector v of length n
+!!! Please use the intrinsic fortran sum()  !!!
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: n
+  double precision,intent(in)   :: v(n)
+
+! Local variables
+
+  integer                       :: i
+
+! Output variables
+
+  double precision,intent(out)  :: Tr
+
+  Tr = 0d0
+  do i=1,n
+    Tr = Tr + v(i)
+  enddo
+
+end subroutine trace_vector
+
+!------------------------------------------------------------------------
+function trace_matrix(n,A) result(Tr)
+
+! Calculate the trace of the square matrix A
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: n
+  double precision,intent(in)   :: A(n,n)
+
+! Local variables
+
+  integer                       :: i
+
+! Output variables
+
+  double precision              :: Tr
+
+  Tr = 0d0
+  do i=1,n
+    Tr = Tr + A(i,i)
+  enddo
+
+end function trace_matrix
+
+!------------------------------------------------------------------------
+subroutine compute_error(nData,Mean,Var,Error)
+
+! Calculate the statistical error
+
+  implicit none
+
+! Input variables
+
+  double precision,intent(in)   :: nData,Mean(3)
+
+! Output variables
+
+  double precision,intent(out)  :: Error(3)
+  double precision,intent(inout):: Var(3)
+  
+  Error = sqrt((Var-Mean**2/nData)/nData/(nData-1d0))
+
+end subroutine compute_error
+
+!------------------------------------------------------------------------
+subroutine identity_matrix(N,A)
+
+! Set the matrix A to the identity matrix
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: N
+
+! Local viaruabkes
+
+  integer                       :: i
+
+! Output variables
+
+  double precision,intent(out)  :: A(N,N)
+  
+  A = 0d0
+
+  do i=1,N
+    A(i,i) = 1d0
+  enddo
+     
+end subroutine identity_matrix
+
+!------------------------------------------------------------------------
+subroutine prepend(N,M,A,b)
+
+! Prepend the vector b of size N into the matrix A of size NxM
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: N,M
+  double precision,intent(in)   :: b(N)
+
+! Local viaruabkes
+
+  integer                       :: i,j
+
+! Output variables
+
+  double precision,intent(out)  :: A(N,M)
+
+
+! print*,'b in append'
+! call matout(N,1,b)
+
+  do i=1,N
+    do j=M-1,1,-1
+      A(i,j+1) = A(i,j)
+    enddo
+    A(i,1) = b(i)
+  enddo
+
+end subroutine prepend
+
+!------------------------------------------------------------------------
+subroutine append(N,M,A,b)
+
+! Append the vector b of size N into the matrix A of size NxM
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: N,M
+  double precision,intent(in)   :: b(N)
+
+! Local viaruabkes
+
+  integer                       :: i,j
+
+! Output variables
+
+  double precision,intent(out)  :: A(N,M)
+
+  do i=1,N
+    do j=2,M
+      A(i,j-1) = A(i,j)
+    enddo
+    A(i,M) = b(i)
+  enddo
+
+end subroutine append
+
+!------------------------------------------------------------------------
+subroutine AtDA(N,A,D,B)
+
+! Perform B = At.D.A where A is a NxN matrix and D is a diagonal matrix given 
+! as a vector of length N
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: N
+  double precision,intent(in)   :: A(N,N),D(N)
+
+! Local viaruabkes
+
+  integer                       :: i,j,k
+
+! Output variables
+
+  double precision,intent(out)  :: B(N,N)
+
+  B = 0d0
+
+  do i=1,N
+    do j=1,N
+      do k=1,N
+        B(i,k) = B(i,k) + A(j,i)*D(j)*A(j,k)
+      enddo
+    enddo
+  enddo
+
+end subroutine AtDA
+
+!------------------------------------------------------------------------
+subroutine ADAt(N,A,D,B)
+
+! Perform B = A.D.At where A is a NxN matrix and D is a diagonal matrix given 
+! as a vector of length N
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: N
+  double precision,intent(in)   :: A(N,N),D(N)
+
+! Local viaruabkes
+
+  integer                       :: i,j,k
+
+! Output variables
+
+  double precision,intent(out)  :: B(N,N)
+
+  B = 0d0
+
+  do i=1,N
+    do j=1,N
+      do k=1,N
+        B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j)
+      enddo
+    enddo
+  enddo
+
+end subroutine ADAt
+!------------------------------------------------------------------------
+subroutine DA(N,D,A)
+
+! Perform A <- D.A where A is a NxN matrix and D is a diagonal matrix given 
+! as a vector of length N
+
+  implicit none
+
+  integer,intent(in)            :: N
+  integer                       :: i,j,k
+  double precision,intent(in)   :: D(N)
+  double precision,intent(inout):: A(N,N)
+
+  do i=1,N
+    do j=1,N
+      A(i,j) = D(i)*A(i,j)
+    enddo
+  enddo
+
+end subroutine DA
+
+!------------------------------------------------------------------------
+subroutine AD(N,A,D)
+
+! Perform A <- A.D where A is a NxN matrix and D is a diagonal matrix given 
+! as a vector of length N
+
+  implicit none
+
+  integer,intent(in)            :: N
+  integer                       :: i,j,k
+  double precision,intent(in)   :: D(N)
+  double precision,intent(inout):: A(N,N)
+
+  do i=1,N
+    do j=1,N
+      A(i,j) = A(i,j)*D(j)
+    enddo
+  enddo
+
+end subroutine AD
+
+!------------------------------------------------------------------------
+
+subroutine print_warning(message)
+
+! Print warning
+
+  implicit none
+
+  character(len=*),intent(in)            :: message
+
+  write(*,*) message
+
+end subroutine print_warning
+
+!------------------------------------------------------------------------
+recursive function fac(n) result(fact)
+
+  implicit none
+  integer :: fact
+  integer, intent(in) :: n
+
+  if (n == 0) then
+     fact = 1
+  else
+     fact = n * fac(n-1)
+  end if
+
+end function fac
+
+!------------------------------------------------------------------------
+
+function dfac(n) result(fact)
+
+  implicit none
+  double precision :: fact
+  integer, intent(in) :: n
+  integer             :: fac
+
+  fact = dble(fac(n))
+
+end function dfac
+
+!------------------------------------------------------------------------
+
+function NormCoeff(alpha,a)
+  
+! Compute normalization coefficients for cartesian gaussians
+
+  implicit none
+
+! Input variables
+
+  double precision,intent(in)   :: alpha
+  integer,intent(in)            :: a(3)
+
+! local variable
+  double precision              :: pi,dfa(3),dfac
+  integer                       :: atot
+
+! Output variable
+  double precision NormCoeff
+
+  pi = 4d0*atan(1d0)
+  atot = a(1) + a(2) + a(3)
+
+  dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1)))
+  dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2)))
+  dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3)))
+
+
+  NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot
+  NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3))
+  NormCoeff = sqrt(NormCoeff)
+
+end function NormCoeff

src/xcDFT/wrap_lapack.f90

@@ -0,0 +1,207 @@
+!subroutine eigenvalues_non_symmetric_matrix(N,A,e)
+!
+!! Diagonalize a square matrix
+!
+!  implicit none
+!
+!! Input variables
+!
+!  integer,intent(in)            :: N
+!  double precision,intent(inout):: A(N,N)
+!  double precision,intent(out)  :: e(N)
+!
+!! Local variables
+!
+!  integer                       :: lwork,info
+!  double precision,allocatable  :: work(:)
+!
+!! Memory allocation
+!
+!  allocate(eRe(N),eIm(N),work(3*N))
+!  lwork = size(work)
+!
+!  call DGEEV('N','N',N,A,N, eRe, eIm, 0d0,1, VR,LDVR, WORK, LWORK, INFO )
+!
+!  if(info /= 0) then 
+!    print*,'Problem in diagonalize_matrix (dseev)!!'
+!    stop
+!  endif
+!
+!end subroutine eigenvalues_non_symmetric_matrix
+
+subroutine diagonalize_matrix(N,A,e)
+
+! Diagonalize a square matrix
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: N
+  double precision,intent(inout):: A(N,N)
+  double precision,intent(out)  :: e(N)
+
+! Local variables
+
+  integer                       :: lwork,info
+  double precision,allocatable  :: work(:)
+
+! Memory allocation
+
+  allocate(work(3*N))
+  lwork = size(work)
+
+  call dsyev('V','U',N,A,N,e,work,lwork,info)
+ 
+  if(info /= 0) then 
+    print*,'Problem in diagonalize_matrix (dsyev)!!'
+  endif
+
+end subroutine diagonalize_matrix
+
+subroutine svd(N,A,U,D,Vt)
+
+  ! Compute A = U.D.Vt
+  ! Dimension of A is NxN
+
+  implicit none
+
+  integer, intent(in)             :: N
+  double precision,intent(in)     :: A(N,N)
+  double precision,intent(out)    :: U(N,N)
+  double precision,intent(out)    :: Vt(N,N)
+  double precision,intent(out)    :: D(N)
+  double precision,allocatable    :: work(:)
+  integer                         :: info,lwork
+
+  double precision,allocatable    :: scr(:,:)
+
+  allocate (scr(N,N))
+
+  scr(:,:) = A(:,:)
+
+  ! Find optimal size for temporary arrays
+
+  allocate(work(1))
+
+  lwork = -1
+  call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info)
+  lwork = int(work(1))
+
+  deallocate(work)
+
+  allocate(work(lwork))
+
+  call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info)
+
+  deallocate(work,scr)
+
+  if (info /= 0) then
+    print *,  info, ': SVD failed'
+    stop
+  endif
+
+end
+
+subroutine inverse_matrix(N,A,B)
+
+! Returns the inverse of the square matrix A in B
+
+  implicit none
+
+  integer,intent(in)             :: N
+  double precision, intent(in)   :: A(N,N)
+  double precision, intent(out)  :: B(N,N)
+
+  integer                        :: info,lwork
+  integer, allocatable           :: ipiv(:)
+  double precision,allocatable   :: work(:)
+
+  allocate (ipiv(N),work(N*N))
+  lwork = size(work)
+
+  B(1:N,1:N) = A(1:N,1:N)
+
+  call dgetrf(N,N,B,N,ipiv,info)
+
+  if (info /= 0) then
+
+    print*,info
+    stop 'error in inverse (dgetrf)!!'
+
+  endif
+
+  call dgetri(N,B,N,ipiv,work,lwork,info)
+
+  if (info /= 0) then
+
+    print *,  info
+    stop 'error in inverse (dgetri)!!'
+
+  endif
+
+  deallocate(ipiv,work)
+
+end subroutine inverse_matrix
+
+subroutine linear_solve(N,A,b,x,rcond)
+
+! Solve the linear system A.x = b where A is a NxN matrix
+! and x and x are vectors of size N
+
+  implicit none
+
+  integer,intent(in)             :: N
+  double precision,intent(in)    :: A(N,N),b(N),rcond
+  double precision,intent(out)   :: x(N)
+
+  integer                        :: info,lwork
+  double precision               :: ferr,berr
+  integer,allocatable            :: ipiv(:),iwork(:)
+  double precision,allocatable   :: AF(:,:),work(:)
+
+  lwork = 3*N
+  allocate(AF(N,N),ipiv(N),work(lwork),iwork(N))
+
+  call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info)
+
+! if (info /= 0) then
+
+!   print *,  info
+!   stop 'error in linear_solve (dsysvx)!!'
+
+! endif
+
+end subroutine linear_solve
+
+subroutine easy_linear_solve(N,A,b,x)
+
+! Solve the linear system A.x = b where A is a NxN matrix
+! and x and x are vectors of size N
+
+  implicit none
+
+  integer,intent(in)             :: N
+  double precision,intent(in)    :: A(N,N),b(N)
+  double precision,intent(out)   :: x(N)
+
+  integer                        :: info,lwork
+  integer,allocatable            :: ipiv(:)
+  double precision,allocatable   :: work(:)
+
+  allocate(ipiv(N),work(N*N))
+  lwork = size(work)
+
+  x = b
+
+  call dsysv('U',N,1,A,N,ipiv,x,N,work,lwork,info)
+
+  if (info /= 0) then
+
+    print *,  info
+    stop 'error in linear_solve (dsysv)!!'
+
+  endif
+
+end subroutine easy_linear_solve
+