Merge branch 'rel1_57_15-local' into hse-dev-local

change HSE implementation to be more similar to PBE0 in the
general kpoint case

Conflicts:
	src/ExchangeOperator.C
	src/XMLGFPreprocessor.C
	src/notes
	src/release.C

git-svn-id: http://qboxcode.org/svn/qb/branches/hse-dev@1413 cba15fb0-1239-40c8-b417-11db7ca47a34

src/ExchangeOperator.C

@@ -103,8 +103,6 @@ ExchangeOperator::ExchangeOperator( Sample& s, double HFCoeff, const Interaction
     // allocate memory for |q+G| and related quantities
     qpG21_.resize(ngloc);
     qpG22_.resize(ngloc);
-    qpG2i1_.resize(ngloc);
-    qpG2i2_.resize(ngloc);
     int_pot1_.resize(ngloc);
     int_pot2_.resize(ngloc);
   }
@@ -547,26 +545,26 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
                        vbasis_->gx(ig+ngloc*2));
 
             // compute G+q for each qi and find the value of the
-            // correction term: sum_(G,q) exp(-rcut_^2*|G+q|^2)/|G+q|^2
-            // => compute the square norm and inverse of q1+G
+            // correction term: sum_(G,q) exp(-rcut_^2*|G+q|^2) V(|G+q|)
+            // => compute the square norm of q1+G
             qpG21_[ig]  = ( G + q1 ) * ( G + q1 );
-            qpG2i1_[ig] = ( qpG21_[ig] > 0.0 ) ? 1.0 / qpG21_[ig] : 0.0;
 
-            // => compute the square norm and inverse of q2+G
+            // => compute the square norm of q2+G
             qpG22_[ig]  = ( G + q2 ) * ( G + q2 );
-            qpG2i2_[ig] = ( qpG22_[ig] > 0.0 ) ? 1.0 / qpG22_[ig] : 0.0;
 
-            // for Coulomb potential use 1/|G+q|^2
+            // for Coulomb potential V(|G+q|) = 1/|G+q|^2
             if ( coulomb_ )
             {
-              int_pot1_[ig] = qpG2i1_[ig];
-              int_pot2_[ig] = qpG2i2_[ig];
+              int_pot1_[ig] = ( qpG21_[ig] > 0.0 ) ? 1.0 / qpG21_[ig] : 0.0;
+              int_pot2_[ig] = ( qpG22_[ig] > 0.0 ) ? 1.0 / qpG22_[ig] : 0.0;
             }
-            // otherwise use given function of |G+q|^2
+            // otherwise use given function V(|G+q|^2)
             else
             {
-              int_pot1_[ig] = interaction_potential_(qpG21_[ig]);
-              int_pot2_[ig] = interaction_potential_(qpG22_[ig]);
+              int_pot1_[ig] = ( qpG21_[ig] > 0.0 ) ?
+                  interaction_potential_(qpG21_[ig]) : 0;
+              int_pot2_[ig] = ( qpG22_[ig] > 0.0 ) ?
+                  interaction_potential_(qpG22_[ig]) : 0;
             }
 
             // if iKpi=0 (first k point)
@@ -653,24 +651,22 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
                   // Add the values of |rho1(G)|^2/|G+q1|^2
                   // and |rho2(G)|^2/|G+q2|^2 to the exchange energy.
                   // This does not take the point G=q=0 into account
-                  // as qpG2i = 0.
+                  // where int_pot is set to 0.
                   const double t1 = norm(rhog1_[ig]) * int_pot1_[ig];
                   const double t2 = norm(rhog2_[ig]) * int_pot2_[ig];
-                  if ( qpG21_[ig] > 0.0 )
-                    ex_ki_i_kj_j += t1;
-                  if ( qpG22_[ig] > 0.0 )
-                    ex_ki_i_kj_j += t2;
+                  ex_ki_i_kj_j += t1;
+                  ex_ki_i_kj_j += t2;
 
                   if ( dwf )
                   {
-                    // compute rhog1_[G]/|G+q1|^2 and rhog2_[G]/|G+q1|^2
+                    // compute rhog1_[G]*V(|G+q1|) and rhog2_[G]*V(|G+q1|)
                     rhog1_[ig] *= int_pot1_[ig];
                     rhog2_[ig] *= int_pot2_[ig];
                   }
                 }
                 if ( dwf )
                 {
-                  // Backtransform rhog[G]/|q+G|^2
+                  // Backtransform rhog[G]*V(|q+G|)
                   vft_->backward(&rhog1_[0], &rhor1_[0]);
                   vft_->backward(&rhog2_[0], &rhor2_[0]);
                 }
@@ -865,9 +861,9 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
       {
         double div_corr = 0.0;
 
-        const double div_corr_1 = 2.0 * exfac * numerical_correction[iKpi];
+        const double div_corr_1 = exfac * numerical_correction[iKpi] * occ_ki_[i];
         div_corr += div_corr_1;
-        const double e_div_corr_1 = - 0.5 * div_corr_1 * occ_ki_[i];
+        const double e_div_corr_1 = -div_corr_1;
         exchange_sum += e_div_corr_1 * wf.weight(iKpi);
         // add here contributions to stress from div_corr_1;
 
@@ -937,8 +933,7 @@ double ExchangeOperator::compute_exchange_for_general_case_( Sample* s,
         } // if quad_correction
 
         // contribution of divergence corrections to forces on wave functions
-        // wave functions treated explicitly for screened Coulomb potential
-        if (dwf) // and coulomb_)
+        if (dwf)
         {
           // sum the partial contributions to the correction for state i
           gcontext_.dsum('C', 1, 1, &div_corr, 1);

src/ExchangeOperator.h

@@ -80,8 +80,6 @@ class ExchangeOperator
   // G vectors
   valarray<double> qpG21_;
   valarray<double> qpG22_;
-  valarray<double> qpG2i1_;
-  valarray<double> qpG2i2_;
   valarray<double> int_pot1_;
   valarray<double> int_pot2_;
 

src/Makefile

@@ -33,6 +33,7 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
         XCOperator.o ExchangeOperator.o Bisection.o KPGridConnectivity.o \
 	XCPotential.o LDAFunctional.o \
         PBEFunctional.o BLYPFunctional.o B3LYPFunctional.o  \
+	ExponentialIntegral.o HSEFunctional.o \
 	NonLocalPotential.o SampleReader.o StructuredDocumentHandler.o \
 	SampleHandler.o AtomSetHandler.o WavefunctionHandler.o \
         XMLGFPreprocessor.o Base64Transcoder.o \

src/XMLGFPreprocessor.C

@@ -24,6 +24,7 @@
 #include <cstdlib>
 #include <cstdio>
 #include <cassert>
+#include <unistd.h> // close
 
 #include <vector>
 #include <valarray>

src/notes

 --------------------------------------------------------------------------------
+rel1_57_15
+--------------------------------------------------------------------------------
+r1409: ExchangeOperator.C: fixed coeffs of divergent terms for empty states
+r1407: XMLGFPreprocessor.C: add unistd header for close function
+--------------------------------------------------------------------------------
 rel1_57_14
 --------------------------------------------------------------------------------
 r1396: util/kpgen/kpgen.C: fixed odd numbered case and added test.

src/release.C

@@ -19,5 +19,5 @@
 #include "release.h"
 std::string release(void)
 {
-  return std::string("1.57.14");
+  return std::string("1.57.15");
 }