merged trunk into efield branch

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

src/Bisection.C

@@ -403,11 +403,11 @@ void Bisection::compute_transform(const SlaterDet& sd, int maxsweep, double tol)
 
   // diagonalize projectors
   int nsweep = jade(maxsweep,tol,amat_,*u_,adiag_);
-  if ( nsweep >= maxsweep )
-    if ( gcontext_.onpe0() )
-      cout << "Bisection::compute_transform: nsweep=" << nsweep
-           << " maxsweep=" << maxsweep << " tol=" << tol << endl;
-
+#ifdef TIMING
+  if ( gcontext_.onpe0() )
+    cout << "Bisection::compute_transform: nsweep=" << nsweep
+         << " maxsweep=" << maxsweep << " tol=" << tol << endl;
+#endif
 
 #ifdef TIMING
   tmap["jade"].stop();
@@ -706,13 +706,13 @@ void Bisection::trim_amat(const vector<double>& occ)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-bool Bisection::overlap(int i, int j)
+bool Bisection::overlap(int i, int j) const
 {
   return overlap(localization_,i,j);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-bool Bisection::overlap(vector<long int> loc_, int i, int j)
+bool Bisection::overlap(const vector<long int>& loc_, int i, int j) const
 {
   // overlap: return true if the functions i and j overlap according
   // to the localization vector loc_
@@ -739,25 +739,28 @@ bool Bisection::overlap(vector<long int> loc_, int i, int j)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-double Bisection::pair_fraction(void)
+double Bisection::pair_fraction(void) const
 {
   // pair_fraction: return fraction of pairs having non-zero overlap
+  // count pairs (i,j) having non-zero overlap for i != j only
   int sum = 0;
-  for ( int i = 0; i < nst_; i++ )
+  for ( int i = 1; i < nst_; i++ )
   {
     int count = 0;
-    for ( int j = 0; j < nst_; j++ )
+    for ( int j = i+1; j < nst_; j++ )
     {
       if ( overlap(i,j) )
         count++;
     }
     sum += count;
   }
-  return ((double) sum)/(nst_*nst_);
+  // add overlap with self: (i,i)
+  sum += nst_;
+  return ((double) sum)/((nst_*(nst_+1))/2);
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-double Bisection::size(int i)
+double Bisection::size(int i) const
 {
   // size: return fraction of the domain on which state i is non-zero
   long int loc = localization_[i];
@@ -780,7 +783,7 @@ double Bisection::size(int i)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-double Bisection::total_size(void)
+double Bisection::total_size(void) const
 {
   // total_size: return normalized sum of sizes
   double sum = 0.0;

src/Bisection.h

@@ -79,13 +79,14 @@ class Bisection
 
     int nmat(void) const { return nmat_; }
     long int localization(int i) const { return localization_[i]; }
-    std::vector<long int> localization(void) const { return localization_; }
-    bool overlap(int i, int j);
-    bool overlap(std::vector<long int> loc, int i, int j);
+    const std::vector<long int>& localization(void) const
+    { return localization_; }
+    bool overlap(int i, int j) const;
+    bool overlap(const std::vector<long int>& loc, int i, int j) const;
     const DoubleMatrix& u(void) const { return *u_; }
-    double pair_fraction(void);
-    double size(int i);
-    double total_size(void);
+    double pair_fraction(void) const;
+    double size(int i) const;
+    double total_size(void) const;
     ~Bisection();
 };
 #endif

src/CPSampleStepper.C

@@ -219,8 +219,11 @@ void CPSampleStepper::step(int niter)
         cout << "  <temp_ion> " << mdionic_stepper->temp() << " </temp_ion>\n";
         cout << "  <eta_ion> " << mdionic_stepper->eta() << " </eta_ion>\n";
       }
-      cout << "  <econst> " << enthalpy+ekin_ion+ekin_e << " </econst>\n";
-      cout << "  <ekin_ec> " << energy+ekin_ion+2*ekin_e << " </ekin_ec>\n";
+      double econst = energy + ekin_ion + ekin_e;
+      if ( mdionic_stepper )
+        econst += mdionic_stepper->ekin_stepper();
+      cout << "  <econst> " << econst << " </econst>\n";
+      cout << "  <ekin_ec> " << econst + ekin_e << " </ekin_ec>\n";
     }
 
     if ( compute_stress )

src/jade.C

@@ -198,6 +198,7 @@ int jade(int maxsweep, double tol, vector<DoubleMatrix*> a,
       // skip the pair if one or both of the vectors is a dummy vector
       // i.e. a vector having jglobal==-1
 
+#pragma omp parallel for
       for ( int k = 0; k < a.size(); k++ )
       {
         for ( int ipair = 0; ipair < nploc; ipair++ )
@@ -228,6 +229,7 @@ int jade(int maxsweep, double tol, vector<DoubleMatrix*> a,
       tm_comm.stop();
       // apq now contains the matrix elements
 
+#pragma omp parallel for
       for ( int ipair = 0; ipair < nploc; ipair++ )
       {
         if ( jglobal[top[ipair]] >= 0 &&
@@ -328,6 +330,7 @@ int jade(int maxsweep, double tol, vector<DoubleMatrix*> a,
           drot(&mloc,up,&one,uq,&one,&c,&s);
 
           // new value of off-diag element apq
+          double diag_change_ipair = 0.0;
           for ( int k = 0; k < a.size(); k++ )
           {
             const int iapq = 3*ipair + k*3*nploc;
@@ -337,8 +340,10 @@ int jade(int maxsweep, double tol, vector<DoubleMatrix*> a,
 
             // accumulate change in sum of squares of diag elements
             // note negative sign: decrease in offdiag is increase in diag
-            diag_change -= 2.0 * ( apqnew*apqnew - apq[iapq]*apq[iapq] );
+            diag_change_ipair -= 2.0 * ( apqnew*apqnew - apq[iapq]*apq[iapq] );
           }
+#pragma omp critical
+          diag_change -= diag_change_ipair;
         }
       } // for ipair
 

src/notes

 --------------------------------------------------------------------------------
-<<<<<<< .working
-rel1_56_2_efield1
---------------------------------------------------------------------------------
-rel1_56_1
-ExchangeOperator.C: Fixed permutation of occupation numbers when using bisection
---------------------------------------------------------------------------------
-rel1_56_0
-=======
-rel1_60_5
+Planned merging of oncv branch into trunk. Modifications.
+Species.h: Must check use of namespace to define proj_data object.
+Species.h: ztot_ zcore_ are double. Related to the definition given
+in the species that may not be an integer charge?
+Species.C: use of "or" keyword in test is a C++-11 feature. Revert to || for
+consistency with the rest of the code.
+EnergyFunctional.h: ChargeDensity cd_ not const. Due to core charge update? 
+ChargeDensity.C::update_rhor: use of transform algorithm. Add and subtract
+core charge can be replaced by use of temporary array, simpler loop.
+sample.xsd: remove definition of doubleListType which appears in species.xsd
+--------------------------------------------------------------------------------
+rel1_60_9
+r1656 ExchangeOperator.C: changed default bisection tolerance to 1.0
+--------------------------------------------------------------------------------
+rel1_60_8
+r1655 ExchangeOperator print exch time only if DEBUG defined
+r1652 Bisection.C print nsweep only if TIMING defined
+r1651 ExchangeOperator.C added omp directives in jade.C
+r1649 ExchangeOperator.C, Bisection.C added BISECTION_NSWEEP debug option
+--------------------------------------------------------------------------------
+rel1_60_7
+--------------------------------------------------------------------------------
+r1645: Bisection.C: modified pair_fraction to reduce calls to overlap.
+r1644: optimization of Bisection::overlap speeds up exchange on large problems
+r1644: added detailed timers to ExchangeOperator.C
+--------------------------------------------------------------------------------
+rel1_60_6
 --------------------------------------------------------------------------------
 r1605: added -vlocal option to the plot command. Plot Vloc+VHart+Vxc.
 r1570: cleanup testXMLGFPreprocessor.C

util/qbox_xyz.py

+#!/usr/bin/python
+# qbox_xyz_cell.py: extract atomic positions with cell info
+# from Qbox output
+# use: qbox_xyz_cell.py file.r
+
+import xml.sax
+import sys
+import math
+
+if len(sys.argv) != 2:
+  print "use: ",sys.argv[0]," file.r"
+  sys.exit()
+
+infile = sys.argv[1]
+a0=0.529177
+
+# Qbox output handler to extract and process data
+class QboxOutputHandler(xml.sax.handler.ContentHandler):
+  def __init__(self):
+    self.step = 0
+    self.inAtomset = 0
+    self.inAtom = 0
+    self.inPosition = 0
+
+  def startElement(self, name, attributes):
+    if name == "atomset":
+      self.tau=[]
+      self.atomname=[]
+      self.inAtomset = 1
+    elif (name == "unit_cell") & self.inAtomset:
+      self.cell_a = attributes["a"]
+      self.cell_b = attributes["b"]
+      self.cell_c = attributes["c"]
+    elif (name == "atom") & self.inAtomset:
+      self.atomname.append(attributes["name"])
+      self.inAtom = 1
+    elif (name == "position") & self.inAtom:
+        self.buffer = ""
+        self.inPosition = 1
+
+  def characters(self, data):
+    if self.inPosition:
+      self.buffer += data
+
+  def endElement(self, name):
+    if (name == "atom") and self.inAtomset:
+      self.inAtom = 0
+    if (name == "position") & self.inAtom:
+      pos = self.buffer.split()
+      x = a0*float(pos[0])
+      y = a0*float(pos[1])
+      z = a0*float(pos[2])
+      self.tau.append([x,y,z])
+      self.inPosition = 0
+    elif name == "atomset":
+      self.step += 1
+      print len(self.tau)
+      avec = self.cell_a.split()
+      bvec = self.cell_b.split()
+      cvec = self.cell_c.split()
+      print self.step,\
+      '%.6f'%(a0*float(avec[0])),\
+      '%.6f'%(a0*float(avec[1])),\
+      '%.6f'%(a0*float(avec[2])),\
+      '%.6f'%(a0*float(bvec[0])),\
+      '%.6f'%(a0*float(bvec[1])),\
+      '%.6f'%(a0*float(bvec[2])),\
+      '%.6f'%(a0*float(cvec[0])),\
+      '%.6f'%(a0*float(cvec[1])),\
+      '%.6f'%(a0*float(cvec[2]))
+      for i in range(len(self.tau)):
+        print self.atomname[i],'%.6f'%self.tau[i][0],\
+                               '%.6f'%self.tau[i][1],\
+                               '%.6f'%self.tau[i][2]
+      self.inAtomset = 0
+
+parser = xml.sax.make_parser()
+handler = QboxOutputHandler()
+parser.setContentHandler(handler)
+parser.parse(infile)