Merge branch 'develop'

src/BOSampleStepper.C

@@ -207,6 +207,9 @@ void BOSampleStepper::step(int niter)
   // GS-only calculation:
   const bool gs_only = !atoms_move && !cell_moves;
 
+  const double force_tol = s_.ctrl.force_tol;
+  const double stress_tol = s_.ctrl.stress_tol;
+
   Timer tm_iter;
 
   Preconditioner prec(wf,ef_,s_.ctrl.ecutprec);
@@ -353,7 +356,8 @@ void BOSampleStepper::step(int niter)
     cout << "<net_charge> " << atoms.nel()-wf.nel() << " </net_charge>\n";
 
   // Next line: special case of niter=0: compute GS only
-  for ( int iter = 0; iter < max(niter,1); iter++ )
+  bool iter_done = false;
+  for ( int iter = 0; iter < max(niter,1) && !iter_done; iter++ )
   {
     // ionic iteration
 
@@ -363,6 +367,8 @@ void BOSampleStepper::step(int niter)
       cout << "<iteration count=\"" << iter+1 << "\">\n";
 
     // compute energy and ionic forces using existing wavefunction
+    double maxforce = 0.0;
+    double maxstress = 0.0;
 
     if ( !gs_only )
     {
@@ -380,6 +386,20 @@ void BOSampleStepper::step(int niter)
       tmap["energy"].stop();
       double enthalpy = ef_.enthalpy();
 
+      if ( force_tol > 0.0 )
+      {
+        maxforce = 0.0;
+        for ( int is = 0; is < fion.size(); is++ )
+          for ( int i = 0; i < fion[is].size(); i++ )
+            maxforce = max(maxforce, fabs(fion[is][i]));
+      }
+
+      if ( stress_tol > 0.0 )
+      {
+        for ( int i = 0; i < sigma.size(); i++ )
+          maxstress = max(maxstress, fabs(sigma[i]));
+      }
+
       if ( onpe0 )
       {
         cout << cd_;
@@ -1220,6 +1240,23 @@ void BOSampleStepper::step(int niter)
     if ( atoms_move )
       s_.constraints.update_constraints(dt);
 
+    // if using force_tol or stress_tol, check if maxforce and maxstress
+    // within tolerance
+    if ( force_tol > 0.0 )
+    {
+      if ( onpe0 )
+        cout << "  maxforce: " << scientific
+             << setprecision(4) << maxforce << fixed << endl;
+      iter_done |= ( maxforce < force_tol );
+    }
+    if ( stress_tol > 0.0 )
+    {
+      if ( onpe0 )
+        cout << "  maxstress: " << scientific
+             << setprecision(4) << maxstress << fixed << endl;
+      iter_done |= ( maxstress < stress_tol );
+    }
+
     // print iteration time
     double time = tm_iter.real();
     double tmin = time;
@@ -1235,6 +1272,7 @@ void BOSampleStepper::step(int niter)
            << endl;
       cout << "</iteration>" << endl;
     }
+
   } // for iter
 
   if ( atoms_move )

src/CGIonicStepper.C

@@ -25,7 +25,7 @@ CGIonicStepper::CGIonicStepper(Sample& s) : IonicStepper(s),
   cgopt_(CGOptimizer(3*natoms_))
 {
   cgopt_.set_alpha_start(1.0);
-  cgopt_.set_alpha_max(5.0);
+  cgopt_.set_alpha_max(50.0);
   cgopt_.set_beta_max(10.0);
 #ifdef DEBUG
   if ( s.ctxt_.onpe0() )
@@ -63,7 +63,7 @@ void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
 
   // check largest displacement
   // max_disp: largest acceptable displacement
-  const double max_disp = 0.05;
+  const double max_disp = 0.2;
   double largest_disp = 0.0;
   for ( int i = 0; i < xp.size(); i++ )
     largest_disp = max(largest_disp,fabs(xp[i]-x[i]));
@@ -74,8 +74,6 @@ void CGIonicStepper::compute_r(double e0, const vector<vector<double> >& f0)
     // rescale displacement and reset the CG optimizer
     double fac = max_disp/largest_disp;
     xp = x + fac * (xp - x);
-
-    // reduce alpha starting value in CG optmizer
     cgopt_.set_alpha_start(fac*cgopt_.alpha_start());
     cgopt_.reset();
   }

src/Control.h

@@ -72,6 +72,8 @@ struct Control
   double btHF;
 
   double scf_tol;
+  double force_tol;
+  double stress_tol;
 
   D3vector e_field;
   std::string polarization;

src/LineMinimizer.C

@@ -189,7 +189,7 @@ double LineMinimizer::next_alpha(double alpha, double f, double fp)
   {
     // we are already in bracketing mode
     nstep_++;
-    if ( nstep_ > nstep_max_ )
+    if ( nstep_max_ > 0 && nstep_ > nstep_max_ )
     {
       if ( debug_print )
         cout << "LineMinimizer: fail, nstep_max" << endl;

src/LineMinimizer.h

@@ -53,6 +53,7 @@ class LineMinimizer
   void set_sigma2(double s) { sigma2_ = s; }
   void set_alpha_start(double a) { alpha_start_ = a; }
   void set_alpha_max(double a) { alpha_max_ = a; }
+  void set_nstep_max(int n) { nstep_max_ = n; }
   void set_debug_print(void) { debug_print = true; }
 
   double next_alpha(double alpha, double f, double fp);

src/Makefile

@@ -58,6 +58,9 @@ OBJECTS=qb.o AtomSet.o Atom.o Species.o \
         Function3d.o Function3dHandler.o \
         $(PLTOBJECTS)
 CXXFLAGS += -DTARGET='"$(TARGET)"'
+# to include VERSION info in release string, use:
+# env VERSION=`git describe` make
+CXXFLAGS += -DVERSION='"$(VERSION)"'
  $(EXEC):    $(OBJECTS)
 	$(LD) $(DFLAGS) -o $@ $^ $(LDFLAGS)
  lib:   $(OBJECTS)
@@ -447,6 +450,9 @@ FermiTemp.o: Wavefunction.h Control.h
 FoldInWsCmd.o: UserInterface.h Sample.h AtomSet.h Context.h blacs.h Atom.h
 FoldInWsCmd.o: D3vector.h UnitCell.h D3tensor.h blas.h ConstraintSet.h
 FoldInWsCmd.o: ExtForceSet.h Wavefunction.h Control.h
+ForceTol.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
+ForceTol.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
+ForceTol.o: Control.h
 FourierTransform.o: FourierTransform.h Timer.h Basis.h D3vector.h UnitCell.h
 FourierTransform.o: blas.h
 FourierTransform.o: Timer.h
@@ -729,6 +735,9 @@ StrainCmd.o: ExtForceSet.h Wavefunction.h Control.h
 Stress.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
 Stress.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
 Stress.o: Control.h
+StressTol.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h
+StressTol.o: UnitCell.h D3tensor.h blas.h ConstraintSet.h ExtForceSet.h
+StressTol.o: Wavefunction.h Control.h
 StructureFactor.o: StructureFactor.h Basis.h D3vector.h UnitCell.h
 StructuredDocumentHandler.o: StructuredDocumentHandler.h StrX.h
 StructuredDocumentHandler.o: StructureHandler.h
@@ -803,7 +812,6 @@ Xc.o: Control.h
 isodate.o: isodate.h
 jacobi.o: blacs.h Context.h Matrix.h blas.h
 jade.o: blacs.h Context.h Matrix.h blas.h Timer.h
-kpgen.o: D3vector.h
 qb.o: isodate.h release.h qbox_xmlns.h uuid_str.h Context.h blacs.h
 qb.o: UserInterface.h Sample.h AtomSet.h Atom.h D3vector.h UnitCell.h
 qb.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
@@ -818,10 +826,11 @@ qb.o: FourierTransform.h StrainCmd.h TorsionCmd.h BisectionCmd.h Bisection.h
 qb.o: SlaterDet.h Basis.h Matrix.h AlphaPBE0.h AlphaRSH.h AtomsDyn.h
 qb.o: BetaRSH.h BlHF.h BtHF.h Cell.h CellDyn.h CellLock.h CellMass.h
 qb.o: ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h Debug.h Dspin.h Ecut.h
-qb.o: Ecutprec.h Ecuts.h Efield.h Polarization.h Emass.h ExtStress.h
-qb.o: FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h MuRSH.h Nempty.h NetCharge.h
-qb.o: Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h Thermostat.h ThTemp.h
-qb.o: ThTime.h ThWidth.h Vext.h ExternalPotential.h WfDiag.h WfDyn.h Xc.h
+qb.o: Ecutprec.h Ecuts.h Efield.h ForceTol.h Polarization.h Emass.h
+qb.o: ExtStress.h FermiTemp.h IterCmd.h IterCmdPeriod.h Dt.h MuRSH.h Nempty.h
+qb.o: NetCharge.h Nrowmax.h Nspin.h RefCell.h ScfTol.h Stress.h StressTol.h
+qb.o: Thermostat.h ThTemp.h ThTime.h ThWidth.h Vext.h ExternalPotential.h
+qb.o: WfDiag.h WfDyn.h Xc.h
 qbox_xmlns.o: qbox_xmlns.h
 release.o: release.h
 sinft.o: sinft.h
@@ -865,6 +874,7 @@ testWavefunction.o: SlaterDet.h Basis.h Matrix.h Timer.h
 testXCFunctional.o: LDAFunctional.h XCFunctional.h PBEFunctional.h Timer.h
 testXMLGFPreprocessor.o: Context.h blacs.h Matrix.h XMLGFPreprocessor.h
 test_fftw.o: Timer.h readTSC.h
+test_sym.o: Basis.h D3vector.h UnitCell.h
 test_vext.o: Function3d.h D3vector.h
 testjacobi.o: Timer.h Context.h blacs.h Matrix.h jacobi.h
 testjade.o: Timer.h Context.h blacs.h Matrix.h jade.h

src/PlotCmd.C

@@ -146,7 +146,7 @@ int PlotCmd::action(int argc, char **argv)
     }
     else if ( !strcmp(argv[iarg],"-spin") )
     {
-      if ( !(plot_density || plot_wf || plot_wfs) )
+      if ( !(plot_density || plot_wf || plot_wfs || plot_vlocal) )
       {
         if ( ui->onpe0() )
           cout << usage << endl;

src/qb.C

@@ -94,6 +94,7 @@ using namespace std;
 #include "Ecutprec.h"
 #include "Ecuts.h"
 #include "Efield.h"
+#include "ForceTol.h"
 #include "Polarization.h"
 #include "Emass.h"
 #include "ExtStress.h"
@@ -109,6 +110,7 @@ using namespace std;
 #include "RefCell.h"
 #include "ScfTol.h"
 #include "Stress.h"
+#include "StressTol.h"
 #include "Thermostat.h"
 #include "ThTemp.h"
 #include "ThTime.h"
@@ -155,7 +157,11 @@ int main(int argc, char **argv, char **envp)
   cout << "                   I http://qboxcode.org      I\n";
   cout << "                   ============================\n\n";
   cout << "\n";
-  cout << "<release> " << release() << " " << TARGET << " </release>" << endl;
+  cout << "<release> " << release() << " " << TARGET;
+#ifdef VERSION
+  cout << " " << VERSION;
+#endif
+  cout << " </release>" << endl;
 
   // Identify executable name, checksum, size and link date
   if ( getlogin() != 0 )
@@ -281,6 +287,7 @@ int main(int argc, char **argv, char **envp)
   ui.addVar(new Ecutprec(s));
   ui.addVar(new Ecuts(s));
   ui.addVar(new Efield(s));
+  ui.addVar(new ForceTol(s));
   ui.addVar(new Polarization(s));
   ui.addVar(new Emass(s));
   ui.addVar(new ExtStress(s));
@@ -296,6 +303,7 @@ int main(int argc, char **argv, char **envp)
   ui.addVar(new RefCell(s));
   ui.addVar(new ScfTol(s));
   ui.addVar(new Stress(s));
+  ui.addVar(new StressTol(s));
   ui.addVar(new Thermostat(s));
   ui.addVar(new ThTemp(s));
   ui.addVar(new ThTime(s));

src/testCGOptimizer.C

@@ -54,6 +54,7 @@ int main(int argc, char** argv)
 
   CGOptimizer cgop(n);
   cgop.set_alpha_start(0.01);
+  cgop.set_beta_max(20);
   cgop.set_debug_print();
 
   valarray<double> x(n), xp(n), g(n);