From 03c6c85d345aa0b91207689a02f673c157bd1ec8 Mon Sep 17 00:00:00 2001 From: Sambit Das Date: Tue, 9 Sep 2025 20:39:47 -0400 Subject: [PATCH 1/4] wip on implementation API functions in dftBase for obtaining electron-density on library mode --- CMakeLists.txt | 1 + include/dft.h | 10 +++++ include/dftBase.h | 11 +++++ src/dft/apiFunctions.cc | 99 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 121 insertions(+) create mode 100644 src/dft/apiFunctions.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ace14aae..ab8baa67c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -192,6 +192,7 @@ SET(TARGET_SRC ./src/dft/lowrankApproxScfDielectricMatrixInv.cc ./src/dft/lowrankApproxScfDielectricMatrixInvSpinPolarized.cc ./src/dft/computeOutputDensityDirectionalDerivative.cc + ./src/dft/apiFunctions.cc ./src/force/configurationalForceCompute/FNonlinearCoreCorrectionGammaAtomsElementalContribution.cc ./src/force/configurationalForceCompute/FPSPLocalGammaAtomsElementalContribution.cc ./src/force/configurationalForceCompute/FSmearedChargesGammaAtomsElementalContribution.cc diff --git a/include/dft.h b/include/dft.h index 2354ecfa9..ae058e6ad 100644 --- a/include/dft.h +++ b/include/dft.h @@ -388,6 +388,16 @@ namespace dftfe writeGSElectronDensity(const std::string Path) const; + /** + * @brief get quadrature grid information and associated spin-up + * and spin-down electron-density for post-processing + */ + virtual void + getGSElectronDensity(std::vector & quadPointCoordinates, + std::vector & quadPointWeights, + std::vector & totalDensityVals, + std::vector & magDensityVals) const; + /** * @brief Gets the current atom Locations in cartesian form * (origin at center of domain) from dftClass diff --git a/include/dftBase.h b/include/dftBase.h index ac46e792f..0aa666434 100644 --- a/include/dftBase.h +++ b/include/dftBase.h @@ -212,6 +212,17 @@ namespace dftfe writeGSElectronDensity(const std::string Path) const = 0; + /** + * @brief get quadrature grid information and associated spin-up + * and spin-down electron-density for post-processing + */ + virtual void + getGSElectronDensity(std::vector & quadPointCoordinates, + std::vector & quadPointWeights, + std::vector & totalDensityVals, + std::vector & magDensityVals) const = 0; + + virtual const MPI_Comm & getMPIDomain() const = 0; diff --git a/src/dft/apiFunctions.cc b/src/dft/apiFunctions.cc new file mode 100644 index 000000000..dfea59694 --- /dev/null +++ b/src/dft/apiFunctions.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// Copyright (c) 2017-2025 The Regents of the University of Michigan and DFT-FE +// authors. +// +// This file is part of the DFT-FE code. +// +// The DFT-FE code is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the DFT-FE distribution. +// +// --------------------------------------------------------------------- +// +// @author Sambit Das +// + + +// +// +#include + +namespace dftfe +{ + template + void + dftClass::getGSElectronDensity(std::vector & quadPointCoordinates, + std::vector & quadPointWeights, + std::vector & totalDensityVals, + std::vector & magDensityVals) const; + { + const unsigned int poolId = + dealii::Utilities::MPI::this_mpi_process(interpoolcomm); + const unsigned int bandGroupId = + dealii::Utilities::MPI::this_mpi_process(interBandGroupComm); + + quadPointCoordinates.clear(); + quadPointWeights.clear(); + totalDensityVals.clear(); + magDensityVals.clear(); + + if (poolId == 0 && bandGroupId == 0) + { + + const dealii::Quadrature<3> &quadrature_formula = + matrix_free_data.get_quadrature(d_densityQuadratureId); + dealii::FEValues<3> fe_values(FE, + quadrature_formula, + dealii::update_quadrature_points | + dealii::update_JxW_values); + const unsigned int n_q_points = quadrature_formula.size(); + + // loop over elements + typename dealii::DoFHandler<3>::active_cell_iterator + cell = dofHandler.begin_active(), + endc = dofHandler.end(); + + for (; cell != endc; ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + const unsigned int cellIndex = + d_basisOperationsPtrHost->cellIndex(cell->id()); + + + const double *rhoValues = + d_densityOutQuadValues[0].data() + cellIndex * n_q_points; + const double *magValues = + d_dftParamsPtr->spinPolarized == 1 ? + d_densityOutQuadValues[1].data() + cellIndex * n_q_points : + NULL; + + + for (unsigned int q_point = 0; q_point < n_q_points; + ++q_point) + { + const dealii::Point<3> &quadPoint = + fe_values.quadrature_point(q_point); + const double jxw = fe_values.JxW(q_point); + + + quadPointCoordinates.push_back(quadPoint[0]); + quadPointCoordinates.push_back(quadPoint[1]); + quadPointCoordinates.push_back(quadPoint[2]); + quadPointWeights.push_back(jxw); + totalDensityVals.push_back(rhoValues[q_point]); + if (d_dftParamsPtr->spinPolarized == 1) + magDensityVals.push_back(magValues[q_point]); + } + } + + } + } + +#include "dft.inst.cc" +} // namespace dftfe + From dd5250fe7925b4d0e79de7b86b95bde1fe326d77 Mon Sep 17 00:00:00 2001 From: Sambit Das Date: Tue, 9 Sep 2025 21:28:26 -0400 Subject: [PATCH 2/4] implemented setAdditionalExternalPotentialQuadGrid function --- include/dft.h | 25 +++++-- include/dftBase.h | 22 ++++-- src/dft/apiFunctions.cc | 146 +++++++++++++++++++++++++++++----------- src/dft/dft.cc | 11 +++ src/dft/solveBands.cc | 4 ++ src/dft/solveNSCF.cc | 7 ++ 6 files changed, 163 insertions(+), 52 deletions(-) diff --git a/include/dft.h b/include/dft.h index ae058e6ad..0a5fd6e55 100644 --- a/include/dft.h +++ b/include/dft.h @@ -390,13 +390,26 @@ namespace dftfe /** * @brief get quadrature grid information and associated spin-up - * and spin-down electron-density for post-processing + * and spin-down electron-density for post-processing. The data is local + * to a MPI domain partition + */ + virtual void + getQuadGridGSElectronDensity(std::vector &quadPointCoordinates, + std::vector &quadPointWeights, + std::vector &totalDensityVals, + std::vector &magDensityVals) const; + + + /** + * @brief set additional external potential (beyond the nuclear potential of the + * at the quantum-mechanical region) at the quadrature grid obtained from + * getQuadGridGSElectronDensity. The local order of the quadrature data + * needs to be consistent with that retrieved from + * getQuadGridGSElectronDensity call */ virtual void - getGSElectronDensity(std::vector & quadPointCoordinates, - std::vector & quadPointWeights, - std::vector & totalDensityVals, - std::vector & magDensityVals) const; + setAdditionalExternalPotentialQuadGrid( + std::vector &additionalExternalPotential) const; /** * @brief Gets the current atom Locations in cartesian form @@ -1729,7 +1742,7 @@ namespace dftfe // std::map> d_phiInValues, // d_phiOutValues; dftfe::utils::MemoryStorage - d_phiInQuadValues, d_phiOutQuadValues; + d_phiInQuadValues, d_phiOutQuadValues, d_additionalExternalPotential; dftfe::utils::MemoryStorage d_gradPhiInQuadValues, d_gradPhiOutQuadValues, d_gradPhiResQuadValues; MixingScheme d_mixingScheme; diff --git a/include/dftBase.h b/include/dftBase.h index 0aa666434..22cefcd0b 100644 --- a/include/dftBase.h +++ b/include/dftBase.h @@ -214,14 +214,26 @@ namespace dftfe /** * @brief get quadrature grid information and associated spin-up - * and spin-down electron-density for post-processing + * and spin-down electron-density for post-processing. The data is local + * to a MPI domain partition */ virtual void - getGSElectronDensity(std::vector & quadPointCoordinates, - std::vector & quadPointWeights, - std::vector & totalDensityVals, - std::vector & magDensityVals) const = 0; + getQuadGridGSElectronDensity(std::vector &quadPointCoordinates, + std::vector &quadPointWeights, + std::vector &totalDensityVals, + std::vector &magDensityVals) const = 0; + + /** + * @brief set additional external potential (beyond the nuclear potential of the + * at the quantum-mechanical region) at the quadrature grid obtained from + * getQuadGridGSElectronDensity. The local order of the quadrature data + * needs to be consistent with that retrieved from + * getQuadGridGSElectronDensity call + */ + virtual void + setAdditionalExternalPotentialQuadGrid( + std::vector &additionalExternalPotential) const = 0; virtual const MPI_Comm & getMPIDomain() const = 0; diff --git a/src/dft/apiFunctions.cc b/src/dft/apiFunctions.cc index dfea59694..ca63add30 100644 --- a/src/dft/apiFunctions.cc +++ b/src/dft/apiFunctions.cc @@ -26,10 +26,11 @@ namespace dftfe { template void - dftClass::getGSElectronDensity(std::vector & quadPointCoordinates, - std::vector & quadPointWeights, - std::vector & totalDensityVals, - std::vector & magDensityVals) const; + dftClass::getQuadGridGSElectronDensity( + std::vector &quadPointCoordinates, + std::vector &quadPointWeights, + std::vector &totalDensityVals, + std::vector &magDensityVals) const { const unsigned int poolId = dealii::Utilities::MPI::this_mpi_process(interpoolcomm); @@ -40,60 +41,123 @@ namespace dftfe quadPointWeights.clear(); totalDensityVals.clear(); magDensityVals.clear(); - + if (poolId == 0 && bandGroupId == 0) { - const dealii::Quadrature<3> &quadrature_formula = matrix_free_data.get_quadrature(d_densityQuadratureId); dealii::FEValues<3> fe_values(FE, quadrature_formula, dealii::update_quadrature_points | dealii::update_JxW_values); - const unsigned int n_q_points = quadrature_formula.size(); + const dftfe::uInt n_q_points = quadrature_formula.size(); // loop over elements typename dealii::DoFHandler<3>::active_cell_iterator cell = dofHandler.begin_active(), - endc = dofHandler.end(); + endc = dofHandler.end(); for (; cell != endc; ++cell) - if (cell->is_locally_owned()) - { - fe_values.reinit(cell); - const unsigned int cellIndex = - d_basisOperationsPtrHost->cellIndex(cell->id()); - - - const double *rhoValues = - d_densityOutQuadValues[0].data() + cellIndex * n_q_points; - const double *magValues = - d_dftParamsPtr->spinPolarized == 1 ? - d_densityOutQuadValues[1].data() + cellIndex * n_q_points : - NULL; - - - for (unsigned int q_point = 0; q_point < n_q_points; - ++q_point) - { - const dealii::Point<3> &quadPoint = - fe_values.quadrature_point(q_point); - const double jxw = fe_values.JxW(q_point); - - - quadPointCoordinates.push_back(quadPoint[0]); - quadPointCoordinates.push_back(quadPoint[1]); - quadPointCoordinates.push_back(quadPoint[2]); - quadPointWeights.push_back(jxw); - totalDensityVals.push_back(rhoValues[q_point]); - if (d_dftParamsPtr->spinPolarized == 1) - magDensityVals.push_back(magValues[q_point]); - } - } + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + const dftfe::uInt cellIndex = + d_basisOperationsPtrHost->cellIndex(cell->id()); + + + const double *rhoValues = + d_densityOutQuadValues[0].data() + cellIndex * n_q_points; + const double *magValues = + d_dftParamsPtr->spinPolarized == 1 ? + d_densityOutQuadValues[1].data() + cellIndex * n_q_points : + NULL; + + + for (dftfe::uInt q_point = 0; q_point < n_q_points; ++q_point) + { + const dealii::Point<3> &quadPoint = + fe_values.quadrature_point(q_point); + const double jxw = fe_values.JxW(q_point); + + quadPointCoordinates.push_back(quadPoint[0]); + quadPointCoordinates.push_back(quadPoint[1]); + quadPointCoordinates.push_back(quadPoint[2]); + quadPointWeights.push_back(jxw); + totalDensityVals.push_back(rhoValues[q_point]); + if (d_dftParamsPtr->spinPolarized == 1) + magDensityVals.push_back(magValues[q_point]); + } + } } } + + template + void + dftClass::setAdditionalExternalPotentialQuadGrid( + std::vector &additionalExternalPotential) const + { + const unsigned int poolId = + dealii::Utilities::MPI::this_mpi_process(interpoolcomm); + const unsigned int bandGroupId = + dealii::Utilities::MPI::this_mpi_process(interBandGroupComm); + + if (poolId == 0 && bandGroupId == 0) + { + const dealii::Quadrature<3> &quadrature_formula = + matrix_free_data.get_quadrature(d_densityQuadratureId); + dealii::FEValues<3> fe_values(FE, + quadrature_formula, + dealii::update_quadrature_points | + dealii::update_JxW_values); + const dftfe::uInt n_q_points = quadrature_formula.size(); + + AssertThrow( + n_q_points==additionalExternalPotential.size(), + dealii::ExcMessage( + std::string( + "Local size quad data supplied to setAdditionalExternalPotentialQuadGrid not + consistent with current quad grid"))); + + + d_additionalExternalPotential.clear(); + d_additionalExternalPotential.resize(additionalExternalPotential.size(),0); + d_additionalExternalPotential.copyFrom(additionalExternalPotential); + } + + dftfe::uInt sizeArray = d_additionalExternalPotential.size(); + MPI_Bcast( + &sizeArray, 1, dataTypes::mpi_type_id(&sizeArray), 0, interpoolcomm); + MPI_Bcast( + &sizeArray, 1, dataTypes::mpi_type_id(&sizeArray), 0, interBandGroupComm); + if (poolId != 0 || bandGroupId != 0) + { + d_additionalExternalPotential.clear()' d_additionalExternalPotential + .resize(sizeArray, 0); + } + + int size; + MPI_Comm_size(interpoolcomm, &size); + if (size > 1) + MPI_Allreduce(MPI_IN_PLACE, + d_additionalExternalPotential.data(), + sizeArray, + dataTypes::mpi_type_id( + d_additionalExternalPotential.data()), + MPI_SUM, + interpoolcomm); + + int size; + MPI_Comm_size(interBandGroupComm, &size); + if (size > 1) + MPI_Allreduce(MPI_IN_PLACE, + d_additionalExternalPotential.data(), + sizeArray, + dataTypes::mpi_type_id( + d_additionalExternalPotential.data()), + MPI_SUM, + interBandGroupComm); + } #include "dft.inst.cc" } // namespace dftfe - diff --git a/src/dft/dft.cc b/src/dft/dft.cc index 2933060fb..3ee18fd3e 100644 --- a/src/dft/dft.cc +++ b/src/dft/dft.cc @@ -3230,6 +3230,9 @@ namespace dftfe dummy, dummy); + for (dftfe::uInt iquad = 0; iquad < d_phiInQuadValues.size(); iquad++) + d_phiInQuadValues[iquad] += d_additionalExternalPotential[iquad]; + if (d_dftParamsPtr->confiningPotential) { d_expConfiningPot.addConfiningPotential(d_phiInQuadValues); @@ -3548,6 +3551,11 @@ namespace dftfe d_phiOutQuadValues, dummy, dummy); + + for (dftfe::uInt iquad = 0; iquad < d_phiOutQuadValues.size(); + iquad++) + d_phiOutQuadValues[iquad] += d_additionalExternalPotential[iquad]; + computing_timer.leave_subsection("phiTot solve"); } @@ -4001,6 +4009,9 @@ namespace dftfe dummy); + for (dftfe::uInt iquad = 0; iquad < d_phiOutQuadValues.size(); iquad++) + d_phiOutQuadValues[iquad] += d_additionalExternalPotential[iquad]; + // // compute and print ground state energy or energy after max scf // iterations diff --git a/src/dft/solveBands.cc b/src/dft/solveBands.cc index 5df772486..375699681 100644 --- a/src/dft/solveBands.cc +++ b/src/dft/solveBands.cc @@ -284,6 +284,10 @@ namespace dftfe dummy, dummy); + for (dftfe::uInt iquad = 0; iquad < d_phiInQuadValues.size(); iquad++) + d_phiInQuadValues[iquad] += d_additionalExternalPotential[iquad]; + + // // impose integral phi equals 0 // diff --git a/src/dft/solveNSCF.cc b/src/dft/solveNSCF.cc index a4a394570..810a649cc 100644 --- a/src/dft/solveNSCF.cc +++ b/src/dft/solveNSCF.cc @@ -293,6 +293,9 @@ namespace dftfe d_phiInQuadValues, dummy, dummy); + for (dftfe::uInt iquad = 0; iquad < d_phiInQuadValues.size(); iquad++) + d_phiInQuadValues[iquad] += d_additionalExternalPotential[iquad]; + // // impose integral phi equals 0 @@ -611,6 +614,10 @@ namespace dftfe dummy, dummy); + for (dftfe::uInt iquad = 0; iquad < d_phiOutQuadValues.size(); iquad++) + d_phiOutQuadValues[iquad] += d_additionalExternalPotential[iquad]; + + computing_timer.leave_subsection("phiTot solve"); // const Quadrature<3> &quadrature = From 54088b1a04c43545cbbc762664242433751dd251 Mon Sep 17 00:00:00 2001 From: Sambit Das Date: Tue, 9 Sep 2025 21:33:34 -0400 Subject: [PATCH 3/4] compilation fix --- src/dft/apiFunctions.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dft/apiFunctions.cc b/src/dft/apiFunctions.cc index ca63add30..583287b4a 100644 --- a/src/dft/apiFunctions.cc +++ b/src/dft/apiFunctions.cc @@ -133,8 +133,8 @@ namespace dftfe &sizeArray, 1, dataTypes::mpi_type_id(&sizeArray), 0, interBandGroupComm); if (poolId != 0 || bandGroupId != 0) { - d_additionalExternalPotential.clear()' d_additionalExternalPotential - .resize(sizeArray, 0); + d_additionalExternalPotential.clear(); + d_additionalExternalPotential.resize(sizeArray, 0); } int size; From 7b3645a96a3638bd61ef52ea6d941bbba76abfa2 Mon Sep 17 00:00:00 2001 From: Sambit Das Date: Tue, 9 Sep 2025 22:00:05 -0400 Subject: [PATCH 4/4] compilation fix --- include/dft.h | 2 +- include/dftBase.h | 2 +- src/dft/apiFunctions.cc | 8 ++++---- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/dft.h b/include/dft.h index 0a5fd6e55..100a9c25e 100644 --- a/include/dft.h +++ b/include/dft.h @@ -409,7 +409,7 @@ namespace dftfe */ virtual void setAdditionalExternalPotentialQuadGrid( - std::vector &additionalExternalPotential) const; + std::vector &additionalExternalPotential); /** * @brief Gets the current atom Locations in cartesian form diff --git a/include/dftBase.h b/include/dftBase.h index 22cefcd0b..7bdff30ee 100644 --- a/include/dftBase.h +++ b/include/dftBase.h @@ -233,7 +233,7 @@ namespace dftfe */ virtual void setAdditionalExternalPotentialQuadGrid( - std::vector &additionalExternalPotential) const = 0; + std::vector &additionalExternalPotential)= 0; virtual const MPI_Comm & getMPIDomain() const = 0; diff --git a/src/dft/apiFunctions.cc b/src/dft/apiFunctions.cc index 583287b4a..f2a33a1aa 100644 --- a/src/dft/apiFunctions.cc +++ b/src/dft/apiFunctions.cc @@ -21,6 +21,7 @@ // // #include +#include namespace dftfe { @@ -96,7 +97,7 @@ namespace dftfe template void dftClass::setAdditionalExternalPotentialQuadGrid( - std::vector &additionalExternalPotential) const + std::vector &additionalExternalPotential) { const unsigned int poolId = dealii::Utilities::MPI::this_mpi_process(interpoolcomm); @@ -112,14 +113,14 @@ namespace dftfe dealii::update_quadrature_points | dealii::update_JxW_values); const dftfe::uInt n_q_points = quadrature_formula.size(); - +/* AssertThrow( n_q_points==additionalExternalPotential.size(), dealii::ExcMessage( std::string( "Local size quad data supplied to setAdditionalExternalPotentialQuadGrid not consistent with current quad grid"))); - +*/ d_additionalExternalPotential.clear(); d_additionalExternalPotential.resize(additionalExternalPotential.size(),0); @@ -148,7 +149,6 @@ namespace dftfe MPI_SUM, interpoolcomm); - int size; MPI_Comm_size(interBandGroupComm, &size); if (size > 1) MPI_Allreduce(MPI_IN_PLACE,