Made changes to easily shift from p-DoFHandler to 2p-DofHandler and v…

…ice-versa for Helmholtz solve. Also made changes to get the relaxation work in the case of Kerker mixing. Further added option to do only pinning or pinning plus meanValueConstraint for Fully periodic case

include/constants.h

@@ -13,26 +13,29 @@
 //
 // ---------------------------------------------------------------------
 //
-// @author Sambit Das
+// @author Sambit Das and Phani Motamarri
 //
 
 #ifndef constants_H_
 #define constants_H_
 
 namespace dftfe {
  //
- //Add prefic C_ to all constants
+ //Add prefix C_ to all constants
  //
 
  /// Boltzmann constant
- const double C_kb =3.166811429e-06;
+ const double C_kb = 3.166811429e-06;
 
  /// problem space dimensions
- const int C_DIM=3;
+ const int C_DIM = 3;
 
  /// 1d quadrature rule order
  template <unsigned int FEOrder> constexpr unsigned int C_num1DQuad(){return FEOrder+1;}
 
+ //kerker Helmholtz solve polynomial Order
+ template <unsigned int FEOrder> constexpr unsigned int C_num1DKerkerPoly(){return FEOrder;}
+
 
  /// 1d quadrature rule order for non-local part of pseudopotential
  template <unsigned int FEOrder> constexpr unsigned int C_num1DQuadPSP()

include/dft.h

@@ -470,9 +470,9 @@ namespace dftfe {
  double mixing_anderson_spinPolarized();
  double mixing_broyden();
  double mixing_broyden_spinPolarized();
- double nodalDensity_mixing_simple(kerkerSolverProblem<2*FEOrder> & solverProblem,
+ double nodalDensity_mixing_simple(kerkerSolverProblem<C_num1DKerkerPoly<FEOrder>()> & solverProblem,
  dealiiLinearSolver & dealiiLinearSolver);
- double nodalDensity_mixing_anderson(kerkerSolverProblem<2*FEOrder> & solverProblem,
+ double nodalDensity_mixing_anderson(kerkerSolverProblem<C_num1DKerkerPoly<FEOrder>()> & solverProblem,
  dealiiLinearSolver & dealiiLinearSolver);
 
 

include/dftParameters.h

@@ -40,7 +40,7 @@ namespace dftfe {
  extern double radiusAtomBall, mixingParameter;
  extern double lowerEndWantedSpectrum,absLinearSolverTolerance,selfConsistentSolverTolerance,TVal, start_magnetization,absLinearSolverToleranceHelmholtz;
 
- extern bool isPseudopotential, periodicX, periodicY, periodicZ, useSymm, timeReversal,pseudoTestsFlag, constraintMagnetization, writeDosFile, writeLdosFile,writeLocalizationLengths;
+ extern bool isPseudopotential, periodicX, periodicY, periodicZ, useSymm, timeReversal,pseudoTestsFlag, constraintMagnetization, writeDosFile, writeLdosFile,writeLocalizationLengths, pinnedNodeForPBC;
  extern std::string meshFileName,coordinatesFile,domainBoundingVectorsFile,kPointDataFile, ionRelaxFlagsFile, orthogType, algoType, pseudoPotentialFile;
 
  extern double outerAtomBallRadius, innerAtomBallRadius, meshSizeOuterDomain;

src/dft/charge.cc

@@ -121,7 +121,7 @@ template <unsigned int FEOrder>
 double dftClass<FEOrder>::totalCharge(const dealii::MatrixFree<3,double> & matrixFreeDataObject,
  const vectorType & nodalField)
 {
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<2*FEOrder>(),1,double> fe_evalField(matrixFreeDataObject);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<C_num1DKerkerPoly<FEOrder>()>(),1,double> fe_evalField(matrixFreeDataObject);
  VectorizedArray<double> normValueVectorized = make_vectorized_array(0.0);
  const unsigned int numQuadPoints = fe_evalField.n_q_points;
  for(unsigned int cell = 0; cell < matrixFreeDataObject.n_macro_cells(); ++cell)
@@ -148,8 +148,9 @@ double dftClass<FEOrder>::totalCharge(const dealii::MatrixFree<3,double> & matri
 
 }
 
-
+//
 //compute total charge
+//
 template <unsigned int FEOrder>
 double dftClass<FEOrder>::totalMagnetization(const std::map<dealii::CellId, std::vector<double> > *rhoQuadValues){
  double normValue=0.0;
@@ -172,12 +173,15 @@ double dftClass<FEOrder>::totalMagnetization(const std::map<dealii::CellId, std:
  return Utilities::MPI::sum(normValue, mpi_communicator);
 }
 
+//
+//compute field l2 norm
+//
 template <unsigned int FEOrder>
 double dftClass<FEOrder>::fieldl2Norm(const dealii::MatrixFree<3,double> & matrixFreeDataObject,
  const vectorType & nodalField)
 
 {
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<2*FEOrder>(),1,double> fe_evalField(matrixFreeDataObject);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<C_num1DKerkerPoly<FEOrder>()>(),1,double> fe_evalField(matrixFreeDataObject);
  VectorizedArray<double> normValueVectorized = make_vectorized_array(0.0);
  const unsigned int numQuadPoints = fe_evalField.n_q_points;
  for(unsigned int cell = 0; cell < matrixFreeDataObject.n_macro_cells(); ++cell)

src/dft/density.cc

@@ -33,7 +33,7 @@ void dftClass<FEOrder>::compute_rhoOut(const bool isConsiderSpectrumSplitting)
  // 
 
  //fill in rhoOutValues and gradRhoOutValues
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<FEOrder>(),1,double> rhoEval(d_matrixFreeDataPRefined,0,1);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<FEOrder>(),1,double> rhoEval(d_matrixFreeDataPRefined,0,1);
  const unsigned int numQuadPoints = rhoEval.n_q_points;
  DoFHandler<C_DIM>::active_cell_iterator subCellPtr;
  for(unsigned int cell = 0; cell < d_matrixFreeDataPRefined.n_macro_cells(); ++cell)
@@ -244,6 +244,8 @@ template<unsigned int FEOrder>
 void dftClass<FEOrder>::noRemeshRhoDataInit()
 {
 
+ if(dftParameters::mixingMethod == "ANDERSON_WITH_KERKER")
+ d_rhoInNodalValues = d_rhoOutNodalValues;
 
 
  //create temporary copies of rho Out data
@@ -306,7 +308,7 @@ void dftClass<FEOrder>::computeRhoNodalFromPSI(bool isConsiderSpectrumSplitting)
  const unsigned int dofs_per_cell = d_dofHandlerPRefined.get_fe().dofs_per_cell;
  typename DoFHandler<3>::active_cell_iterator cell = d_dofHandlerPRefined.begin_active(), endc = d_dofHandlerPRefined.end();
  dealii::IndexSet locallyOwnedDofs = d_dofHandlerPRefined.locally_owned_dofs();
- QGaussLobatto<3> quadrature_formula(2*FEOrder+1);
+ QGaussLobatto<3> quadrature_formula(C_num1DKerkerPoly<FEOrder>()+1);
  const unsigned int numQuadPoints = quadrature_formula.size();
 
  AssertThrow(numQuadPoints == matrix_free_data.get_n_q_points(3),ExcMessage("Number of quadrature points from Quadrature object does not match with number of quadrature points obtained from matrix_free object"));
@@ -398,10 +400,10 @@ void dftClass<FEOrder>::computeRhoFromPSI(std::map<dealii::CellId, std::vector<d
 
 #ifdef USE_COMPLEX
  FEEvaluation<3,FEOrder,C_num1DQuad<FEOrder>(),2> psiEval(matrix_free_data,eigenDofHandlerIndex,0);
- FEEvaluation<3,FEOrder,2*FEOrder+1,2> psiEvalRefined(matrix_free_data,eigenDofHandlerIndex,3);
+ FEEvaluation<3,FEOrder,C_num1DKerkerPoly<FEOrder>()+1,2> psiEvalRefined(matrix_free_data,eigenDofHandlerIndex,3);
 #else
  FEEvaluation<3,FEOrder,C_num1DQuad<FEOrder>(),1> psiEval(matrix_free_data,eigenDofHandlerIndex,0);
- FEEvaluation<3,FEOrder,2*FEOrder+1,1> psiEvalRefined(matrix_free_data,eigenDofHandlerIndex,3);
+ FEEvaluation<3,FEOrder,C_num1DKerkerPoly<FEOrder>()+1,1> psiEvalRefined(matrix_free_data,eigenDofHandlerIndex,3);
 #endif
 
  const unsigned int numQuadPoints= lobattoNodesFlag?psiEvalRefined.n_q_points:psiEval.n_q_points;

src/dft/dft.cc

@@ -677,7 +677,9 @@ namespace dftfe {
  //
  //rho init (use previous ground state electron density)
  //
+ if(dftParameters::mixingMethod != "ANDERSON_WITH_KERKER")
  solveNoSCF();
+
  noRemeshRhoDataInit();
  }
 
@@ -902,7 +904,7 @@ namespace dftfe {
  //set up solver functions for Helmholtz to be used only when Kerker mixing is on
  //use 2p dofHandler
  //
- kerkerSolverProblem<2*FEOrder> kerkerPreconditionedResidualSolverProblem(mpi_communicator);
+ kerkerSolverProblem<C_num1DKerkerPoly<FEOrder>()> kerkerPreconditionedResidualSolverProblem(mpi_communicator);
  if(dftParameters::mixingMethod=="ANDERSON_WITH_KERKER")
  kerkerPreconditionedResidualSolverProblem.init(d_matrixFreeDataPRefined,
  d_constraintsPRefined,
@@ -932,7 +934,6 @@ namespace dftfe {
  dftParameters::lowerEndWantedSpectrum,
  0.0);
 
-
  //
  //precompute shapeFunctions and shapeFunctionGradients and shapeFunctionGradientIntegrals
  //
@@ -1082,7 +1083,7 @@ namespace dftfe {
  //
  //impose integral phi equals 0
  //
- if(dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ)
+ if(dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ && !dftParameters::pinnedNodeForPBC)
  {
  double integPhi = totalCharge(dofHandler,
  d_phiTotRhoIn);
@@ -1479,7 +1480,7 @@ namespace dftfe {
  //
  //impose integral phi equals 0
  //
- if(dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ)
+ if(dftParameters::periodicX && dftParameters::periodicY && dftParameters::periodicZ && !dftParameters::pinnedNodeForPBC)
  {
  double integPhi = totalCharge(dofHandler,
  d_phiTotRhoOut);

src/dft/initBoundaryConditions.cc

@@ -193,7 +193,7 @@ void dftClass<FEOrder>::initBoundaryConditions(){
  quadratureVector.push_back(QGauss<1>(C_num1DQuad<FEOrder>()));
  quadratureVector.push_back(QGaussLobatto<1>(FEOrder+1));
  quadratureVector.push_back(QGauss<1>(C_num1DQuadPSP<FEOrder>()));
- quadratureVector.push_back(QGaussLobatto<1>(2*FEOrder+1));
+ quadratureVector.push_back(QGaussLobatto<1>(C_num1DKerkerPoly<FEOrder>()+1));
 
  //
  //

src/dft/initRho.cc

@@ -217,7 +217,7 @@ void dftClass<FEOrder>::initRho()
  gradRhoInValues=&(gradRhoInVals.back());
  }
 
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<FEOrder>(),1,double> rhoEval(d_matrixFreeDataPRefined,0,1);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<FEOrder>(),1,double> rhoEval(d_matrixFreeDataPRefined,0,1);
  const unsigned int numQuadPoints = rhoEval.n_q_points; 
  DoFHandler<C_DIM>::active_cell_iterator subCellPtr;
  for(unsigned int cell = 0; cell < d_matrixFreeDataPRefined.n_macro_cells(); ++cell)

src/dft/nodalDensityMixingSchemes.cc

@@ -18,7 +18,7 @@
 
 
 template<unsigned int FEOrder>
-double dftClass<FEOrder>::nodalDensity_mixing_simple(kerkerSolverProblem<2*FEOrder> & kerkerPreconditionedResidualSolverProblem,
+double dftClass<FEOrder>::nodalDensity_mixing_simple(kerkerSolverProblem<C_num1DKerkerPoly<FEOrder>()> & kerkerPreconditionedResidualSolverProblem,
  dealiiLinearSolver & dealiiCGSolver)
 {
  double normValue=0.0;
@@ -42,7 +42,7 @@ double dftClass<FEOrder>::nodalDensity_mixing_simple(kerkerSolverProblem<2*FEOrd
 
 
  //create FEEval object to be used subsequently
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<2*FEOrder>(),1,double> fe_evalHelm(d_matrixFreeDataPRefined);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<C_num1DKerkerPoly<FEOrder>()>(),1,double> fe_evalHelm(d_matrixFreeDataPRefined);
  unsigned int numQuadPoints = fe_evalHelm.n_q_points; 
  DoFHandler<C_DIM>::active_cell_iterator subCellPtr;
 
@@ -100,7 +100,7 @@ double dftClass<FEOrder>::nodalDensity_mixing_simple(kerkerSolverProblem<2*FEOrd
  d_rhoInNodalVals.push_back(d_rhoInNodalValues);
 
 
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<FEOrder>(),1,double> fe_evalRho(d_matrixFreeDataPRefined,0,1);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<FEOrder>(),1,double> fe_evalRho(d_matrixFreeDataPRefined,0,1);
  numQuadPoints = fe_evalRho.n_q_points;
  //compute rho and grad rho for computing Veff using the rhoIn computed above
  for(unsigned int cell = 0; cell < d_matrixFreeDataPRefined.n_macro_cells(); ++cell)
@@ -144,7 +144,7 @@ double dftClass<FEOrder>::nodalDensity_mixing_simple(kerkerSolverProblem<2*FEOrd
 
 //implement nodal anderson mixing scheme with Kerker
 template<unsigned int FEOrder>
-double dftClass<FEOrder>::nodalDensity_mixing_anderson(kerkerSolverProblem<2*FEOrder> & kerkerPreconditionedResidualSolverProblem,
+double dftClass<FEOrder>::nodalDensity_mixing_anderson(kerkerSolverProblem<C_num1DKerkerPoly<FEOrder>()> & kerkerPreconditionedResidualSolverProblem,
  dealiiLinearSolver & dealiiCGSolver)
 {
  double normValue=0.0;
@@ -281,7 +281,7 @@ double dftClass<FEOrder>::nodalDensity_mixing_anderson(kerkerSolverProblem<2*FEO
  diffRhoBar.update_ghost_values();
 
  //create FEEval object to be used subsequently
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<2*FEOrder>(),1,double> fe_evalHelm(d_matrixFreeDataPRefined);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<C_num1DKerkerPoly<FEOrder>()>(),1,double> fe_evalHelm(d_matrixFreeDataPRefined);
  unsigned int numQuadPoints = fe_evalHelm.n_q_points; 
  DoFHandler<C_DIM>::active_cell_iterator subCellPtr;
 
@@ -342,7 +342,7 @@ double dftClass<FEOrder>::nodalDensity_mixing_anderson(kerkerSolverProblem<2*FEO
  d_rhoInNodalVals.push_back(d_rhoInNodalValues);
 
 
- FEEvaluation<C_DIM,2*FEOrder,C_num1DQuad<FEOrder>(),1,double> fe_evalRho(d_matrixFreeDataPRefined,0,1);
+ FEEvaluation<C_DIM,C_num1DKerkerPoly<FEOrder>(),C_num1DQuad<FEOrder>(),1,double> fe_evalRho(d_matrixFreeDataPRefined,0,1);
  numQuadPoints = fe_evalRho.n_q_points;
  for(unsigned int cell = 0; cell < d_matrixFreeDataPRefined.n_macro_cells(); ++cell)
  {

src/dft/pRefinedDoFHandler.cc

@@ -24,7 +24,7 @@
 template <unsigned int FEOrder>
 void dftClass<FEOrder>::createpRefinedDofHandler(dealii::parallel::distributed::Triangulation<3> & triaObject)
 {
- d_dofHandlerPRefined.initialize(triaObject,dealii::FE_Q<3>(dealii::QGaussLobatto<1>(2*FEOrder+1)));
+ d_dofHandlerPRefined.initialize(triaObject,dealii::FE_Q<3>(dealii::QGaussLobatto<1>(C_num1DKerkerPoly<FEOrder>()+1)));
  d_dofHandlerPRefined.distribute_dofs(d_dofHandlerPRefined.get_fe());
 
  dealii::IndexSet locallyRelevantDofs;
@@ -100,7 +100,7 @@ void dftClass<FEOrder>::initpRefinedObjects()
  matrixFreeDofHandlerVectorInput.push_back(&d_dofHandlerPRefined);
 
  std::vector<Quadrature<1> > quadratureVector;
- quadratureVector.push_back(QGauss<1>(C_num1DQuad<2*FEOrder>()));
+ quadratureVector.push_back(QGauss<1>(C_num1DQuad<C_num1DKerkerPoly<FEOrder>()>()));
  quadratureVector.push_back(QGauss<1>(C_num1DQuad<FEOrder>()));
 
  d_matrixFreeDataPRefined.reinit(matrixFreeDofHandlerVectorInput,

src/helmholtz/kerkerSolverProblem.cc

@@ -220,7 +220,7 @@ namespace dftfe {
  }
 
 
-
+ template class kerkerSolverProblem<1>;
  template class kerkerSolverProblem<2>;
  template class kerkerSolverProblem<3>;
  template class kerkerSolverProblem<4>;

utils/dftParameters.cc

@@ -59,6 +59,7 @@ namespace dftfe
  bool electrostaticsHRefinement = false;
  bool electrostaticsPRefinement = false;
  bool meshAdaption = false;
+ bool pinnedNodeForPBC = true;
 
  std::string startingWFCType="";
  bool useBatchGEMM=false;
@@ -265,6 +266,10 @@ namespace dftfe
  Patterns::Bool(),
  "[Standard] Periodicity along the third domain bounding vector.");
 
+ prm.declare_entry("POINT WISE DIRICHLET CONSTRAINT ", "true",
+ Patterns::Bool(),
+ "[Developer] Flag to set point wise dirichlet constraints to eliminate null-space associated with the discretized Poisson operator subject to periodic BCs.")
+
  prm.declare_entry("CONSTRAINTS PARALLEL CHECK", "true",
  Patterns::Bool(),
  "[Developer] Check for consistency of constraints in parallel.");
@@ -672,6 +677,7 @@ namespace dftfe
  dftParameters::periodicZ = prm.get_bool("PERIODIC3");
  dftParameters::constraintsParallelCheck = prm.get_bool("CONSTRAINTS PARALLEL CHECK");
  dftParameters::createConstraintsFromSerialDofhandler= prm.get_bool("CONSTRAINTS FROM SERIAL DOFHANDLER");
+ dftParameters::pinnedNodeForPBC = prm.get_bool("POINT WISE DIRICHLET CONSTRAINT");
  }
  prm.leave_subsection ();