Merge branch 'Rob-green_spherical' into Rob-green_metalnp

Author: Roberto Di Remigio

src/bin/plot_green_spherical.cpp

@@ -7,24 +7,56 @@
 #include <Eigen/Core>
 
 #include "CollocationIntegrator.hpp"
+#include "DerivativeTypes.hpp"
 #include "GeneralPurpose.hpp"
 #include "LoggerInterface.hpp"
 #include "OneLayerTanh.hpp"
+#include "PhysicalConstants.hpp"
 #include "SphericalDiffuse.hpp"
 #include "UniformDielectric.hpp"
 
-std::vector<double> generate_grid(double, double, size_t);
+/*! \brief Generate uniform grid of points in interval
+ *  \param[in] min_val lower bound of interval
+ *  \param[in] max_val upper bound of interval
+ *  \param[in] nPoints number of points in interval
+ */
+std::vector<double> generate_grid(double min_val, double max_val, size_t nPoints);
+
+/*! \brief Plot Green's function
+ *  The dielectric sphere is centered in (25.0, 25.0, 0.0) with radius 10.0
+ *  and width 5.0 The dielectric profile is parametrized by a tanh sigmoid
+ *  with epsInside = 78.39 (water) and epsOutside = 1.0 (vacuum)
+ *  The maximum value of angular momentum is set to 150
+ *  All distances in Angstrom.
+ *  The source charge is at (26.0, 25.0, 0.0) i.e. inside the dielectric sphere
+ */
 void case1();
+
+/*! \brief Plot Green's function
+ *  The dielectric sphere is centered in (25.0, 25.0, 0.0) with radius 10.0
+ *  and width 5.0 The dielectric profile is parametrized by a tanh sigmoid
+ *  with epsInside = 78.39 (water) and epsOutside = 1.0 (vacuum)
+ *  The maximum value of angular momentum is set to 150
+ *  All distances in Angstrom.
+ *  The source charge is at (1.0, 1.0, 0.0) i.e. outside the dielectric sphere
+ */
 void case2();
+
+/*! \brief Plot Green's function
+ *  The dielectric sphere is centered in (25.0, 25.0, 0.0) with radius 10.0
+ *  and width 5.0 The dielectric profile is parametrized by a tanh sigmoid
+ *  with epsInside = 78.39 (water) and epsOutside = 1.0 (vacuum)
+ *  The maximum value of angular momentum is set to 150
+ *  All distances in Angstrom.
+ *  The source charge is at (15.0, 25.0, 0.0) i.e. in the transition layer
+ */
 void case3();
-void case4();
 
 int main()
 {
     case1();
-//    case2();
-//    case3();
-//    case4();
+    case2();
+    case3();
 
     return EXIT_SUCCESS;
 }
@@ -42,76 +74,167 @@ std::vector<double> generate_grid(double min_val, double max_val, size_t nPoints
 void case1()
 {
     size_t nPoints = 100;
-    double epsInside = 80.0;
-    double epsOutside = 2.0;
+    double epsInside = 78.39;
+    double epsOutside = 1.0;
+    double angToBohr = angstromToBohr(2010);
     Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
-    double sphereRadius = 10.0;
-    double width = 5.0;
+    sphereCenter << 25.0, 25.0, 0.0;
+    sphereCenter *= angToBohr;
+    double sphereRadius = 10.0 * angToBohr;
+    double width = 5.0 * angToBohr;
     int maxL = 150;
 
     Eigen::Vector3d source = Eigen::Vector3d::Zero();
+    source << 26.0, 25.0, 0.0;
+    source *= angToBohr;
     Eigen::Vector3d probe = Eigen::Vector3d::Zero();
 
-    double min = -50.0;
-    double max = 50.0;
+    // Conversion is done later, so that x, y values are printed in Angstrom
+    double min = 0.0;
+    double max = 40.0;
     std::vector<double> grid = generate_grid(min, max, nPoints);
     double delta = 0.1;
 
     std::ofstream out;
 
     SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
     LOG(gf);
+    gf.toFile("CASE1");
+    UniformDielectric<AD_directional, CollocationIntegrator> gf_i(epsInside);
+    UniformDielectric<AD_directional, CollocationIntegrator> gf_o(epsOutside);
 
     out.open("gf_spherical_CASE1.dat");
     out.precision(16);
     for (size_t i = 0; i < nPoints; ++i) {
-	source << 1.0, 1.0, 0.0;
         for (size_t j = 0; j < nPoints; ++j) {
-            probe << (grid[i] + delta), (grid[j] + delta), 0.0;
+            probe << (grid[i] + delta), (grid[j] + delta), delta;
+	    probe *= angToBohr;
             out << '\t' << i
-		        << '\t' << j
-		        << '\t' << probe(0)
-                << '\t' << probe(1)
+		<< '\t' << j
+		<< '\t' << grid[i] + delta
+                << '\t' << grid[j] + delta
                 << '\t' << gf.Coulomb(source, probe)
                 << '\t' << gf.imagePotential(source, probe)
                 << '\t' << gf.kernelS(source, probe)
-		        << '\t' << (1.0 / gf.coefficientCoulomb(source, probe)) << std::endl;
+		<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe))
+		<< '\t' << gf_i.kernelS(source, probe)
+		<< '\t' << gf_o.kernelS(source, probe)
+		<< std::endl;
         }
+	out << std::endl;
     }
     out.close();
     LOG_TIME;
+}
 
-    /*
-    UniformDielectric<AD_directional, CollocationIntegrator> gf_inside(epsInside);
-    out.open("gf_uniform_inside_CASE1.dat");
-    out << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
-    out << "gf_value" << '\t' << "derivativeProbe" << std::endl;
+void case2()
+{
+    size_t nPoints = 100;
+    double epsInside = 78.39;
+    double epsOutside = 1.0;
+    double angToBohr = angstromToBohr(2010);
+    Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
+    sphereCenter << 25.0, 25.0, 0.0;
+    sphereCenter *= angToBohr;
+    double sphereRadius = 10.0 * angToBohr;
+    double width = 5.0 * angToBohr;
+    int maxL = 150;
+
+    Eigen::Vector3d source = Eigen::Vector3d::Zero();
+    source << 1.0, 1.0, 0.0;
+    source *= angToBohr;
+    Eigen::Vector3d probe = Eigen::Vector3d::Zero();
+
+    // Conversion is done later, so that x, y values are printed in Angstrom
+    double min = 0.0;
+    double max = 40.0;
+    std::vector<double> grid = generate_grid(min, max, nPoints);
+    double delta = 0.1;
+
+    std::ofstream out;
+
+    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
+    LOG(gf);
+    gf.toFile("CASE2");
+    UniformDielectric<AD_directional, CollocationIntegrator> gf_i(epsInside);
+    UniformDielectric<AD_directional, CollocationIntegrator> gf_o(epsOutside);
+
+    out.open("gf_spherical_CASE2.dat");
+    out.precision(16);
     for (size_t i = 0; i < nPoints; ++i) {
         for (size_t j = 0; j < nPoints; ++j) {
-            probe << x[i], 0.0, z[j];
-            out << '\t' << x[i]
-                << '\t' << z[j]
-                << '\t' << gf_inside.kernelS(source, probe)
-                << '\t' << gf_inside.derivativeProbe(Eigen::Vector3d::UnitZ(), source, probe) << std::endl;
+            probe << (grid[i] + delta), (grid[j] + delta), delta;
+	    probe *= angToBohr;
+            out << '\t' << i
+		<< '\t' << j
+		<< '\t' << grid[i] + delta
+                << '\t' << grid[j] + delta
+                << '\t' << gf.Coulomb(source, probe)
+                << '\t' << gf.imagePotential(source, probe)
+                << '\t' << gf.kernelS(source, probe)
+		<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe))
+		<< '\t' << gf_i.kernelS(source, probe)
+		<< '\t' << gf_o.kernelS(source, probe)
+		<< std::endl;
         }
+	out << std::endl;
     }
-    out.precision(16);
     out.close();
+    LOG_TIME;
+}
+
+void case3()
+{
+    size_t nPoints = 100;
+    double epsInside = 78.39;
+    double epsOutside = 1.0;
+    double angToBohr = angstromToBohr(2010);
+    Eigen::Vector3d sphereCenter = Eigen::Vector3d::Zero();
+    sphereCenter << 25.0, 25.0, 0.0;
+    sphereCenter *= angToBohr;
+    double sphereRadius = 10.0 * angToBohr;
+    double width = 5.0 * angToBohr;
+    int maxL = 150;
+
+    Eigen::Vector3d source = Eigen::Vector3d::Zero();
+    source << 15.0, 25.0, 0.0;
+    source *= angToBohr;
+    Eigen::Vector3d probe = Eigen::Vector3d::Zero();
+
+    // Conversion is done later, so that x, y values are printed in Angstrom
+    double min = 0.0;
+    double max = 40.0;
+    std::vector<double> grid = generate_grid(min, max, nPoints);
+    double delta = 0.1;
+
+    std::ofstream out;
+
+    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
+    LOG(gf);
+    gf.toFile("CASE3");
+    UniformDielectric<AD_directional, CollocationIntegrator> gf_i(epsInside);
+    UniformDielectric<AD_directional, CollocationIntegrator> gf_o(epsOutside);
 
-    UniformDielectric<AD_directional, CollocationIntegrator> gf_outside(epsOutside);
-    out.open("gf_uniform_outside_CASE1.dat");
-    out << "#" << '\t' << "x_2" << '\t' << "z_2" << '\t';
-    out << "gf_value" << '\t' << "derivativeProbe" << std::endl;
+    out.open("gf_spherical_CASE3.dat");
     out.precision(16);
     for (size_t i = 0; i < nPoints; ++i) {
         for (size_t j = 0; j < nPoints; ++j) {
-            probe << x[i], 0.0, z[j];
-            out << '\t' << x[i]
-                << '\t' << z[j]
-                << '\t' << gf_outside.kernelS(source, probe)
-                << '\t' << gf_outside.derivativeProbe(Eigen::Vector3d::UnitZ(), source, probe) << std::endl;
+            probe << (grid[i] + delta), (grid[j] + delta), delta;
+	    probe *= angToBohr;
+            out << '\t' << i
+		<< '\t' << j
+		<< '\t' << grid[i] + delta
+                << '\t' << grid[j] + delta
+                << '\t' << gf.Coulomb(source, probe)
+                << '\t' << gf.imagePotential(source, probe)
+                << '\t' << gf.kernelS(source, probe)
+		<< '\t' << (1.0 / gf.coefficientCoulomb(source, probe))
+		<< '\t' << gf_i.kernelS(source, probe)
+		<< '\t' << gf_o.kernelS(source, probe)
+		<< std::endl;
         }
+	out << std::endl;
     }
     out.close();
-    */
+    LOG_TIME;
 }

src/green/SphericalDiffuse.hpp

@@ -191,6 +191,16 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
     std::tuple<double, double> epsilon(const Eigen::Vector3d & point) const {
         return this->profile_((point + this->origin_).norm());
     }
+    void toFile(const std::string & prefix = "") {
+	std::string tmp;
+	prefix.empty() ? tmp = prefix : tmp = prefix + "-";
+        writeToFile(zetaC_, tmp + "zetaC.dat");
+        writeToFile(omegaC_, tmp + "omegaC.dat");
+	for (int L = 1; L <= maxLGreen_; ++L) {
+	    writeToFile(zeta_[L], tmp + "zeta_" + std::to_string(L) + ".dat");
+	    writeToFile(omega_[L], tmp + "omega_" + std::to_string(L) + ".dat");
+	}
+    }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW /* See http://eigen.tuxfamily.org/dox/group__TopicStructHavingEigenMembers.html */
 private:
     using StateType = interfaces::StateType;
@@ -274,14 +284,12 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
         LOG("Computing first radial solution L = " + std::to_string(maxLC_));
         timerON("computeZeta for coefficient");
         zetaC_ = RadialFunction<StateType, LnTransformedRadial, Zeta>(maxLC_, r_0_, r_infinity_, eval_, params_);
-        writeToFile(zetaC_, "zetaC.dat");
         timerOFF("computeZeta for coefficient");
         LOG("DONE: Computing first radial solution L = " + std::to_string(maxLC_));
 
         LOG("Computing second radial solution L = " + std::to_string(maxLC_));
         timerON("computeOmega for coefficient");
         omegaC_ = RadialFunction<StateType, LnTransformedRadial, Omega>(maxLC_, r_0_, r_infinity_, eval_, params_);
-        writeToFile(omegaC_, "omegaC.dat");
         timerOFF("computeOmega for coefficient");
         LOG("Computing second radial solution L = " + std::to_string(maxLC_));
         LOG("DONE: Computing coefficient for the separation of the Coulomb singularity");
@@ -369,10 +377,14 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
         double gr12 = 0.0;
         if (r1 < r2) {
             gr12 = std::exp(zeta1 - zeta2) * (2*L +1) / denominator;
-            gr12 = (gr12 - std::pow(r1/r2, L) / (r2 * Cr12) ) * pl_x ;
+	    double f_L = r1 / r2;
+	    for (int i = 1; i < L; ++i) { f_L *= r1 / r2; }
+            gr12 = (gr12 - f_L / (r2 * Cr12) ) * pl_x ;
         } else {
             gr12 = std::exp(omega1 - omega2) * (2*L +1) / denominator;
-            gr12 = (gr12 - std::pow(r2/r1, L) / (r1 * Cr12) ) * pl_x ;
+	    double f_L = r2 / r1;
+	    for (int i = 1; i < L; ++i) { f_L *= r2 / r1; }
+            gr12 = (gr12 - f_L / (r1 * Cr12) ) * pl_x ;
         }
 
         return gr12;
@@ -422,11 +434,15 @@ class SphericalDiffuse : public GreensFunction<Numerical, IntegratorPolicy, Prof
         double denominator = (d_zeta2 - d_omega2) * std::pow(r2, 2) * eps_r2;
 
         if (r1 < r2) {
+	    double f_L = r1 / r2;
+	    for (int i = 1; i < maxLC_; ++i) { f_L *= r1 / r2; }
             tmp = std::exp(zeta1 - zeta2) * (2*maxLC_ +1) / denominator;
-            coeff = std::pow(r1/r2, maxLC_) / (tmp * r2);
+            coeff = f_L / (tmp * r2);
         } else {
+	    double f_L = r2 / r1;
+	    for (int i = 1; i < maxLC_; ++i) { f_L *= r2 / r1; }
             tmp = std::exp(omega1 - omega2) * (2*maxLC_ +1) / denominator;
-            coeff = std::pow(r2/r1, maxLC_) / (tmp * r1);
+            coeff = f_L / (tmp * r1);
         }
 
         return coeff;

tests/bi_operators/bi_operators_collocation.cpp

@@ -52,17 +52,10 @@
 struct CollocationIntegratorTest {
     double radius;
     double epsilon;
-    Eigen::Vector3d source, probe, sourceNormal, probeNormal;
     GePolCavity cavity;
     CollocationIntegratorTest() { SetUp(); }
     void SetUp() {
 	    epsilon = 80.0;
-        source = Eigen::Vector3d::Random();
-        sourceNormal = source + Eigen::Vector3d::Random();
-        sourceNormal.normalize();
-        probe = Eigen::Vector3d::Random();
-        probeNormal = probe + Eigen::Vector3d::Random();
-        probeNormal.normalize();
 
 	    Molecule molec = dummy<0>(1.44 / convertBohrToAngstrom);
         double area = 10.0;

tests/bi_operators/bi_operators_numerical.cpp

@@ -59,7 +59,6 @@ struct NumericalIntegratorTest {
     Eigen::Vector3d euler;
     double eps;
     double kappa;
-    Eigen::Vector3d source, probe, sourceNormal, probeNormal;
     GePolCavity cavity;
     NumericalIntegratorTest() { SetUp(); }
     void SetUp() {
@@ -69,12 +68,6 @@ struct NumericalIntegratorTest {
         euler << 0.0, 0.0, 0.0;
         eps = 80.0;
         kappa = 1.5;
-        source = Eigen::Vector3d::Random();
-        sourceNormal = source + Eigen::Vector3d::Random();
-        sourceNormal.normalize();
-        probe = Eigen::Vector3d::Random();
-        probeNormal = probe + Eigen::Vector3d::Random();
-        probeNormal.normalize();
 
         Molecule molec = dummy<0>(1.44 / convertBohrToAngstrom);
         double area = 10.0;