Skip to content

Commit

Permalink
Changes done - Added reinit in PoinLinearSolverFunction , added local…
Browse files Browse the repository at this point in the history 
…ly_relevant_dofs as parameter to ConstraintsLocalDealii , Optimized ConstraintsInternal functions , optimized adaptive quadrature, removed boost::legendre from SphericalHarmonicFunctions.cpp
  • Loading branch information
Avirup Sircar committed Sep 20, 2024
1 parent 6bf7925 commit 20f7f6e
Show file tree
Hide file tree
Showing 52 changed files with 1,091 additions and 267 deletions.
54 changes: 40 additions & 14 deletions analysis/classicalEnrichmentComparison/Forces/KSDFTClassical/TestKohnShamDft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,12 +120,19 @@ T readParameter(std::string ParamFile, std::string param, utils::ConditionalOStr
std::vector<double>
operator()(const std::vector<utils::Point> &points) const
{
double ylm00 = atoms::Clm(0, 0) * atoms::Dm(0) * atoms::Plm(0, 0, 1) * atoms::Qm(0, 0);
std::vector<double> ret(0);
ret.resize(points.size());
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = (*this)(points[i]);
}
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
auto vec = d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "density");
for(auto &enrichmentObjId : vec)
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = ret[i] + std::abs(enrichmentObjId->getValue(points[i], origin) * (1/ylm00));
}
}
return ret;
}
};
Expand Down Expand Up @@ -176,10 +183,19 @@ T readParameter(std::string ParamFile, std::string param, utils::ConditionalOStr
{
std::vector<double> ret(0);
ret.resize(points.size());
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = (*this)(points[i]);
}
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
auto vec = d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "vnuclear");
for (unsigned int i = 0 ; i < points.size() ; i++)
{
for(auto &enrichmentObjId : vec)
{
ret[i] = ret[i] + std::abs(enrichmentObjId->getValue(points[i], origin) *
((*d_b)(points[i])));
}
}
}
return ret;
}
};
Expand Down Expand Up @@ -229,10 +245,19 @@ T readParameter(std::string ParamFile, std::string param, utils::ConditionalOStr
{
std::vector<double> ret(0);
ret.resize(points.size());
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = (*this)(points[i]);
}
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
auto vec = d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "orbital");
for (unsigned int i = 0 ; i < points.size() ; i++)
{
for(auto &enrichmentObjId : vec)
{
double val = enrichmentObjId->getValue(points[i], origin);
ret[i] = ret[i] + std::abs(val * val * (*d_vext)(points[i]));
}
}
}
return ret;
}
};
Expand Down Expand Up @@ -402,6 +427,7 @@ int main(int argc, char** argv)
std::vector<std::vector<dftefe::utils::Point>>
atomCoordinatesVecInGrid(0, std::vector<dftefe::utils::Point>(0,dftefe::utils::Point(dim, 0.0)));

atomCoordinatesVecInGrid.push_back(atomCoordinatesVec);
for(unsigned int perturbDim = 0 ; perturbDim < numDimPerturbed ; perturbDim ++)
{
for(unsigned int perturbAtomId = 0 ; perturbAtomId < numAtomPerturbed ; perturbAtomId++ )
Expand Down Expand Up @@ -805,7 +831,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
{
for(unsigned int perturbAtomId = 0 ; perturbAtomId < numAtomPerturbed ; perturbAtomId++ )
{
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId)*4;
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId)*4 + 1;
rootCout << "The energies calculated by perturbation to atom "<< atomIdPerturbed <<
" along Dim " << dimIdPerturbed <<" are: "<< energyInPerturbIds[index] << ", "<<energyInPerturbIds[index+1]
<<", "<<energyInPerturbIds[index+2]<<", "<<energyInPerturbIds[index+3]<<"\n";
Expand All @@ -819,7 +845,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
{
for(unsigned int perturbAtomId = 0 ; perturbAtomId < numAtomPerturbed ; perturbAtomId++ )
{
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId);
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId) + 1;
rootCout << "The force calculated by perturbation to atom "<< atomIdPerturbed <<
" along Dim " << dimIdPerturbed <<" is: "<< force[numAtomPerturbed*perturbDim + perturbAtomId]<<"\n";
}
Expand Down
77 changes: 56 additions & 21 deletions analysis/classicalEnrichmentComparison/Forces/KSDFTOrthoEFE/TestKohnShamDft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,12 +123,19 @@ T readParameter(std::string ParamFile, std::string param, utils::ConditionalOStr
std::vector<double>
operator()(const std::vector<utils::Point> &points) const
{
double ylm00 = atoms::Clm(0, 0) * atoms::Dm(0) * atoms::Plm(0, 0, 1) * atoms::Qm(0, 0);
std::vector<double> ret(0);
ret.resize(points.size());
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = (*this)(points[i]);
}
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
auto vec = d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "density");
for(auto &enrichmentObjId : vec)
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = ret[i] + std::abs(enrichmentObjId->getValue(points[i], origin) * (1/ylm00));
}
}
return ret;
}
};
Expand Down Expand Up @@ -181,10 +188,19 @@ T readParameter(std::string ParamFile, std::string param, utils::ConditionalOStr
{
std::vector<double> ret(0);
ret.resize(points.size());
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = (*this)(points[i]);
}
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
auto vec = d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "vtotal");
for (unsigned int i = 0 ; i < points.size() ; i++)
{
for(auto &enrichmentObjId : vec)
{
ret[i] = ret[i] + std::abs(enrichmentObjId->getValue(points[i], origin) *
((*d_b)(points[i]) + (*d_rho)(points[i])));
}
}
}
return ret;
}
};
Expand Down Expand Up @@ -235,10 +251,19 @@ T readParameter(std::string ParamFile, std::string param, utils::ConditionalOStr
{
std::vector<double> ret(0);
ret.resize(points.size());
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = (*this)(points[i]);
}
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
auto vec = d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "vnuclear");
for (unsigned int i = 0 ; i < points.size() ; i++)
{
for(auto &enrichmentObjId : vec)
{
ret[i] = ret[i] + std::abs(enrichmentObjId->getValue(points[i], origin) *
((*d_b)(points[i])));
}
}
}
return ret;
}
};
Expand Down Expand Up @@ -288,10 +313,19 @@ T readParameter(std::string ParamFile, std::string param, utils::ConditionalOStr
{
std::vector<double> ret(0);
ret.resize(points.size());
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = (*this)(points[i]);
}
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
auto vec = d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "orbital");
for (unsigned int i = 0 ; i < points.size() ; i++)
{
for(auto &enrichmentObjId : vec)
{
double val = enrichmentObjId->getValue(points[i], origin);
ret[i] = ret[i] + std::abs(val * val * (*d_vext)(points[i]));
}
}
}
return ret;
}
};
Expand Down Expand Up @@ -462,6 +496,7 @@ int main(int argc, char** argv)
std::vector<std::vector<dftefe::utils::Point>>
atomCoordinatesVecInGrid(0, std::vector<dftefe::utils::Point>(0,dftefe::utils::Point(dim, 0.0)));

atomCoordinatesVecInGrid.push_back(atomCoordinatesVec);
for(unsigned int perturbDim = 0 ; perturbDim < numDimPerturbed ; perturbDim ++)
{
for(unsigned int perturbAtomId = 0 ; perturbAtomId < numAtomPerturbed ; perturbAtomId++ )
Expand Down Expand Up @@ -1084,7 +1119,7 @@ int main(int argc, char** argv)
atomSphericalDataContainer,
atomPartitionTolerance,
atomSymbolVec,
atomCoordinatesVec,
atomCoordinatesVecInGrid[perturbId],
"vnuclear",
linAlgOpContext,
comm);
Expand Down Expand Up @@ -1123,7 +1158,7 @@ int main(int argc, char** argv)
double,
Host,
dim>>(
atomCoordinatesVec,
atomCoordinatesVecInGrid[perturbId],
atomChargesVec,
smearedChargeRadiusVec,
numElectrons,
Expand Down Expand Up @@ -1186,7 +1221,7 @@ int main(int argc, char** argv)
double,
Host,
dim>>(
atomCoordinatesVec,
atomCoordinatesVecInGrid[perturbId],
atomChargesVec,
smearedChargeRadiusVec,
numElectrons,
Expand Down Expand Up @@ -1241,7 +1276,7 @@ int main(int argc, char** argv)
{
for(unsigned int perturbAtomId = 0 ; perturbAtomId < numAtomPerturbed ; perturbAtomId++ )
{
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId)*4;
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId)*4 + 1;
rootCout << "The energies calculated by perturbation to atom "<< atomIdPerturbed <<
" along Dim " << dimIdPerturbed <<" are: "<< energyInPerturbIds[index] << ", "<<energyInPerturbIds[index+1]
<<", "<<energyInPerturbIds[index+2]<<", "<<energyInPerturbIds[index+3]<<"\n";
Expand All @@ -1255,7 +1290,7 @@ int main(int argc, char** argv)
{
for(unsigned int perturbAtomId = 0 ; perturbAtomId < numAtomPerturbed ; perturbAtomId++ )
{
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId);
unsigned int index = (numAtomPerturbed*perturbDim + perturbAtomId) + 1;
rootCout << "The force calculated by perturbation to atom "<< atomIdPerturbed <<
" along Dim " << dimIdPerturbed <<" is: "<< force[numAtomPerturbed*perturbDim + perturbAtomId]<<"\n";
}
Expand Down
1 change: 0 additions & 1 deletion analysis/classicalEnrichmentComparison/Forces/KSDFTOrthoEFE/param.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@ smearingTemperature = 500.
fermiEnergyTolerance = 1e-8
fracOccupancyTolerance = 1e-3
eigenSolveResidualTolerance = 1e-3
chebyshevPolynomialDegree = 500
maxChebyshevFilterPass = 20
numWantedEigenvalues = 20
scfDensityResidualNormTolerance = 1e-5
Expand Down
8 changes: 4 additions & 4 deletions analysis/classicalEnrichmentComparison/H.xml
View file

Large diffs are not rendered by default.

30 changes: 23 additions & 7 deletions analysis/classicalEnrichmentComparison/KSDFTClassical/TestKohnShamDft.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,7 @@ double rho1sOrbital(const dftefe::utils::Point &point, const std::vector<dftefe:
std::vector<std::string> d_atomSymbolVec;
std::vector<utils::Point> d_atomCoordinatesVec;
std::vector<double> d_atomChargesVec;
double d_ylm00;

public:
RhoFunction(
Expand All @@ -117,28 +118,27 @@ double rho1sOrbital(const dftefe::utils::Point &point, const std::vector<dftefe:
, d_atomSymbolVec(atomSymbol)
, d_atomCoordinatesVec(atomCoordinates)
, d_atomChargesVec(atomCharges)
, d_ylm00(atoms::Clm(0, 0) * atoms::Dm(0) * atoms::Plm(0, 0, 1) * atoms::Qm(0, 0))
{}

double
operator()(const utils::Point &point) const
{
double ylm00 = atoms::Clm(0, 0) * atoms::Dm(0) * atoms::Plm(0, 0, 1) * atoms::Qm(0, 0);
double retValue = 0;
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
{
utils::Point origin(d_atomCoordinatesVec[atomId]);
for(auto &enrichmentObjId :
d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "density"))
{
retValue = retValue + std::abs(enrichmentObjId->getValue(point, origin) * (1/ylm00));
retValue = retValue + std::abs(enrichmentObjId->getValue(point, origin) * (1/d_ylm00));
}
}
return retValue;
}
std::vector<double>
operator()(const std::vector<utils::Point> &points) const
{
double ylm00 = atoms::Clm(0, 0) * atoms::Dm(0) * atoms::Plm(0, 0, 1) * atoms::Qm(0, 0);
std::vector<double> ret(0);
ret.resize(points.size());
for (size_type atomId = 0 ; atomId < d_atomCoordinatesVec.size() ; atomId++)
Expand All @@ -148,7 +148,7 @@ double rho1sOrbital(const dftefe::utils::Point &point, const std::vector<dftefe:
for(auto &enrichmentObjId : vec)
for (unsigned int i = 0 ; i < points.size() ; i++)
{
ret[i] = ret[i] + std::abs(enrichmentObjId->getValue(points[i], origin) * (1/ylm00));
ret[i] = ret[i] + std::abs(enrichmentObjId->getValue(points[i], origin) * (1/d_ylm00));
}
}
return ret;
Expand Down Expand Up @@ -252,8 +252,8 @@ double rho1sOrbital(const dftefe::utils::Point &point, const std::vector<dftefe:
for(auto &enrichmentObjId :
d_atomSphericalDataContainer->getSphericalData(d_atomSymbolVec[atomId], "orbital"))
{
retValue = retValue + std::abs(enrichmentObjId->getValue(point, origin) * enrichmentObjId->getValue(point, origin) *
(*d_vext)(point));
double val = enrichmentObjId->getValue(point, origin);
retValue = retValue + std::abs(val * val * (*d_vext)(point));
}
}
return retValue;
Expand Down Expand Up @@ -460,7 +460,7 @@ int main(int argc, char** argv)

std::vector<double> smearedChargeRadiusVec(atomCoordinatesVec.size(),rc);

// initialize the basis Manager
// initialize the basis DofHandler

std::shared_ptr<const basis::FEBasisDofHandler<double, Host,dim>> basisDofHandlerTotalPot =
std::make_shared<basis::CFEBasisDofHandlerDealii<double, Host,dim>>(triangulationBase, feOrderElec, comm);
Expand Down Expand Up @@ -625,6 +625,7 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
maxRecursion);

unsigned int nQuad = quadRuleContainerAdaptiveElec->nQuadraturePoints();
unsigned int nQuadMax = nQuad;
int mpierr = utils::mpi::MPIAllreduce<Host>(
utils::mpi::MPIInPlace,
&nQuad,
Expand All @@ -633,6 +634,14 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
utils::mpi::MPISum,
comm);

int mpierr = utils::mpi::MPIAllreduce<Host>(
utils::mpi::MPIInPlace,
&nQuadMax,
1,
utils::mpi::Types<size_type>::getMPIDatatype(),
utils::mpi::MPIMax,
comm);
rootCout << "Maximum Number of quadrature points in a processor: "<< nQuadMax<<"\n";
rootCout << "Number of quadrature points in adaptive quadrature: "<< nQuad<<"\n";

basisAttrMap[basis::BasisStorageAttributes::StoreValues] = true;
Expand Down Expand Up @@ -670,6 +679,13 @@ std::shared_ptr<linearAlgebra::OperatorContext<double,
*feBDElectrostaticsHamiltonian,
numWantedEigenvalues * ksdft::KSDFTDefaults::CELL_BATCH_SIZE);

basisAttrMap[basis::BasisStorageAttributes::StoreValues] = true;
basisAttrMap[basis::BasisStorageAttributes::StoreGradient] = false;
basisAttrMap[basis::BasisStorageAttributes::StoreHessian] = false;
basisAttrMap[basis::BasisStorageAttributes::StoreOverlap] = false;
basisAttrMap[basis::BasisStorageAttributes::StoreGradNiGradNj] = false;
basisAttrMap[basis::BasisStorageAttributes::StoreJxW] = true;

std::shared_ptr<basis::FEBasisDataStorage<double, Host>> feBDTotalChargeRhs =
std::make_shared<basis::CFEBasisDataStorageDealii<double, double, Host,dim>>
(basisDofHandlerTotalPot, quadAttrAdaptive, basisAttrMap);
Expand Down
3 changes: 1 addition & 2 deletions analysis/classicalEnrichmentComparison/KSDFTClassical/param.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ smearingTemperature = 500.
fermiEnergyTolerance = 1e-8
fracOccupancyTolerance = 1e-3
eigenSolveResidualTolerance = 1e-3
chebyshevPolynomialDegree = 1250
maxChebyshevFilterPass = 40
numWantedEigenvalues = 20
scfDensityResidualNormTolerance = 1e-5
Expand All @@ -23,4 +22,4 @@ isAdaptiveAndersonMixingParameter = 0
evaluateEnergyEverySCF = 1
num1DGaussSizeVCorrecPlusPhi = 6
isNumericalNuclearSolve = 1
num1DGaussSizeSmearNucl = 6
num1DGaussSizeSmearNucl = 6
Loading

0 comments on commit 20f7f6e

Please sign in to comment.