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..100a9c25e 100644 --- a/include/dft.h +++ b/include/dft.h @@ -388,6 +388,29 @@ 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. 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 + setAdditionalExternalPotentialQuadGrid( + std::vector &additionalExternalPotential); + /** * @brief Gets the current atom Locations in cartesian form * (origin at center of domain) from dftClass @@ -1719,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 ac46e792f..7bdff30ee 100644 --- a/include/dftBase.h +++ b/include/dftBase.h @@ -212,6 +212,29 @@ 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. 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 = 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)= 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..f2a33a1aa --- /dev/null +++ b/src/dft/apiFunctions.cc @@ -0,0 +1,163 @@ +// --------------------------------------------------------------------- +// +// 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 +#include + +namespace dftfe +{ + template + void + 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); + 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 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(); + + for (; cell != endc; ++cell) + 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 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); + + 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 =