Merged in energyResidual (pull request #579)

Energy residual criteria

Approved-by: Sambit Das
Approved-by: Phani Motamarri

doc/manual/parameters.tex

@@ -2285,6 +2285,23 @@ \subsection{Parameters in section \tt SCF parameters}
 {\it Description:} [Advanced] Boolean parameter specifying whether to compute the total energy at the end of every SCF. Setting it to false can lead to some computational time savings. Default value is false but is internally set to true if VERBOSITY==5
 
 
+{\it Possible values:} A boolean value (true or false)
+\item {\it Parameter name:} {\tt USE ENERGY RESIDUAL METRIC}
+\phantomsection\label{parameters:SCF parameters/USE ENERGY RESIDUAL METRIC}
+\label{parameters:SCF_20parameters/USE_20ENERGY_20RESIDUAL_20METRIC}
+
+
+\index[prmindex]{USE ENERGY RESIDUAL METRIC}
+\index[prmindexfull]{SCF parameters!USE ENERGY RESIDUAL METRIC}
+{\it Value:} false
+
+
+{\it Default:} false
+
+
+{\it Description:} [Advanced] Boolean parameter specifying whether to use the energy residual metric (equation 7.23 of Richard Matrin second edition) for convergence check. Setting it to false can lead to some computational time savings. Default value is false
+
+
 {\it Possible values:} A boolean value (true or false)
 \item {\it Parameter name:} {\tt CONSTRAINT MAGNETIZATION}
 \phantomsection\label{parameters:SCF parameters/CONSTRAINT MAGNETIZATION}

include/dft.h

@@ -1437,7 +1437,9 @@ namespace dftfe
  // std::map<dealii::CellId, std::vector<double>> d_phiInValues,
  // d_phiOutValues;
  dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
- d_phiInQuadValues, d_phiOutQuadValues;
+ d_phiInQuadValues, d_phiOutQuadValues;
+ dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
+ d_gradPhiInQuadValues, d_gradPhiOutQuadValues, d_gradPhiResQuadValues;
  MixingScheme d_mixingScheme;
 
  distributedCPUVec<double> d_rhoInNodalValuesRead, d_rhoOutNodalValuesSplit,

include/dftParameters.h

@@ -50,7 +50,8 @@ namespace dftfe
  double radiusAtomBall, mixingParameter, spinMixingEnhancementFactor;
  bool adaptAndersonMixingParameter;
  double absLinearSolverTolerance, selfConsistentSolverTolerance, TVal,
- start_magnetization, absLinearSolverToleranceHelmholtz;
+ selfConsistentSolverEnergyTolerance, start_magnetization,
+ absLinearSolverToleranceHelmholtz;
 
  bool isPseudopotential, periodicX, periodicY, periodicZ, useSymm,
  timeReversal, pseudoTestsFlag, constraintMagnetization, writeDosFile,
@@ -114,6 +115,7 @@ namespace dftfe
  unsigned int subspaceRotDofsBlockSize;
  unsigned int nbandGrps;
  bool computeEnergyEverySCF;
+ bool useEnergyResidualTolerance;
  unsigned int scalapackParalProcs;
  unsigned int scalapackBlockSize;
  unsigned int natoms;

include/energyCalculator.h

@@ -147,6 +147,50 @@ namespace dftfe
  const bool print,
  const bool smearedNuclearCharges = false);
 
+ double
+ computeEnergyResidual(
+ const std::shared_ptr<
+ dftfe::basis::FEBasisOperations<dataTypes::number,
+ double,
+ dftfe::utils::MemorySpace::HOST>>
+ &basisOperationsPtr,
+ const std::shared_ptr<
+ dftfe::basis::
+ FEBasisOperations<double, double, dftfe::utils::MemorySpace::HOST>>
+ & basisOperationsPtrElectro,
+ const unsigned int densityQuadratureID,
+ const unsigned int densityQuadratureIDElectro,
+ const unsigned int smearedChargeQuadratureIDElectro,
+ const unsigned int lpspQuadratureIDElectro,
+ const std::shared_ptr<excManager> excManagerPtr,
+ const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
+ &phiTotRhoInValues,
+ const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
+ & phiTotRhoOutValues,
+ const distributedCPUVec<double> &phiTotRhoIn,
+ const distributedCPUVec<double> &phiTotRhoOut,
+ const std::vector<
+ dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
+ &densityInValues,
+ const std::vector<
+ dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
+ &densityOutValues,
+ const std::vector<
+ dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
+ &gradDensityInValues,
+ const std::vector<
+ dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
+ & gradDensityOutValues,
+ const std::map<dealii::CellId, std::vector<double>> &rhoCoreValues,
+ const std::map<dealii::CellId, std::vector<double>> &gradRhoCoreValues,
+ const std::map<dealii::CellId, std::vector<double>> &smearedbValues,
+ const std::map<dealii::CellId, std::vector<unsigned int>>
+ & smearedbNonTrivialAtomIds,
+ const std::vector<std::vector<double>> &localVselfs,
+ const std::map<dealii::types::global_dof_index, double>
+ & atomElectrostaticNodeIdToChargeMap,
+ const bool smearedNuclearCharges);
+
 
  void
  computeXCEnergyTermsSpinPolarized(

src/dft/dft.cc

@@ -2300,6 +2300,7 @@ namespace dftfe
  //
  unsigned int scfIter = 0;
  double norm = 1.0;
+ double energyResidual = 1.0;
  d_rankCurrentLRD = 0;
  d_relativeErrorJacInvApproxPrevScfLRD = 100.0;
  // CAUTION: Choosing a looser tolerance might lead to failed tests
@@ -2309,8 +2310,7 @@ namespace dftfe
  pcout << std::endl;
  if (d_dftParamsPtr->verbosity == 0)
  pcout << "Starting SCF iterations...." << std::endl;
- while ((norm > d_dftParamsPtr->selfConsistentSolverTolerance) &&
- (scfIter < d_dftParamsPtr->numSCFIterations))
+ while (!scfConverged && (scfIter < d_dftParamsPtr->numSCFIterations))
  {
  dealii::Timer local_timer(d_mpiCommParent, true);
  if (d_dftParamsPtr->verbosity >= 1)
@@ -2541,7 +2541,10 @@ namespace dftfe
  d_phiTotRhoIn = d_phiTotRhoOut;
  computing_timer.leave_subsection("density mixing");
 
- if (!(norm > d_dftParamsPtr->selfConsistentSolverTolerance))
+ if (!((norm > d_dftParamsPtr->selfConsistentSolverTolerance) ||
+ (d_dftParamsPtr->useEnergyResidualTolerance &&
+ energyResidual >
+ d_dftParamsPtr->selfConsistentSolverEnergyTolerance)))
  scfConverged = true;
 
  if (d_dftParamsPtr->multipoleBoundaryConditions)
@@ -3303,8 +3306,9 @@ namespace dftfe
  //
  // phiTot with rhoOut
  //
- if (d_dftParamsPtr->computeEnergyEverySCF &&
- d_numEigenValuesRR == d_numEigenValues)
+ if ((d_dftParamsPtr->computeEnergyEverySCF &&
+ d_numEigenValuesRR == d_numEigenValues) ||
+ d_dftParamsPtr->useEnergyResidualTolerance)
  {
  if (d_dftParamsPtr->verbosity >= 2)
  pcout
@@ -3406,25 +3410,44 @@ namespace dftfe
  d_phiTotRhoOut,
  d_phiOutQuadValues,
  dummy);
-
-
- //
- // impose integral phi equals 0
- //
- /*
- if(d_dftParamsPtr->periodicX && d_dftParamsPtr->periodicY &&
- d_dftParamsPtr->periodicZ && !d_dftParamsPtr->pinnedNodeForPBC)
- {
- if(d_dftParamsPtr->verbosity>=2)
- pcout<<"Value of integPhiOut:
- "<<totalCharge(d_dofHandlerPRefined,d_phiTotRhoOut);
- }
- */
-
  computing_timer.leave_subsection("phiTot solve");
-
- const dealii::Quadrature<3> &quadrature =
- matrix_free_data.get_quadrature(d_densityQuadratureId);
+ }
+ if (d_dftParamsPtr->useEnergyResidualTolerance)
+ {
+ computing_timer.enter_subsection("Energy residual computation");
+ energyResidual = energyCalc.computeEnergyResidual(
+ d_basisOperationsPtrHost,
+ d_basisOperationsPtrElectroHost,
+ d_densityQuadratureId,
+ d_densityQuadratureIdElectro,
+ d_smearedChargeQuadratureIdElectro,
+ d_lpspQuadratureIdElectro,
+ d_excManagerPtr,
+ d_phiInQuadValues,
+ d_phiOutQuadValues,
+ d_phiTotRhoIn,
+ d_phiTotRhoOut,
+ d_densityInQuadValues,
+ d_densityOutQuadValues,
+ d_gradDensityInQuadValues,
+ d_gradDensityOutQuadValues,
+ d_rhoCore,
+ d_gradRhoCore,
+ d_bQuadValuesAllAtoms,
+ d_bCellNonTrivialAtomIds,
+ d_localVselfs,
+ d_atomNodeIdToChargeMap,
+ d_dftParamsPtr->smearedNuclearCharges);
+ if (d_dftParamsPtr->verbosity >= 1)
+ pcout << "Energy residual : " << energyResidual << std::endl;
+ if (d_dftParamsPtr->reproducible_output)
+ pcout << "Energy residual : " << std::setprecision(4)
+ << energyResidual << std::endl;
+ computing_timer.leave_subsection("Energy residual computation");
+ }
+ if (d_dftParamsPtr->computeEnergyEverySCF &&
+ d_numEigenValuesRR == d_numEigenValues)
+ {
  d_dispersionCorr.computeDispresionCorrection(
  atomLocations, d_domainBoundingVectors);
  const double totalEnergy = energyCalc.computeEnergy(