Skip to content

Commit

Permalink
Removed bugs in KohnShamDFT class
Browse files Browse the repository at this point in the history
  • Loading branch information
Avirup Sircar committed Jun 5, 2024
1 parent 6521062 commit 8a4968b
Show file tree
Hide file tree
Showing 8 changed files with 182 additions and 58 deletions.
1 change: 1 addition & 0 deletions src/ksdft/DensityCalculator.t.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -355,6 +355,7 @@ namespace dftefe
& waveFunc,
quadrature::QuadratureValuesContainer<RealType, memorySpace> &rho)
{
rho.setValue((RealType)0);
const size_type numLocallyOwnedCells = d_feBMPsi->nLocallyOwnedCells();

utils::MemoryTransfer<memorySpace, memorySpace> memoryTransfer;
Expand Down
17 changes: 11 additions & 6 deletions src/ksdft/ElectrostaticAllElectronFE.t.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -590,13 +590,14 @@ namespace dftefe
for (size_type iComp = 0; iComp < d_numComponents; iComp++)
{
size_type quadId = 0;
std::vector<double> jxw = quadRuleContainer->getCellJxW(iCell);
std::vector<RealType> a(
quadRuleContainer->nCellQuadraturePoints(iCell));
for (auto j : quadRuleContainer->getCellRealPoints(iCell))
{
a[quadId] = (RealType)smfunc(j);
quadId = quadId + 1;
totNuclearChargeQuad += a[quadId];
totNuclearChargeQuad += (RealType)smfunc(j) * jxw[quadId];
quadId = quadId + 1;
}
RealType *b = a.data();
d_nuclearChargesDensity
Expand Down Expand Up @@ -830,14 +831,16 @@ namespace dftefe
{
for (size_type iComp = 0; iComp < d_numComponents; iComp++)
{
size_type quadId = 0;
size_type quadId = 0;
std::vector<double> jxw =
quadRuleContainer->getCellJxW(iCell);
std::vector<RealType> a(
quadRuleContainer->nCellQuadraturePoints(iCell));
for (auto j : quadRuleContainer->getCellRealPoints(iCell))
{
a[quadId] = (RealType)(*smfunc)(j);
quadId = quadId + 1;
nuclearChargeQuad += a[quadId];
nuclearChargeQuad += (RealType)(*smfunc)(j)*jxw[quadId];
quadId = quadId + 1;
}
RealType *b = a.data();
nuclearChargeDensity
Expand All @@ -856,7 +859,9 @@ namespace dftefe
d_feBMTotalCharge->getMPIPatternP2P()->mpiCommunicator());

// Scale by 4\pi * d_atomCharges[iAtom]/nuclearChargeQuad
quadrature::scale((RealType)std::abs(4 * utils::mathConstants::pi),
quadrature::scale((RealType)std::abs(4 * utils::mathConstants::pi *
d_atomCharges[iAtom] /
nuclearChargeQuad),
nuclearChargeDensity,
*d_linAlgOpContext);

Expand Down
120 changes: 79 additions & 41 deletions src/ksdft/KohnShamDFT.t.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,41 +88,54 @@ namespace dftefe
}

template <typename RealType, utils::MemorySpace memorySpace>
void
RealType
normalizeDensityQuadData(
quadrature::QuadratureValuesContainer<RealType, memorySpace> &inValues,
const size_type numElectrons,
const size_type numElectrons,
const utils::MemoryStorage<RealType, utils::MemorySpace::HOST> &JxW,
linearAlgebra::LinAlgOpContext<memorySpace> &linAlgOpContext,
const utils::mpi::MPIComm & mpiComm)
const utils::mpi::MPIComm & mpiComm,
bool computeTotalDensity,
bool scaleDensity)
{
RealType totalDensityInQuad = 0.0;

for (size_type iCell = 0; iCell < inValues.nCells(); iCell++)
if (computeTotalDensity || scaleDensity)
{
for (size_type iComp = 0; iComp < inValues.getNumberComponents();
iComp++)
int quadId = 0;
for (size_type iCell = 0; iCell < inValues.nCells(); iCell++)
{
std::vector<RealType> a(inValues.nCellQuadraturePoints(iCell));
inValues.template getCellQuadValues<utils::MemorySpace::HOST>(
iCell, iComp, a.data());
for (auto j : a)
for (size_type iComp = 0;
iComp < inValues.getNumberComponents();
iComp++)
{
totalDensityInQuad += j;
std::vector<RealType> a(
inValues.nCellQuadraturePoints(iCell));
inValues
.template getCellQuadValues<utils::MemorySpace::HOST>(
iCell, iComp, a.data());
for (auto j : a)
{
totalDensityInQuad += j * *(JxW.data() + quadId);
quadId += 1;
}
}
}
utils::mpi::MPIAllreduce<utils::MemorySpace::HOST>(
utils::mpi::MPIInPlace,
&totalDensityInQuad,
1,
utils::mpi::Types<RealType>::getMPIDatatype(),
utils::mpi::MPISum,
mpiComm);
}
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);

if (scaleDensity)
quadrature::scale((RealType)(std::abs((RealType)numElectrons /
totalDensityInQuad)),
inValues,
linAlgOpContext);

return totalDensityInQuad;
}
} // namespace KohnShamDFTInternal

Expand Down Expand Up @@ -238,11 +251,23 @@ namespace dftefe
utils::mpi::MPICommRank(d_mpiCommDomain, &rank);
d_rootCout.setCondition(rank == 0);

//************* CHANGE THIS **********************
utils::MemoryTransfer<utils::MemorySpace::HOST, memorySpace> memTransfer;
auto jxwData = feBDHamiltonian->getJxWInAllCells();
d_jxwDataHost.resize(jxwData.size());
memTransfer.copy(jxwData.size(), d_jxwDataHost.data(), jxwData.data());

// normalize electroncharge density
KohnShamDFTInternal::normalizeDensityQuadData(d_densityInQuadValues,
numElectrons,
*d_linAlgOpContext,
d_mpiCommDomain);
RealType totalDensityInQuad =
KohnShamDFTInternal::normalizeDensityQuadData(d_densityInQuadValues,
numElectrons,
d_jxwDataHost,
*d_linAlgOpContext,
d_mpiCommDomain,
true,
true);

d_rootCout << "Electron density in : " << totalDensityInQuad << "\n";

d_hamitonianKin = std::make_shared<KineticFE<ValueTypeWaveFunctionBasis,
ValueTypeWaveFunctionCoeff,
Expand Down Expand Up @@ -325,12 +350,6 @@ namespace dftefe
linAlgOpContext,
cellBlockSize,
waveFunctionBatchSize);

//************* CHANGE THIS **********************
utils::MemoryTransfer<utils::MemorySpace::HOST, memorySpace> memTransfer;
auto jxwData = feBDHamiltonian->getJxWInAllCells();
d_jxwDataHost.resize(jxwData.size());
memTransfer.copy(jxwData.size(), d_jxwDataHost.data(), jxwData.data());
}

template <typename ValueTypeElectrostaticsCoeff,
Expand Down Expand Up @@ -411,11 +430,18 @@ namespace dftefe
if (scfIter > 0)
{
// normalize electroncharge density each scf
KohnShamDFTInternal::normalizeDensityQuadData(
d_densityInQuadValues,
d_numElectrons,
*d_linAlgOpContext,
d_mpiCommDomain);
RealType totalDensityInQuad =
KohnShamDFTInternal::normalizeDensityQuadData(
d_densityInQuadValues,
d_numElectrons,
d_jxwDataHost,
*d_linAlgOpContext,
d_mpiCommDomain,
true,
false);

d_rootCout << "Electron density in : " << totalDensityInQuad
<< "\n";

d_hamitonianElec->reinitField(d_densityInQuadValues);
d_hamitonianXC->reinitField(d_densityInQuadValues);
Expand Down Expand Up @@ -450,6 +476,18 @@ namespace dftefe
d_kohnShamWaveFunctions,
d_densityOutQuadValues);

RealType totalDensityInQuad =
KohnShamDFTInternal::normalizeDensityQuadData(
d_densityOutQuadValues,
d_numElectrons,
d_jxwDataHost,
*d_linAlgOpContext,
d_mpiCommDomain,
true,
false);

d_rootCout << "Electron density out : " << totalDensityInQuad << "\n";

// check residual in density if else
if (d_evaluateEnergyEverySCF)
{
Expand All @@ -470,7 +508,7 @@ namespace dftefe
d_rootCout << "LDA EXC energy: " << xcEnergy << "\n";

// calculate band energy
RealType bandEnergy;
RealType bandEnergy = 0;
for (size_type i = 0; i < d_occupation.size(); i++)
{
bandEnergy += 2 * d_occupation[i] * d_kohnShamEnergies[i];
Expand Down Expand Up @@ -507,7 +545,7 @@ namespace dftefe
d_rootCout << "LDA EXC energy: " << xcEnergy << "\n";

// calculate band energy
RealType bandEnergy;
RealType bandEnergy = 0;
for (size_type i = 0; i < d_occupation.size(); i++)
{
bandEnergy += 2 * d_occupation[i] * d_kohnShamEnergies[i];
Expand Down
3 changes: 3 additions & 0 deletions src/quadrature/QuadratureValuesContainer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,9 @@ namespace dftefe
const size_type componentId,
const ValueType *values);

void
setValue(const ValueType value);

template <utils::MemorySpace memorySpaceDst>
void
getCellValues(const size_type cellId, ValueType *values) const;
Expand Down
8 changes: 8 additions & 0 deletions src/quadrature/QuadratureValuesContainer.t.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,14 @@ namespace dftefe
d_storage.template copyFrom<memorySpaceSrc>(values, size, 0, offset);
}

template <typename ValueType, utils::MemorySpace memorySpace>
void
QuadratureValuesContainer<ValueType, memorySpace>::setValue(
const ValueType value)
{
d_storage.setValue(value);
}

template <typename ValueType, utils::MemorySpace memorySpace>
template <utils::MemorySpace memorySpaceSrc>
void
Expand Down
37 changes: 35 additions & 2 deletions test/ksdft/src/TestKSAllElectronEigenSolveClassical.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {10, 10, 10};
std::vector<unsigned int> subdivisions = {20, 20, 20};
std::vector<bool> isPeriodicFlags(dim, false);
std::vector<utils::Point> domainVectors(dim, utils::Point(dim, 0.0));

Expand Down Expand Up @@ -274,6 +274,39 @@ int main()
std::shared_ptr<const quadrature::QuadratureRuleContainer> quadRuleContainer =
feBasisData->getQuadratureRuleContainer();

std::vector<double> atomChargesVecScaled(atomCoordinatesVec.size());
for(size_type i = 0 ; i < atomCoordinatesVec.size() ; i++)
{
std::vector<double> chargeDensity(1, 0.0), mpiReducedChargeDensity(1, 0.0);

const utils::SmearChargeDensityFunction smden(atomCoordinatesVec[i],
atomChargesVec[i],
smearedChargeRadiusVec[i]);

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);

atomChargesVecScaled[i] *= 1/std::abs(mpiReducedChargeDensity[0]);
}

quadrature::QuadratureValuesContainer<double, Host>
electronChargeDensity(quadRuleContainer, 1, 0.0);

Expand All @@ -290,7 +323,7 @@ int main()
std::shared_ptr<const utils::ScalarSpatialFunctionReal> smfunc =
std::make_shared<const utils::SmearChargePotentialFunction>(
atomCoordinatesVec,
atomChargesVec,
atomChargesVecScaled,
smearedChargeRadiusVec);

std::shared_ptr<const basis::FEBasisManager
Expand Down
38 changes: 37 additions & 1 deletion test/ksdft/src/TestKSAllElectronEigenSolveOrthoEFE.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -438,6 +438,42 @@ int main()
// std::shared_ptr<basis::FEBasisDataStorage<double, Host>> efeBasisDataAdaptiveHamiltonian = efeBasisDataAdaptive;

std::shared_ptr<const quadrature::QuadratureRuleContainer> quadRuleContainer =
efeBasisDataAdaptive->getQuadratureRuleContainer();

std::vector<double> atomChargesVecScaled(atomCoordinatesVec.size());
for(size_type i = 0 ; i < atomCoordinatesVec.size() ; i++)
{
std::vector<double> chargeDensity(1, 0.0), mpiReducedChargeDensity(1, 0.0);

const utils::SmearChargeDensityFunction smden(atomCoordinatesVec[i],
atomChargesVec[i],
smearedChargeRadiusVec[i]);

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);

atomChargesVecScaled[i] *= 1/std::abs(mpiReducedChargeDensity[0]);
}

quadRuleContainer =
efeBasisDataOrbitalAdaptive->getQuadratureRuleContainer();

quadrature::QuadratureValuesContainer<double, Host>
Expand All @@ -456,7 +492,7 @@ int main()
std::shared_ptr<const utils::ScalarSpatialFunctionReal> smfunc =
std::make_shared<const utils::SmearChargePotentialFunction>(
atomCoordinatesVec,
atomChargesVec,
atomChargesVecScaled,
smearedChargeRadiusVec);

std::shared_ptr<const basis::FEBasisManager
Expand Down
Loading

0 comments on commit 8a4968b

Please sign in to comment.