Modified mesh generator to create consistent refinement across period…

…ic surfaces.

include/dft.h

@@ -137,7 +137,8 @@ namespace dftfe {
  /**
  * @brief Does KSDFT problem pre-processing steps but without remeshing.
  */
- void initNoRemesh(bool flag = true);
+ void initNoRemesh(const bool updateImageKPoints = true,
+ const bool useSingleAtomSolution = false );
 
  /**
  * @brief Selects between only electronic field relaxation or combined electronic and geometry relxation

include/triangulationManager.h

@@ -197,13 +197,13 @@ namespace dftfe {
  parallel::distributed::Triangulation<3> & parallelTriangulationB);
 
  /**
- * @brief generates reset meshes to the last mesh generated by auto mesh approach. This 
+ * @brief generates reset meshes to the last mesh generated by auto mesh approach. This
  * is used in Gaussian update of meshes during electrostatics
  */
  void generateResetMeshes(const std::vector<std::vector<double> > & domainBoundingVectors,
  const bool generateSerialTria,
  const bool generateElectrostaticsTria);
- 
+
 
  /**
  * @brief serialize the triangulations and the associated solution vectors
@@ -277,17 +277,8 @@ namespace dftfe {
  parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho,
  parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp,
  parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce,
- const bool generateElectrostaticsTria);
-
- /**
- * @brief internal function which generates a parallel mesh using a adaptive refinement strategy.
- *
- */
- void generateMesh(parallel::distributed::Triangulation<3>& parallelTriangulation,
- parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho,
- parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp,
- parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce,
- const bool generateElectrostaticsTria);
+ const bool generateElectrostaticsTria,
+ const bool generateSerialTria=false);
 
  /**
  * @brief internal function which generates a coarse mesh which is required for the load function call in
@@ -299,16 +290,38 @@ namespace dftfe {
  /**
  * @brief internal function which sets refinement flags based on a custom created algorithm
  *
+ * @return bool boolean flag is any local cell has refinement flag set
+ */
+ bool refinementAlgorithmA(parallel::distributed::Triangulation<3>& parallelTriangulation,
+ parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho,
+ parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp,
+ parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce,
+ const bool generateElectrostaticsTria,
+ std::vector<unsigned int> & locallyOwnedCellsRefineFlags,
+ std::map<dealii::CellId,unsigned int> & cellIdToCellRefineFlagMapLocal,
+ const bool smoothenCellsOnPeriodicBoundary=false);
+
+ /**
+ * @brief internal function which sets refinement flags to have consistent refinement across periodic
+ * boundary
+ *
+ * @return bool boolean flag is any local cell has refinement flag set
  */
- void refinementAlgorithmA(parallel::distributed::Triangulation<3>& parallelTriangulation,
+ bool consistentPeriodicBoundaryRefinement(parallel::distributed::Triangulation<3>& parallelTriangulation,
  parallel::distributed::Triangulation<3>& electrostaticsTriangulationRho,
  parallel::distributed::Triangulation<3>& electrostaticsTriangulationDisp,
  parallel::distributed::Triangulation<3>& electrostaticsTriangulationForce,
  const bool generateElectrostaticsTria,
  std::vector<unsigned int> & locallyOwnedCellsRefineFlags,
  std::map<dealii::CellId,unsigned int> & cellIdToCellRefineFlagMapLocal);
 
-
+ /**
+ * @brief check that triangulation has consistent refinement across periodic boundary including
+ * for ghost cells
+ *
+ */
+ bool checkPeriodicSurfaceRefinementConsistency(parallel::distributed::Triangulation<3>& parallelTriangulation);
+
 
  /**
  * @brief internal function which refines the serial mesh based on refinement flags from parallel mesh.
@@ -351,7 +364,7 @@ namespace dftfe {
  std::vector<std::vector<double> > d_atomPositions;
  std::vector<std::vector<double> > d_imageAtomPositions;
  std::vector<std::vector<double> > d_domainBoundingVectors;
- const unsigned int d_max_refinement_steps=20;
+ const unsigned int d_max_refinement_steps=40;
 
  //
  //parallel objects

src/dft/dft.cc

@@ -338,7 +338,7 @@ namespace dftfe {
  "to half the number of electrons with a 15 percent buffer to avoid convergence issues in"
  "SCF iterations"<<std::endl;
  }
- d_numEigenValues = (numElectrons/2.0) + std::max(0.15*(numElectrons/2.0),20.0);
+ d_numEigenValues = (numElectrons/2.0) + std::max(0.2*(numElectrons/2.0),20.0);
 
  if(dftParameters::verbosity >= 1)
  {
@@ -645,10 +645,11 @@ namespace dftfe {
  }
 
  template<unsigned int FEOrder>
- void dftClass<FEOrder>::initNoRemesh(bool flag)
+ void dftClass<FEOrder>::initNoRemesh(const bool updateImageKPoints,
+ const bool useSingleAtomSolution)
  {
  computingTimerStandard.enter_section("KSDFT problem initialization");
- if(flag)
+ if(updateImageKPoints)
  initImageChargesUpdateKPoints();
 
  if (dftParameters::isIonOpt)
@@ -658,11 +659,19 @@ namespace dftfe {
  //
  initBoundaryConditions();
 
- //
- //rho init (use previous ground state electron density)
- //
- solveNoSCF();
- noRemeshRhoDataInit();
+ if (useSingleAtomSolution)
+ {
+ readPSI();
+ initRho();
+ }
+ else
+ {
+ //
+ //rho init (use previous ground state electron density)
+ //
+ solveNoSCF();
+ noRemeshRhoDataInit();
+ }
 
  //
  //reinitialize pseudopotential related data structures
@@ -1583,7 +1592,7 @@ namespace dftfe {
  true);
 
  d_groundStateEnergy = totalEnergy;
- 
+
  }
 
  MPI_Barrier(interpoolcomm);

src/dft/moveAtoms.cc

@@ -100,7 +100,7 @@ namespace internal{
 template<unsigned int FEOrder>
 void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<1,3,double> > & globalAtomsDisplacements,
  double maximumForceToBeRelaxed)
- 
+
 {
  const int numberGlobalAtoms = atomLocations.size();
  int numberImageCharges = d_imageIds.size();
@@ -164,7 +164,7 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  0,
  MPI_COMM_WORLD);
 
- 
+
 
 
 
@@ -220,7 +220,7 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  atomLocationsFractional[iAtom][2]=newFracCoord[0];
  atomLocationsFractional[iAtom][3]=newFracCoord[1];
  atomLocationsFractional[iAtom][4]=newFracCoord[2];
- 
+
  }
 
  initImageChargesUpdateKPoints(false);
@@ -262,12 +262,12 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  controlPointLocationsInitialMove.push_back(d_closestTriaVertexToAtomsLocation[iAtom]);
  controlPointDisplacementsInitialMove.push_back(d_dispClosestTriaVerticesToAtoms[iAtom]);
  controlPointDisplacementsCurrentMove.push_back(d_gaussianMovementAtomsNetDisplacements[iAtom]);
- 
+
  }
- 
+
  MPI_Barrier(mpi_communicator);
  d_autoMesh=0;
- 
+
  const bool useHybridMeshUpdateScheme = true;// dftParameters::electrostaticsHRefinement?false:true;
 
  if(!useHybridMeshUpdateScheme)//always remesh
@@ -302,7 +302,7 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  /*d_mesh.resetMesh(d_mesh.getParallelMeshUnmoved(),
  d_mesh.getParallelMeshMoved());*/
 
- 
+
 
  d_mesh.generateResetMeshes(d_domainBoundingVectors,
  dftParameters::useSymm
@@ -349,10 +349,10 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  // MPI_INT,
  // 0,
  // MPI_COMM_WORLD);
- 
+
  if(useGaussian!=1)
  {
- if (!dftParameters::reproducible_output) 
+ if (!dftParameters::reproducible_output)
  pcout << "Auto remeshing and reinitialization of dft problem for new atom coordinates as max net displacement magnitude: "<<maxDispAtom<< " is greater than: "<< break1 << " Bohr..." << std::endl;
  init(0);
 
@@ -368,7 +368,7 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  }
  else
  {
- if (!dftParameters::reproducible_output) 
+ if (!dftParameters::reproducible_output)
  pcout << "Trying to Move using Gaussian with same Gaussian constant for computing the forces: "<<forcePtr->getGaussianGeneratorParameter()<<" as net max displacement magnitude: "<< maxDispAtom<< " is below " << break1 <<" Bohr"<<std::endl;
  if (!dftParameters::reproducible_output)
  pcout << "Max current disp magnitude: "<<maxCurrentDispAtom<<" Bohr"<<std::endl;
@@ -381,7 +381,7 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  factor = 1.25;
  else if(maximumForceToBeRelaxed < 1e-04)
  factor = 1.15;
- 
+
  if (meshQualityMetrics.first || meshQualityMetrics.second > factor*d_autoMeshMaxJacobianRatio)
  d_autoMesh=1;
  MPI_Bcast(&(d_autoMesh),
@@ -437,12 +437,12 @@ void dftClass<FEOrder>::updateAtomPositionsAndMoveMesh(const std::vector<Tensor<
  }
  else
  {
- if (!dftParameters::reproducible_output) 
+ if (!dftParameters::reproducible_output)
  pcout<< " Mesh quality check for Gaussian movement of mesh along with atoms: maximum jacobian ratio after movement: "<< meshQualityMetrics.second<<std::endl;
  if (!dftParameters::reproducible_output)
  pcout << "Now Reinitializing all moved triangulation dependent objects..." << std::endl;
 
- initNoRemesh(false);
+ initNoRemesh(false,maxCurrentDispAtom>0.06?true:false);
  if (!dftParameters::reproducible_output)
  pcout << "...Reinitialization end" << std::endl;
  }

src/solvers/eigenSolvers/chebyshevOrthogonalizedSubspaceIterationSolver.cc

@@ -33,21 +33,21 @@ namespace dftfe{
  else if(upperBoundUnwantedSpectrum > 500 && upperBoundUnwantedSpectrum <= 1000)
  chebyshevOrder = 30;
  else if(upperBoundUnwantedSpectrum > 1000 && upperBoundUnwantedSpectrum <= 1500)
- chebyshevOrder = 34;
+ chebyshevOrder = 39;
  else if(upperBoundUnwantedSpectrum > 1500 && upperBoundUnwantedSpectrum <= 2000)
- chebyshevOrder = 38;
+ chebyshevOrder = 43;
  else if(upperBoundUnwantedSpectrum > 2000 && upperBoundUnwantedSpectrum <= 3000)
- chebyshevOrder = 45;
+ chebyshevOrder = 50;
  else if(upperBoundUnwantedSpectrum > 3000 && upperBoundUnwantedSpectrum <= 4000)
- chebyshevOrder = 53;
+ chebyshevOrder = 58;
  else if(upperBoundUnwantedSpectrum > 4000 && upperBoundUnwantedSpectrum <= 5000)
- chebyshevOrder = 60;
+ chebyshevOrder = 65;
  else if(upperBoundUnwantedSpectrum > 5000 && upperBoundUnwantedSpectrum <= 9000)
- chebyshevOrder = 68;
+ chebyshevOrder = 73;
  else if(upperBoundUnwantedSpectrum > 9000 && upperBoundUnwantedSpectrum <= 14000)
- chebyshevOrder = 94;
+ chebyshevOrder = 100;
  else if(upperBoundUnwantedSpectrum > 14000 && upperBoundUnwantedSpectrum <= 20000)
- chebyshevOrder = 109;
+ chebyshevOrder = 115;
  else if(upperBoundUnwantedSpectrum > 20000 && upperBoundUnwantedSpectrum <= 30000)
  chebyshevOrder = 162;
  else if(upperBoundUnwantedSpectrum > 30000 && upperBoundUnwantedSpectrum <= 50000)