Add force_tol and stress_tol variables

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,18 @@ void BOSampleStepper::step(int niter)
     if ( atoms_move )
       s_.constraints.update_constraints(dt);
 
+    // check if maxforce and maxstress within tolerance
+    if ( onpe0 )
+    {
+      if ( force_tol > 0.0 )
+        cout << "  maxforce: " << scientific
+             << setprecision(4) << maxforce << endl;
+      if ( stress_tol > 0.0 )
+        cout << "  maxstress: " << scientific
+             << setprecision(4) << maxstress << endl;
+    }
+    iter_done = ( maxforce <= force_tol ) && ( maxstress <= stress_tol );
+
     // print iteration time
     double time = tm_iter.real();
     double tmin = time;
@@ -1235,6 +1267,7 @@ void BOSampleStepper::step(int niter)
            << endl;
       cout << "</iteration>" << endl;
     }
+
   } // for iter
 
   if ( atoms_move )

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/Makefile

@@ -447,6 +447,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 +732,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 +809,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 +823,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 +871,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/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"
@@ -281,6 +283,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 +299,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));