Modified KohnShamDFTClass to normalize electron density in each scf

src/ksdft/ElectrostaticAllElectronFE.t.cpp

@@ -584,6 +584,7 @@ namespace dftefe
                                                      d_atomCharges,
                                                      d_smearedChargeRadius);
 
+      RealType totNuclearChargeQuad = 0;
       for (size_type iCell = 0; iCell < quadRuleContainer->nCells(); iCell++)
         {
           for (size_type iComp = 0; iComp < d_numComponents; iComp++)
@@ -595,6 +596,7 @@ namespace dftefe
                 {
                   a[quadId] = (RealType)smfunc(j);
                   quadId    = quadId + 1;
+                  totNuclearChargeQuad += a[quadId];
                 }
               RealType *b = a.data();
               d_nuclearChargesDensity
@@ -604,6 +606,22 @@ namespace dftefe
             }
         }
 
+      utils::mpi::MPIAllreduce<utils::MemorySpace::HOST>(
+        utils::mpi::MPIInPlace,
+        &totNuclearChargeQuad,
+        1,
+        utils::mpi::Types<RealType>::getMPIDatatype(),
+        utils::mpi::MPISum,
+        d_feBMTotalCharge->getMPIPatternP2P()->mpiCommunicator());
+
+      double totalAtomCharges =
+        std::accumulate(d_atomCharges.begin(), d_atomCharges.end(), (double)0);
+      // d_nuclearChargesDensity by totalAtomCharges/totNuclearChargeQuad
+      quadrature::scale((RealType)std::abs(totalAtomCharges /
+                                           totNuclearChargeQuad),
+                        *d_nuclearChargesDensity,
+                        *d_linAlgOpContext);
+
       // create the correction quadValuesContainer for analytical solve
 
       const utils::SmearChargePotentialFunction smfuncPot(
@@ -806,6 +824,7 @@ namespace dftefe
             d_atomCharges[iAtom],
             d_smearedChargeRadius[iAtom]);
 
+          RealType nuclearChargeQuad = 0;
           for (size_type iCell = 0; iCell < quadRuleContainer->nCells();
                iCell++)
             {
@@ -818,6 +837,7 @@ namespace dftefe
                     {
                       a[quadId] = (RealType)(*smfunc)(j);
                       quadId    = quadId + 1;
+                      nuclearChargeQuad += a[quadId];
                     }
                   RealType *b = a.data();
                   nuclearChargeDensity
@@ -827,10 +847,18 @@ namespace dftefe
                 }
             }
 
-          // Scale by 4\pi
-          scale((RealType)(4 * utils::mathConstants::pi),
-                nuclearChargeDensity,
-                *d_linAlgOpContext);
+          utils::mpi::MPIAllreduce<utils::MemorySpace::HOST>(
+            utils::mpi::MPIInPlace,
+            &nuclearChargeQuad,
+            1,
+            utils::mpi::Types<RealType>::getMPIDatatype(),
+            utils::mpi::MPISum,
+            d_feBMTotalCharge->getMPIPatternP2P()->mpiCommunicator());
+
+          // Scale by 4\pi * d_atomCharges[iAtom]/nuclearChargeQuad
+          quadrature::scale((RealType)std::abs(4 * utils::mathConstants::pi),
+                            nuclearChargeDensity,
+                            *d_linAlgOpContext);
 
           std::shared_ptr<
             linearAlgebra::LinearSolverFunction<ValueTypeBasisData,

src/ksdft/KohnShamDFT.h

@@ -210,6 +210,11 @@ namespace dftefe
       MixingScheme<RealType, RealType> d_mixingScheme;
       std::shared_ptr<linearAlgebra::LinAlgOpContext<memorySpace>>
         d_linAlgOpContext;
+      linearAlgebra::Vector<ValueTypeWaveFunctionCoeff, memorySpace>
+        d_lanczosGuess;
+      linearAlgebra::MultiVector<ValueTypeWaveFunctionCoeff, memorySpace>
+                      d_waveFunctionSubspaceGuess;
+      const size_type d_numElectrons;
 
     }; // end of KohnShamDFT
   }    // end of namespace ksdft

src/ksdft/KohnShamDFT.t.cpp

@@ -86,6 +86,44 @@ namespace dftefe
           }
         return std::sqrt(normValue);
       }
+
+      template <typename RealType, utils::MemorySpace memorySpace>
+      void
+      normalizeDensityQuadData(
+        quadrature::QuadratureValuesContainer<RealType, memorySpace> &inValues,
+        const size_type                              numElectrons,
+        linearAlgebra::LinAlgOpContext<memorySpace> &linAlgOpContext,
+        const utils::mpi::MPIComm &                  mpiComm)
+      {
+        RealType totalDensityInQuad = 0.0;
+
+        for (size_type iCell = 0; iCell < inValues.nCells(); iCell++)
+          {
+            for (size_type iComp = 0; iComp < inValues.getNumberComponents();
+                 iComp++)
+              {
+                std::vector<RealType> a(inValues.nCellQuadraturePoints(iCell));
+                inValues.template getCellQuadValues<utils::MemorySpace::HOST>(
+                  iCell, iComp, a.data());
+                for (auto j : a)
+                  {
+                    totalDensityInQuad += j;
+                  }
+              }
+          }
+        utils::mpi::MPIAllreduce<utils::MemorySpace::HOST>(
+          utils::mpi::MPIInPlace,
+          &totalDensityInQuad,
+          1,
+          utils::mpi::Types<RealType>::getMPIDatatype(),
+          utils::mpi::MPISum,
+          mpiComm);
+
+        quadrature::scale((RealType)(std::abs((RealType)numElectrons /
+                                              totalDensityInQuad)),
+                          inValues,
+                          linAlgOpContext);
+      }
     } // namespace KohnShamDFTInternal
 
     template <typename ValueTypeElectrostaticsCoeff,
@@ -176,6 +214,16 @@ namespace dftefe
       , d_waveFunctionBatchSize(waveFunctionBatchSize)
       , d_SCFTol(scfDensityResidualNormTolerance)
       , d_rootCout(std::cout)
+      , d_waveFunctionSubspaceGuess(feBMWaveFn->getMPIPatternP2P(),
+                                    linAlgOpContext,
+                                    numWantedEigenvalues,
+                                    0.0,
+                                    1.0)
+      , d_lanczosGuess(feBMWaveFn->getMPIPatternP2P(),
+                       linAlgOpContext,
+                       0.0,
+                       1.0)
+      , d_numElectrons(numElectrons)
     {
       utils::throwException(electronChargeDensityInput.getNumberComponents() ==
                               1,
@@ -189,6 +237,13 @@ namespace dftefe
       int rank;
       utils::mpi::MPICommRank(d_mpiCommDomain, &rank);
       d_rootCout.setCondition(rank == 0);
+
+      // normalize electroncharge density
+      KohnShamDFTInternal::normalizeDensityQuadData(d_densityInQuadValues,
+                                                    numElectrons,
+                                                    *d_linAlgOpContext,
+                                                    d_mpiCommDomain);
+
       d_hamitonianKin = std::make_shared<KineticFE<ValueTypeWaveFunctionBasis,
                                                    ValueTypeWaveFunctionCoeff,
                                                    memorySpace,
@@ -238,22 +293,12 @@ namespace dftefe
 
       // call the eigensolver
 
-      linearAlgebra::MultiVector<ValueTypeWaveFunctionCoeff, memorySpace>
-        waveFunctionSubspaceGuess(feBMWaveFn->getMPIPatternP2P(),
-                                  linAlgOpContext,
-                                  numWantedEigenvalues,
-                                  0.0,
-                                  1.0);
-
-      linearAlgebra::Vector<ValueTypeWaveFunctionCoeff, memorySpace>
-        lanczosGuess(feBMWaveFn->getMPIPatternP2P(), linAlgOpContext, 0.0, 1.0);
-
-      lanczosGuess.updateGhostValues();
-      feBMWaveFn->getConstraints().distributeParentToChild(lanczosGuess, 1);
+      d_lanczosGuess.updateGhostValues();
+      feBMWaveFn->getConstraints().distributeParentToChild(d_lanczosGuess, 1);
 
-      waveFunctionSubspaceGuess.updateGhostValues();
+      d_waveFunctionSubspaceGuess.updateGhostValues();
       feBMWaveFn->getConstraints().distributeParentToChild(
-        waveFunctionSubspaceGuess, numWantedEigenvalues);
+        d_waveFunctionSubspaceGuess, numWantedEigenvalues);
 
       // form the kohn sham operator
       d_ksEigSolve = std::make_shared<
@@ -265,8 +310,8 @@ namespace dftefe
         eigenSolveResidualTolerance,
         chebyshevPolynomialDegree,
         maxChebyshevFilterPass,
-        waveFunctionSubspaceGuess,
-        lanczosGuess,
+        d_waveFunctionSubspaceGuess,
+        d_lanczosGuess,
         waveFunctionBatchSize,
         MContextForInv,
         MInvContext);
@@ -365,6 +410,13 @@ namespace dftefe
           // reinit the components of hamiltonian
           if (scfIter > 0)
             {
+              // normalize electroncharge density each scf
+              KohnShamDFTInternal::normalizeDensityQuadData(
+                d_densityInQuadValues,
+                d_numElectrons,
+                *d_linAlgOpContext,
+                d_mpiCommDomain);
+
               d_hamitonianElec->reinitField(d_densityInQuadValues);
               d_hamitonianXC->reinitField(d_densityInQuadValues);
               std::vector<HamiltonianPtrVariant> hamiltonianComponentsVec{
@@ -431,7 +483,8 @@ namespace dftefe
               d_rootCout << "Ground State Energy: " << totalEnergy << "\n";
             }
 
-          d_rootCout << "Density Residual Norm : " << norm << "\n";
+          if (scfIter > 0)
+            d_rootCout << "Density Residual Norm : " << norm << "\n";
           scfIter += 1;
         }
 

src/ksdft/MixingScheme.t.cpp

@@ -30,6 +30,7 @@ namespace dftefe
     MixingScheme<ValueTypeMixingVariable, ValueTypeWeights>::MixingScheme(
       const utils::mpi::MPIComm &mpiComm)
       : d_mpiComm(mpiComm)
+      , d_anyMixingParameterAdaptive(false)
     {}
 
     template <typename ValueTypeMixingVariable, typename ValueTypeWeights>
@@ -102,7 +103,7 @@ namespace dftefe
       cDensity.resize(c.size());
       std::fill(cDensity.begin(), cDensity.end(), 0.0);
 
-      size_type N             = inHist.size() - 1;
+      int       N             = inHist.size() - 1;
       size_type numQuadPoints = 0;
       if (N > 0)
         numQuadPoints = inHist[0].size();
@@ -117,10 +118,10 @@ namespace dftefe
           for (size_type iQuad = 0; iQuad < numQuadPoints; iQuad++)
             {
               ValueTypeMixingVariable Fn = residualHist[N][iQuad];
-              for (size_type m = 0; m < N; m++)
+              for (int m = 0; m < N; m++)
                 {
                   ValueTypeMixingVariable Fnm = residualHist[N - 1 - m][iQuad];
-                  for (size_type k = 0; k < N; k++)
+                  for (int k = 0; k < N; k++)
                     {
                       ValueTypeMixingVariable Fnk =
                         residualHist[N - 1 - k][iQuad];
@@ -204,16 +205,16 @@ namespace dftefe
     {
       // initialize data structures
       // assumes rho is a mixing variable
-      size_type N = d_variableHistoryIn[mixingVariable::rho].size() - 1;
+      int N = d_variableHistoryIn[mixingVariable::rho].size() - 1;
       if (N > 0)
         {
           size_type NRHS = 1, lda = N, ldb = N;
           std::vector<linearAlgebra::blasLapack::LapackInt> ipiv(N);
           d_A.resize(lda * N);
           d_c.resize(ldb * NRHS);
-          for (size_type i = 0; i < lda * N; i++)
+          for (int i = 0; i < lda * N; i++)
             d_A[i] = 0.0;
-          for (size_type i = 0; i < ldb * NRHS; i++)
+          for (int i = 0; i < ldb * NRHS; i++)
             d_c[i] = 0.0;
 
           for (const auto &key : mixingVariablesList)
@@ -231,7 +232,7 @@ namespace dftefe
             N, NRHS, &d_A[0], lda, &ipiv[0], &d_c[0], ldb, linAlgOpContextHost);
         }
       d_cFinal = 1.0;
-      for (size_type i = 0; i < N; i++)
+      for (int i = 0; i < N; i++)
         d_cFinal -= d_c[i];
       computeAdaptiveAndersonMixingParameter();
     }
@@ -375,7 +376,7 @@ namespace dftefe
           ValueType varInBar =
             d_cFinal * d_variableHistoryIn[mixingVariableName][N][iQuad];
 
-          for (size_type i = 0; i < N; i++)
+          for (int i = 0; i < N; i++)
             {
               varResidualBar +=
                 d_c[i] *

test/ksdft/src/CMakeLists.txt

@@ -18,17 +18,17 @@ add_subdirectory(${MAIN_PROJECT_DIR}/src/ksdft ${MAIN_PROJECT_DIR}/test/ksdft/li
 if(ENABLE_MPI)
 	add_compile_definitions(DFTEFE_WITH_MPI)
 
-	add_executable(TestHamiltonianClassical TestHamiltonianClassical.cpp )
-	target_link_libraries(TestHamiltonianClassical PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms dft-efe-ksdft)
-	set_target_properties(TestHamiltonianClassical PROPERTIES OUTPUT_NAME "TestHamiltonianClassical" SUFFIX ".x")
+	# add_executable(TestHamiltonianClassical TestHamiltonianClassical.cpp )
+	# target_link_libraries(TestHamiltonianClassical PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms dft-efe-ksdft)
+	# set_target_properties(TestHamiltonianClassical PROPERTIES OUTPUT_NAME "TestHamiltonianClassical" SUFFIX ".x")
 
-	add_executable(TestKSAllElectronEigenSolveClassical TestKSAllElectronEigenSolveClassical.cpp )
-	target_link_libraries(TestKSAllElectronEigenSolveClassical PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms dft-efe-ksdft)
-	set_target_properties(TestKSAllElectronEigenSolveClassical PROPERTIES OUTPUT_NAME "TestKSAllElectronEigenSolveClassical" SUFFIX ".x")
+	# add_executable(TestKSAllElectronEigenSolveClassical TestKSAllElectronEigenSolveClassical.cpp )
+	# target_link_libraries(TestKSAllElectronEigenSolveClassical PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms dft-efe-ksdft)
+	# set_target_properties(TestKSAllElectronEigenSolveClassical PROPERTIES OUTPUT_NAME "TestKSAllElectronEigenSolveClassical" SUFFIX ".x")
 
-	add_executable(TestKSAllElectronEigenSolveOrthoEFE TestKSAllElectronEigenSolveOrthoEFE.cpp )
-	target_link_libraries(TestKSAllElectronEigenSolveOrthoEFE PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms dft-efe-ksdft)
-	set_target_properties(TestKSAllElectronEigenSolveOrthoEFE PROPERTIES OUTPUT_NAME "TestKSAllElectronEigenSolveOrthoEFE" SUFFIX ".x")
+	# add_executable(TestKSAllElectronEigenSolveOrthoEFE TestKSAllElectronEigenSolveOrthoEFE.cpp )
+	# target_link_libraries(TestKSAllElectronEigenSolveOrthoEFE PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms dft-efe-ksdft)
+	# set_target_properties(TestKSAllElectronEigenSolveOrthoEFE PROPERTIES OUTPUT_NAME "TestKSAllElectronEigenSolveOrthoEFE" SUFFIX ".x")
 
 	add_executable(TestKSDFTClassical TestKSDFTClassical.cpp )
 	target_link_libraries(TestKSDFTClassical PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms dft-efe-ksdft)

test/ksdft/src/TestKSAllElectronEigenSolveClassical.cpp

@@ -132,7 +132,7 @@ int main()
   // Set up Triangulation
     std::shared_ptr<basis::TriangulationBase> triangulationBase =
         std::make_shared<basis::TriangulationDealiiParallel<dim>>(comm);
-  std::vector<unsigned int>         subdivisions = {30, 30, 30};
+  std::vector<unsigned int>         subdivisions = {10, 10, 10};
   std::vector<bool>                 isPeriodicFlags(dim, false);
   std::vector<utils::Point> domainVectors(dim, utils::Point(dim, 0.0));
 
@@ -274,39 +274,6 @@ int main()
   std::shared_ptr<const quadrature::QuadratureRuleContainer> quadRuleContainer =  
                 feBasisData->getQuadratureRuleContainer();
 
-  std::vector<double> chargeDensity(atomCoordinatesVec.size(), 0.0), mpiReducedChargeDensity(chargeDensity.size(), 0.0);
-
-  const utils::SmearChargeDensityFunction smden(atomCoordinatesVec,
-                                                atomChargesVec,
-                                                smearedChargeRadiusVec);
-
-    double charge = 0;
-    for(size_type i = 0 ; i < quadRuleContainer->nCells() ; i++)
-    {
-      std::vector<double> JxW = quadRuleContainer->getCellJxW(i);
-      size_type quadId = 0;
-      for (auto j : quadRuleContainer->getCellRealPoints(i))
-      {
-        charge += smden(j) * JxW[quadId];
-        quadId = quadId + 1;
-      }
-    }
-    chargeDensity[0] = charge;
-
-  utils::mpi::MPIAllreduce<Host>(
-        chargeDensity.data(),
-        mpiReducedChargeDensity.data(),
-        chargeDensity.size(),
-        utils::mpi::MPIDouble,
-        utils::mpi::MPISum,
-        comm);
-
-  for(size_type i = 0 ; i < atomCoordinatesVec.size() ; i++)
-  {
-    atomChargesVec[i] *= 1/std::abs(mpiReducedChargeDensity[i]);
-  }
-  std::cout << "Total nuclear charge in system: " << mpiReducedChargeDensity[0] << std::endl;
-
    quadrature::QuadratureValuesContainer<double, Host> 
       electronChargeDensity(quadRuleContainer, 1, 0.0);
 
@@ -368,7 +335,7 @@ int main()
                                                   feBasisData,   
                                                   feBasisData,
                                                   feBasisData, 
-                                                  feBasisData                                                                                            
+                                                  feBasisData,                                                                                          
                                                   *externalPotentialFunction,
                                                   linAlgOpContext,
                                                   50);
@@ -555,7 +522,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
   std::cout << "elec energy: " << elecEnergy << "\n";
 
   // calculate band energy
-  double bandEnergy;
+  double bandEnergy = 0;
   for(size_type i = 0 ; i < occupation.size(); i++)
   {
     bandEnergy += 2 * occupation[i] * kohnShamEnergies[i];