clean up NuclearDensityFunctor

src/apps/moldft/preal.cc

@@ -96,27 +96,6 @@ class AtomicBasisFunctor : public FunctionFunctorInterface<Q,3> {
     }
 };
 
-// Functor for the nuclear charge density
-class NuclearDensityFunctor : public FunctionFunctorInterface<double,3> {
-  Molecule molecule;
-  std::vector<coord_3d> specialpts;
-public:
-  NuclearDensityFunctor(const Molecule& molecule) : 
-    molecule(molecule), specialpts(molecule.get_all_coords_vec()) {}
-
-  double operator()(const Vector<double,3>& r) const {
-    return molecule.mol_nuclear_charge_density(r[0], r[1], r[2]);
-  }
-
-  std::vector<coord_3d> special_points() const{
-    return specialpts;
-  }
-
-  Level special_level() {
-    return 15;
-  }
-};
-
 // Functor for the nuclear potential
 class MolecularPotentialFunctor : public FunctionFunctorInterface<double,3> {
 private:

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

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

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: [email protected]
+  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 <madness/chem/potentialmanager.h>
+
+namespace madness {
+
+NuclearDensityFunctor::NuclearDensityFunctor(
+    const Molecule &atoms,
+    const BoundaryConditions<3> &bc,
+    const Tensor<double> &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<coord_3d> 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

src/madness/chem/potentialmanager.h

@@ -106,46 +106,40 @@ class CoreOrbitalDerivativeFunctor : public FunctionFunctorInterface<double,3> {
     };
 };
 
+/// 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<double,3> {
 private:
-  const Molecule& molecule;
-  const double R; // Needed only for the periodic case
-  std::vector<coord_3d> specialpt;
-  const int maxR;
+  const Molecule& atoms;
+  BoundaryConditions<3> bc_;
+  Tensor<double> cell;
+  std::vector<coord_3d> 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<double>& cell = FunctionDefaults<3>::get_cell(),
+                        int special_level = 15,
+                        double rscale = 1.0);
+
+  double operator()(const coord_3d& x) const;
+
+  std::vector<coord_3d> special_points() const;
 
-  std::vector<coord_3d> special_points() const {return specialpt;}
+  Level special_level();
 
-  Level special_level() {
-    return 50;
-  }
+  NuclearDensityFunctor& set_rscale(double rscale);
 
 };