Merged in fixtotalDOS (pull request #535)

cleanups to total density of states computation

Approved-by: Phani Motamarri

doc/manual/parameters.tex

@@ -2100,7 +2100,7 @@ \subsection{Parameters in section \tt Post-processing Options}
 {\it Default:} false
 
 
-{\it Description:} [Standard] Computes density of states using Lorentzians. Uses specified Temperature for SCF as the broadening parameter. Outputs a file name 'dosData.out' containing two columns with first column indicating the energy in eV and second column indicating the density of states
+{\it Description:} [Standard] Computes density of states using Lorentzians. Uses specified Temperature for SCF as the broadening parameter. Outputs a file name 'dosData.out' containing two columns with first column indicating the energy in eV relative to the Fermi energy and second column indicating the density of states. In case of collinear spin polarization, the second and third columns indicate the spin-up and spin-down density of states.
 
 
 {\it Possible values:} A boolean value (true or false)

include/dft.h

@@ -964,6 +964,7 @@ namespace dftfe
  */
  void
  compute_tdos(const std::vector<std::vector<double>> &eigenValuesInput,
+ const unsigned int highestStateOfInterest,
  const std::string & fileName);
 
  void

src/dft/dft.cc

@@ -1746,7 +1746,9 @@ namespace dftfe
  writeGSElectronDensity("densityQuadData.txt");
 
  if (d_dftParamsPtr->writeDosFile)
- compute_tdos(eigenValues, "dosData.out");
+ compute_tdos(eigenValues,
+ d_dftParamsPtr->highestStateOfInterestForChebFiltering,
+ "dosData.out");
 
  if (d_dftParamsPtr->writeLdosFile)
  compute_ldos(eigenValues, "ldosData.out");

src/dft/dos.cc

@@ -140,34 +140,47 @@ namespace dftfe
 
 
 
- // compute fermi energy
+ // compute tdos
  template <unsigned int FEOrder, unsigned int FEOrderElectro>
  void
  dftClass<FEOrder, FEOrderElectro>::compute_tdos(
  const std::vector<std::vector<double>> &eigenValuesInput,
+ const unsigned int highestStateOfInterest,
  const std::string & dosFileName)
  {
  computing_timer.enter_subsection("DOS computation");
+
+ // from 0th kpoint and 0th spin index
+ int indexFermiEnergy = -1.0;
+ for (int i = 0; i < d_numEigenValues; ++i)
+ if (eigenValuesInput[0][i] >= fermiEnergy)
+ {
+ if (i > indexFermiEnergy)
+ {
+ indexFermiEnergy = i;
+ break;
+ }
+ }
+
+ const unsigned int highestStateOfInterestUsed =
+ (highestStateOfInterest == 0) ? indexFermiEnergy : highestStateOfInterest;
+
+ // from 0th spin as this is only to get a printing range
  std::vector<double> eigenValuesAllkPoints;
  for (int kPoint = 0; kPoint < d_kPointWeights.size(); ++kPoint)
- {
- for (int statesIter = 0; statesIter < eigenValuesInput[0].size();
- ++statesIter)
- {
- eigenValuesAllkPoints.push_back(
- eigenValuesInput[kPoint][statesIter]);
- }
- }
+ for (int statesIter = 0; statesIter < highestStateOfInterestUsed;
+ ++statesIter)
+ eigenValuesAllkPoints.push_back(eigenValuesInput[kPoint][statesIter]);
 
  std::sort(eigenValuesAllkPoints.begin(), eigenValuesAllkPoints.end());
 
- double totalEigenValues = eigenValuesAllkPoints.size();
- double intervalSize = 0.001;
- double sigma = C_kb * d_dftParamsPtr->TVal;
- double lowerBoundEpsilon = 1.5 * eigenValuesAllkPoints[0];
- double upperBoundEpsilon =
+ const double totalEigenValues = eigenValuesAllkPoints.size();
+ const double intervalSize = 0.001;
+ const double sigma = C_kb * d_dftParamsPtr->TVal;
+ const double lowerBoundEpsilon = 1.5 * eigenValuesAllkPoints[0];
+ const double upperBoundEpsilon =
  eigenValuesAllkPoints[totalEigenValues - 1] * 1.5;
- unsigned int numberIntervals =
+ const unsigned int numberIntervals =
  std::ceil((upperBoundEpsilon - lowerBoundEpsilon) / intervalSize);
 
  std::vector<double> densityOfStates, densityOfStatesUp, densityOfStatesDown;
@@ -187,7 +200,7 @@ namespace dftfe
  ++spinType)
  {
  for (unsigned int statesIter = 0;
- statesIter < d_numEigenValues;
+ statesIter < highestStateOfInterestUsed;
  ++statesIter)
  {
  double term1 =
@@ -198,14 +211,30 @@ namespace dftfe
  double denom = term1 * term1 + sigma * sigma;
  if (spinType == 0)
  densityOfStatesUp[epsInt] +=
- (sigma / M_PI) * (1.0 / denom);
+ d_kPointWeights[kPoint] * (sigma / M_PI) *
+ (1.0 / denom) / d_domainVolume;
  else
  densityOfStatesDown[epsInt] +=
- (sigma / M_PI) * (1.0 / denom);
+ d_kPointWeights[kPoint] * (sigma / M_PI) *
+ (1.0 / denom) / d_domainVolume;
  }
  }
  }
  }
+
+ MPI_Allreduce(MPI_IN_PLACE,
+ &densityOfStatesUp[0],
+ densityOfStatesUp.size(),
+ dataTypes::mpi_type_id(&densityOfStatesUp[0]),
+ MPI_SUM,
+ interpoolcomm);
+
+ MPI_Allreduce(MPI_IN_PLACE,
+ &densityOfStatesDown[0],
+ densityOfStatesDown.size(),
+ dataTypes::mpi_type_id(&densityOfStatesDown[0]),
+ MPI_SUM,
+ interpoolcomm);
  }
  else
  {
@@ -215,17 +244,36 @@ namespace dftfe
  double epsValue = lowerBoundEpsilon + epsInt * intervalSize;
  for (int kPoint = 0; kPoint < d_kPointWeights.size(); ++kPoint)
  {
- for (unsigned int statesIter = 0; statesIter < d_numEigenValues;
+ for (unsigned int statesIter = 0;
+ statesIter < highestStateOfInterestUsed;
  ++statesIter)
  {
  double term1 =
  (epsValue - eigenValuesInput[kPoint][statesIter]);
  double denom = term1 * term1 + sigma * sigma;
- densityOfStates[epsInt] +=
- 2.0 * (sigma / M_PI) * (1.0 / denom);
+ densityOfStates[epsInt] += 2.0 * d_kPointWeights[kPoint] *
+ (sigma / M_PI) * (1.0 / denom) /
+ d_domainVolume;
  }
  }
  }
+
+ MPI_Allreduce(MPI_IN_PLACE,
+ &densityOfStates[0],
+ densityOfStates.size(),
+ dataTypes::mpi_type_id(&densityOfStates[0]),
+ MPI_SUM,
+ interpoolcomm);
+ }
+
+
+ if (d_dftParamsPtr->reproducible_output && d_dftParamsPtr->verbosity == 0)
+ {
+ pcout << "Writing tdos File..." << std::endl;
+ if (d_dftParamsPtr->spinPolarized)
+ pcout << "epsValue SpinUpDos SpinDownDos" << std::endl;
+ else
+ pcout << "epsValue Dos" << std::endl;
  }
 
  if (dealii::Utilities::MPI::this_mpi_process(d_mpiCommParent) == 0)
@@ -240,20 +288,60 @@ namespace dftfe
  for (unsigned int epsInt = 0; epsInt < numberIntervals;
  ++epsInt)
  {
- double epsValue = lowerBoundEpsilon + epsInt * intervalSize;
+ double epsValue =
+ lowerBoundEpsilon + epsInt * intervalSize - fermiEnergy;
  outFile << std::setprecision(18) << epsValue * 27.21138602
  << " " << densityOfStatesUp[epsInt] << " "
  << densityOfStatesDown[epsInt] << std::endl;
+ if (d_dftParamsPtr->reproducible_output &&
+ d_dftParamsPtr->verbosity == 0)
+ {
+ double epsValueTrunc =
+ std::floor(1000000000 *
+ (lowerBoundEpsilon +
+ epsInt * intervalSize - fermiEnergy) *
+ 27.21138602) /
+ 1000000000.0;
+ double dosSpinUpTrunc =
+ std::floor(1000000000 * densityOfStatesUp[epsInt]) /
+ 1000000000.0;
+
+ double dosSpinDownTrunc =
+ std::floor(1000000000 * densityOfStatesDown[epsInt]) /
+ 1000000000.0;
+
+
+ pcout << std::fixed << std::setprecision(8)
+ << epsValueTrunc << " " << dosSpinUpTrunc << " "
+ << dosSpinDownTrunc << std::endl;
+ }
  }
  }
  else
  {
  for (unsigned int epsInt = 0; epsInt < numberIntervals;
  ++epsInt)
  {
- double epsValue = lowerBoundEpsilon + epsInt * intervalSize;
+ double epsValue =
+ lowerBoundEpsilon + epsInt * intervalSize - fermiEnergy;
  outFile << std::setprecision(18) << epsValue * 27.21138602
  << " " << densityOfStates[epsInt] << std::endl;
+
+ if (d_dftParamsPtr->reproducible_output &&
+ d_dftParamsPtr->verbosity == 0)
+ {
+ double epsValueTrunc =
+ std::floor(1000000000 *
+ (lowerBoundEpsilon +
+ epsInt * intervalSize - fermiEnergy) *
+ 27.21138602) /
+ 1000000000.0;
+ double dosTrunc =
+ std::floor(1000000000 * densityOfStates[epsInt]) /
+ 1000000000.0;
+ pcout << std::fixed << std::setprecision(8)
+ << epsValueTrunc << " " << dosTrunc << std::endl;
+ }
  }
  }
  }