From 107eb5a6b77d074882e839685462c1446919b873 Mon Sep 17 00:00:00 2001 From: Eduard Valeyev Date: Thu, 8 Feb 2024 09:01:50 -0500 Subject: [PATCH] clean up NuclearDensityFunctor --- src/apps/moldft/preal.cc | 21 ------- src/apps/moldft/testperiodicdft.cc | 2 +- src/madness/chem/CMakeLists.txt | 1 + src/madness/chem/potentialmanager.cc | 93 ++++++++++++++++++++++++++++ src/madness/chem/potentialmanager.h | 64 +++++++++---------- 5 files changed, 124 insertions(+), 57 deletions(-) create mode 100644 src/madness/chem/potentialmanager.cc diff --git a/src/apps/moldft/preal.cc b/src/apps/moldft/preal.cc index c365f2158eb..f7d491942e6 100644 --- a/src/apps/moldft/preal.cc +++ b/src/apps/moldft/preal.cc @@ -96,27 +96,6 @@ class AtomicBasisFunctor : public FunctionFunctorInterface { } }; -// Functor for the nuclear charge density -class NuclearDensityFunctor : public FunctionFunctorInterface { - Molecule molecule; - std::vector specialpts; -public: - NuclearDensityFunctor(const Molecule& molecule) : - molecule(molecule), specialpts(molecule.get_all_coords_vec()) {} - - double operator()(const Vector& r) const { - return molecule.mol_nuclear_charge_density(r[0], r[1], r[2]); - } - - std::vector special_points() const{ - return specialpts; - } - - Level special_level() { - return 15; - } -}; - // Functor for the nuclear potential class MolecularPotentialFunctor : public FunctionFunctorInterface { private: diff --git a/src/apps/moldft/testperiodicdft.cc b/src/apps/moldft/testperiodicdft.cc index be73d367a2a..cc0fc1bae94 100644 --- a/src/apps/moldft/testperiodicdft.cc +++ b/src/apps/moldft/testperiodicdft.cc @@ -940,7 +940,7 @@ int main(int argc, char** argv) { aobasis.read_file("sto-3g"); // Nuclear potential - real_function_3d vnuc = real_factory_3d(world).functor(real_functor_3d(new NuclearDensityFunctor(molecule, L))).truncate_mode(0).truncate_on_project(); + real_function_3d vnuc = real_factory_3d(world).functor(real_functor_3d(new NuclearDensityFunctor(molecule))).truncate_mode(0).truncate_on_project(); double nuclear_charge=vnuc.trace(); if (world.rank() == 0) print("total nuclear charge", nuclear_charge); vnuc = -1.0*make_coulomb_potential(world, vnuc); diff --git a/src/madness/chem/CMakeLists.txt b/src/madness/chem/CMakeLists.txt index fd0553707b2..a4715b97f1c 100644 --- a/src/madness/chem/CMakeLists.txt +++ b/src/madness/chem/CMakeLists.txt @@ -88,6 +88,7 @@ set(MADCHEM_SOURCES oep.cc pcm.cc pointgroupsymmetry.cc + potentialmanager.cc polynomial.cc SCF.cc SCFOperators.cc diff --git a/src/madness/chem/potentialmanager.cc b/src/madness/chem/potentialmanager.cc new file mode 100644 index 00000000000..7eb707e40d7 --- /dev/null +++ b/src/madness/chem/potentialmanager.cc @@ -0,0 +1,93 @@ +/* + This file is part of MADNESS. + + Copyright (C) 2007,2010 Oak Ridge National Laboratory + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + For more information please contact: + + Robert J. Harrison + Oak Ridge National Laboratory + One Bethel Valley Road + P.O. Box 2008, MS-6367 + + email: harrisonrj@ornl.gov + tel: 865-241-3937 + fax: 865-572-0680 + + $Id$ +*/ + +/// \file moldft/potentialmanager.cc +/// \brief Definition of commonly used density/potential related classes and functions + +#include + +namespace madness { + +NuclearDensityFunctor::NuclearDensityFunctor( + const Molecule &atoms, + const BoundaryConditions<3> &bc, + const Tensor &cell, + int special_level, double rscale) + : atoms(atoms), bc_(bc), cell(cell), special_level_(special_level), + special_points_(this->atoms.get_all_coords_vec()), rscale(rscale) { + if (bc_.is_periodic_any()) { + MADNESS_ASSERT(cell.ndim() == 2 && cell.dim(0) == 3 && cell.dim(1) == 2); + this->maxR = 1; + } else { + this->maxR = 0; + } +} + +double NuclearDensityFunctor::operator()(const coord_3d &x) const { + const auto tol = + 6.0*atoms.smallest_length_scale(); + double sum = 0.0; + for (int i=-maxR; i<=+maxR; i++) { + const double ext_x = cell(0,1) - cell(0,0); + double xx = x[0]+i*ext_x; + if (xx < cell(0,1) + tol && xx > cell(0,0) - tol) { + for (int j=-maxR; j<=+maxR; j++) { + const double ext_y = cell(1,1) - cell(1,0); + double yy = x[1]+j*ext_y; + if (yy < cell(1,1) + tol && yy > cell(1,0) - tol) { + for (int k=-maxR; k<=+maxR; k++) { + const double ext_z = cell(2,1) - cell(2,0); + double zz = x[2]+k*ext_z; + if (zz < cell(2,1) + tol && zz > cell(2,0) - tol) { + sum += atoms.nuclear_charge_density( + xx, yy, zz, rscale); + } + } + } + } + } + } + return sum; +} + +std::vector NuclearDensityFunctor::special_points() const {return special_points_;} + +Level NuclearDensityFunctor::special_level() { + return special_level_; +} + +NuclearDensityFunctor &NuclearDensityFunctor::set_rscale(double rscale) { + this->rscale = rscale; + return *this; +} + +} // namespace madness diff --git a/src/madness/chem/potentialmanager.h b/src/madness/chem/potentialmanager.h index 22ad158783c..5603ff0b695 100644 --- a/src/madness/chem/potentialmanager.h +++ b/src/madness/chem/potentialmanager.h @@ -106,46 +106,40 @@ class CoreOrbitalDerivativeFunctor : public FunctionFunctorInterface { }; }; +/// Default functor for the nuclear charge density + +/// This assumes the default nuclear model optimized to produce potential +/// close that of a point nucleus (smoothed Coulomb potential). The model +/// is class NuclearDensityFunctor : public FunctionFunctorInterface { private: - const Molecule& molecule; - const double R; // Needed only for the periodic case - std::vector specialpt; - const int maxR; + const Molecule& atoms; + BoundaryConditions<3> bc_; + Tensor cell; + std::vector special_points_; + int maxR; + int special_level_ = 15; + double rscale = 1.0; public: - explicit NuclearDensityFunctor(const Molecule& molecule) - : molecule(molecule), R(0), specialpt(molecule.get_all_coords_vec()), maxR(0) - {} - NuclearDensityFunctor(const Molecule& molecule, double R) - : molecule(molecule), R(R), specialpt(molecule.get_all_coords_vec()), maxR(1) - {} - - double operator()(const coord_3d& x) const { - double big = 2*R + 6.0*molecule.smallest_length_scale(); - double sum = 0.0; - for (int i=-maxR; i<=+maxR; i++) { - double xx = x[0]+i*R; - if (xx < big && xx > -big) { - for (int j=-maxR; j<=+maxR; j++) { - double yy = x[1]+j*R; - if (yy < big && yy > -big) { - for (int k=-maxR; k<=+maxR; k++) { - double zz = x[2]+k*R; - if (zz < big && zz > -big) - sum += molecule.nuclear_charge_density(x[0]+i*R, x[1]+j*R, x[2]+k*R); - } - } - } - } - } - return sum; - } + /// Generic constructor, can handle open and periodic boundaries + /// \param molecule atoms + /// \param bc boundary conditions + /// \param cell simulation cell (unit cell, if periodic) + /// \param special_level the initial refinement level + /// \param rscale setting `rscale>1` will make a nucleus larger by a factor of \p rscale (in other words, `rcut` is multiplied by the inverse of by this) + NuclearDensityFunctor(const Molecule& atoms, + const BoundaryConditions<3>& bc = FunctionDefaults<3>::get_bc(), + const Tensor& cell = FunctionDefaults<3>::get_cell(), + int special_level = 15, + double rscale = 1.0); + + double operator()(const coord_3d& x) const; + + std::vector special_points() const; - std::vector special_points() const {return specialpt;} + Level special_level(); - Level special_level() { - return 50; - } + NuclearDensityFunctor& set_rscale(double rscale); };