Merge branch 'master' of https://github.com/pfloos/QuAcK

input/dft

@@ -6,24 +6,24 @@
 # GGA          = 2: B88,G96,PBE
 # MGGA         = 3:
 # Hybrid       = 4 HF,B3LYP,PBE
-1 S51
+1 S51  
 # correlation rung: 
 # Hartree      = 0: H
 # LDA          = 1: PW92,VWN3,VWN5,eVWN5
 # GGA          = 2: LYP,PBE
 # MGGA         = 3:
 # Hybrid       = 4: HF,B3LYP,PBE 
-1 VWN5
+1 VWN5 
 # quadrature grid SG-n
   1
 # Number of states in ensemble (nEns)
-2
+4
 # occupation numbers
-1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
-1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
  
-1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
  
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
@@ -36,8 +36,8 @@
 F
 # Parameters for CC weight-dependent exchange functional
 4
-0.642674 -0.07818 -0.0280307 0.00144198
+-0.718713,-0.133321,0.226288,-0.250718
-0.254939 -0.0893405 0.00765581 0.
+-0.525899,0.687216,-0.13866,-0.0226579
 0.0 0.0 0.0 0.0
 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2
 1

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
- 
+
 H 0. 0. 0.
 H 0. 0. 0.741

src/eDFT/CC_lda_exchange_derivative_discontinuity.f90

@@ -1,4 +1,6 @@
-subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,ExDD)
+subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,&
+                                                    doNcentered,kappa,ExDD)
+
 
 ! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative
 
@@ -16,6 +18,7 @@ subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weig
   double precision,intent(in)   :: rhow(nGrid)
   integer,intent(in)            :: Cx_choice
   logical,intent(in)            :: doNcentered
+  double precision,intent(in)   :: kappa(nEns)
 
 ! Local variables
 
@@ -153,7 +156,13 @@ subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weig
   do iEns=1,nEns
     do jEns=2,nEns
 
-        ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns)
+      if(doNcentered) then
+
+         ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - kappa(iEns)*wEns(jEns))*dExdw(jEns)
+      else
+
+         ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns)
+      end if
 
     end do
   end do

src/eDFT/UKS.f90

@@ -352,6 +352,9 @@ subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,t
   end do
   write(*,*)'------------------------------------------------------------------------------------------'
 
+!  print*,'Ensemble energy:',Ew + ENuc,'au'
+  
+
 !------------------------------------------------------------------------
 ! End of SCF loop
 !------------------------------------------------------------------------

src/eDFT/exchange_derivative_discontinuity.f90

@@ -1,5 +1,5 @@
 subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,drhow,&
-                                                          Cx_choice,doNcentered,ExDD)
+                                                          Cx_choice,doNcentered,kappa,ExDD)
 
 ! Compute the exchange part of the derivative discontinuity
 
@@ -20,11 +20,12 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,we
   double precision,intent(in)   :: drhow(ncart,nGrid)
   integer,intent(in)            :: Cx_choice
   logical,intent(in)            :: doNcentered
+  double precision,intent(in)   :: kappa(nEns)
 
-! Local variables
+!Local variables
 
 
-! Output variables
+!Output variables
 
   double precision,intent(out)  :: ExDD(nEns)
 
@@ -41,7 +42,7 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,we
     case(1) 
 
       call lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),&
-                                                              rhow(:),Cx_choice,doNcentered,ExDD(:))
+                                                              rhow(:),Cx_choice,doNcentered,kappa,ExDD(:))
 !   GGA functionals
 
     case(2) 

src/eDFT/exchange_energy.f90

@@ -49,7 +49,7 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,
 
     case(2) 
 
-      call gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
+      call gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex)
 
 !   MGGA functionals
 

src/eDFT/exchange_potential.f90

@@ -59,7 +59,8 @@ subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weig
 
     case(2) 
 
-      call gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+      call gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,&
+                                  Cx_choice,Fx)
 
 !   MGGA functionals
 

src/eDFT/gga_exchange_energy.f90

@@ -1,4 +1,4 @@
-subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
+subroutine gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex)
 
 ! Select GGA exchange functional for energy calculation
 
@@ -11,11 +11,15 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
   integer,intent(in)            :: DFA
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   double precision,intent(in)   :: rho(nGrid)
+  integer,intent(in)            :: Cx_choice
   double precision,intent(in)   :: drho(ncart,nGrid)
 
+
 ! Output variables
 
   double precision              :: Ex
@@ -34,6 +38,11 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
 
       call PBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
 
+    case (4)
+
+      call CC_B88_gga_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,&
+                                      Cx_choice,Ex)
+
     case default
 
       call print_warning('!!! GGA exchange energy not available !!!')

src/eDFT/gga_exchange_potential.f90

@@ -1,4 +1,5 @@
-subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
+subroutine gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,&
+                                  rho,drho,Cx_choice,Fx)
 
 ! Select GGA exchange functional for potential calculation
 
@@ -10,6 +11,8 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh
   integer,intent(in)            :: DFA
   integer,intent(in)            :: nEns
   double precision,intent(in)   :: wEns(nEns)
+  integer,intent(in)            :: nCC
+  double precision,intent(in)   :: aCC(nCC,nEns-1)
   integer,intent(in)            :: nGrid
   double precision,intent(in)   :: weight(nGrid)
   integer,intent(in)            :: nBas
@@ -17,6 +20,7 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh
   double precision,intent(in)   :: dAO(3,nBas,nGrid)
   double precision,intent(in)   :: rho(nGrid)
   double precision,intent(in)   :: drho(3,nGrid)
+  integer,intent(in)            :: Cx_choice
 
 ! Output variables
 
@@ -38,6 +42,11 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh
     
       call PBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
 
+    case (4)
+
+      call CC_B88_gga_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,&
+                                         Cx_choice,Fx)
+
     case default
 
       call print_warning('!!! GGA exchange potential not available !!!')

src/eDFT/individual_energy.f90

@@ -150,21 +150,22 @@ subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nC
 
   do ispin=1,nspin 
     call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,nCC,aCC,nGrid,weight, &
-                                                        rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,ExDD(ispin,:))
+                                           rhow(:,ispin),drhow(:,:,ispin),Cx_choice,&
+                                           doNcentered,kappa,ExDD(ispin,:))
   end do
 
   call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD)
 
 ! Scaling derivative discontinuity for N-centered ensembles
 
-  if(doNcentered) then
+!  if(doNcentered) then
 
-     do iEns=1,nEns
+!     do iEns=1,nEns
-       ExDD(:,iEns) = (1d0 - kappa(iEns))*ExDD(:,iEns)
+!       ExDD(:,iEns) = (1d0 - kappa(iEns))*ExDD(:,iEns)
-       EcDD(:,iEns) = (1d0 - kappa(iEns))*EcDD(:,iEns)
+!       EcDD(:,iEns) = (1d0 - kappa(iEns))*EcDD(:,iEns)
-     end do
+!     end do
 
-  end if
+!  end if
 
 !------------------------------------------------------------------------
 ! Total energy
@@ -187,6 +188,8 @@ subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nC
     end do
 
   end if
+
+print*,'LZshift=',sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:))
   
 ! do iEns=1,nEns
 !   E(iEns)   = sum(ET(:,iEns)) + sum(EV(:,iEns))                     & 

src/eDFT/lda_exchange_derivative_discontinuity.f90

@@ -1,5 +1,5 @@
 subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,&
-                                                              Cx_choice,doNcentered,ExDD)
+                                                        Cx_choice,doNcentered,kappa,ExDD)
 
 ! Compute the exchange LDA part of the derivative discontinuity
 
@@ -19,6 +19,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,wei
   double precision,intent(in)   :: rhow(nGrid)
   integer,intent(in)            :: Cx_choice
   logical,intent(in)            :: doNcentered
+  double precision,intent(in)   :: kappa(nEns)
 
 ! Local variables
 
@@ -38,7 +39,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,wei
     case (2)
 
       call CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,&
-                                                     Cx_choice,doNcentered,ExDD)
+                                                     Cx_choice,doNcentered,kappa,ExDD)
 
     case default
 

src/eDFT/read_options_dft.f90

@@ -117,7 +117,11 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,
         case ('PBE')
   
           x_DFA = 3
-  
+ 
+        case ('CC-B88')
+
+          x_DFA = 4
+
         case default
  
           call print_warning('!!! GGA exchange functional not available !!!')