[WIP] Tessellationless integration grid #65
Closed
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Conflicts: src/cavity/GePolCavity.cpp src/cavity/GePolCavity.hpp src/pedra/pedra_cavity_interface.F90
Conflicts: doc/gfx/cavity.png doc/gfx/dielectric_profile.png doc/gfx/green.png doc/gfx/tsless.png doc/gfx/utils.png tests/cpcm/cpcm_tsless-NH3.cpp tests/cpcm/cpcm_tsless-point.cpp tests/iefpcm/iefpcm_tsless-NH3.cpp tests/iefpcm/iefpcm_tsless-point.cpp
Code Style CheckCode style violations detected in the following files:
Execute one of the following actions and commit again:
src/cavity/GePolCavity.cpp--- src/cavity/GePolCavity.cpp
+++ src/cavity/GePolCavity.cpp
@@ -81,9 +81,12 @@
int maxsph,
int maxvert) {
- // This is a wrapper for the generatecavity_cpp function defined in the Fortran code PEDRA.
- // Here we allocate the necessary arrays to be passed to PEDRA, in particular we allow
- // for the insertion of additional spheres as in the most general formulation of the
+ // This is a wrapper for the generatecavity_cpp function defined in the Fortran
+ // code PEDRA.
+ // Here we allocate the necessary arrays to be passed to PEDRA, in particular we
+ // allow
+ // for the insertion of additional spheres as in the most general formulation of
+ // the
// GePol algorithm.
double * xtscor = new double[maxts];
@@ -171,38 +174,38 @@
// Go PEDRA, Go!
TIMER_ON("GePolCavity::generatecavity_cpp");
pedra_driver(&maxts,
- &maxsph,
- &maxvert,
- xtscor,
- ytscor,
- ztscor,
- ar,
- xsphcor,
- ysphcor,
- zsphcor,
- rsph,
- &nts,
- &ntsirr,
- &nSpheres_,
- &addedSpheres,
- xe,
- ye,
- ze,
- rin,
- mass,
- &averageArea,
- &probeRadius,
- &minimalRadius,
- &nr_gen,
- &gen1,
- &gen2,
- &gen3,
- nvert,
- vert,
- centr,
- isphe,
- pedra.str().c_str(),
- &len_f_pedra);
+ &maxsph,
+ &maxvert,
+ xtscor,
+ ytscor,
+ ztscor,
+ ar,
+ xsphcor,
+ ysphcor,
+ zsphcor,
+ rsph,
+ &nts,
+ &ntsirr,
+ &nSpheres_,
+ &addedSpheres,
+ xe,
+ ye,
+ ze,
+ rin,
+ mass,
+ &averageArea,
+ &probeRadius,
+ &minimalRadius,
+ &nr_gen,
+ &gen1,
+ &gen2,
+ &gen3,
+ nvert,
+ vert,
+ centr,
+ isphe,
+ pedra.str().c_str(),
+ &len_f_pedra);
TIMER_OFF("GePolCavity::generatecavity_cpp");
// The "intensive" part of updating the spheres related class data members will be
@@ -413,7 +416,7 @@
os << boost::format(" %s") % "Du";
os << boost::format("%9.4f") % (sphereRadius_(idx) * bohrToAngstrom());
os << boost::format(" %1.2f ") % 1.00;
- os << boost::format("%10.6f ") % (sphereCenter_(0, idx) * bohrToAngstrom());
+ os << boost::format("%10.6f ") % (sphereCenter_(0, idx) * bohrToAngstrom());
os << boost::format(" %10.6f ") % (sphereCenter_(1, idx) * bohrToAngstrom());
os << boost::format(" %10.6f ") % (sphereCenter_(2, idx) * bohrToAngstrom());
os << std::endl;
src/cavity/TsLessCavity.cpp--- src/cavity/TsLessCavity.cpp
+++ src/cavity/TsLessCavity.cpp
@@ -41,11 +41,11 @@
namespace pcm {
namespace cavity {
TsLessCavity::TsLessCavity(const Molecule & molec,
- double a,
- double pr,
- double minR,
- double minD,
- int der)
+ double a,
+ double pr,
+ double minR,
+ double minD,
+ int der)
: ICavity(molec),
averageArea_(a),
probeRadius_(pr),
@@ -58,11 +58,11 @@
}
TsLessCavity::TsLessCavity(const Sphere & sph,
- double a,
- double pr,
- double minR,
- double minD,
- int der)
+ double a,
+ double pr,
+ double minR,
+ double minD,
+ int der)
: ICavity(sph),
averageArea_(a),
probeRadius_(pr),
@@ -75,11 +75,11 @@
}
TsLessCavity::TsLessCavity(const std::vector<Sphere> & sph,
- double a,
- double pr,
- double minR,
- double minD,
- int der)
+ double a,
+ double pr,
+ double minR,
+ double minD,
+ int der)
: ICavity(sph),
averageArea_(a),
probeRadius_(pr),
@@ -91,183 +91,219 @@
TIMER_OFF("TsLessCavity::build from list of spheres");
}
-void TsLessCavity::build(int maxts, int maxsph, int maxvert)
-{
- // This is a wrapper for the generatecavity_cpp_ function defined in the Fortran code TsLess.
- // Here we allocate the necessary arrays to be passed to TsLess, in particular we allow
- // for the insertion of additional spheres as in the most general formulation of the
- // GePol algorithm.
-
- int lwork = maxts*maxsph;
- double * xtscor = new double[maxts];
- double * ytscor = new double[maxts];
- double * ztscor = new double[maxts];
- double * ar = new double[maxts];
- double * xsphcor = new double[maxts];
- double * ysphcor = new double[maxts];
- double * zsphcor = new double[maxts];
- double * rsph = new double[maxts];
- double * work = new double[lwork];
-
- // Clean-up possible heap-crap
- std::fill_n(xtscor, maxts, 0.0);
- std::fill_n(ytscor, maxts, 0.0);
- std::fill_n(ztscor, maxts, 0.0);
- std::fill_n(ar, maxts, 0.0);
- std::fill_n(xsphcor, maxts, 0.0);
- std::fill_n(ysphcor, maxts, 0.0);
- std::fill_n(zsphcor, maxts, 0.0);
- std::fill_n(rsph, maxts, 0.0);
- std::fill_n(work, lwork, 0.0);
-
- int nts = 0;
- int ntsirr = 0;
-
- // If there's an overflow in the number of spheres TsLess will die.
- // The maximum number of spheres in PEDRA is set to 200 (primitive+additional)
- // so the integer here declared is just to have enough space C++ side to pass everything back.
- int maxAddedSpheres = 200;
-
- // maximum number of spheres we allow the algorithm to add to our original set.
- // If this number is exceeded, then the algorithm crashes (should look into this...)
- // After the cavity is generated we will update ALL the class data members, both related
- // to spheres and finite elements so that the cavity is fully formed.
-
- Eigen::VectorXd xv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres);
- Eigen::VectorXd yv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres);
- Eigen::VectorXd zv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres);
- Eigen::VectorXd radii_scratch = Eigen::VectorXd::Zero(nSpheres_ +
- maxAddedSpheres); // Not to be confused with the data member inherited from Cavity!!!
-
- for ( int i = 0; i < nSpheres_; ++i ) {
- for ( int j = 0; j < 3; ++j ) {
- xv(i) = sphereCenter_(0, i);
- yv(i) = sphereCenter_(1, i);
- zv(i) = sphereCenter_(2, i);
- }
- radii_scratch(i) = sphereRadius_(i);
+void TsLessCavity::build(int maxts, int maxsph, int maxvert) {
+ // This is a wrapper for the generatecavity_cpp_ function defined in the Fortran
+ // code TsLess.
+ // Here we allocate the necessary arrays to be passed to TsLess, in particular we
+ // allow
+ // for the insertion of additional spheres as in the most general formulation of
+ // the
+ // GePol algorithm.
+
+ int lwork = maxts * maxsph;
+ double * xtscor = new double[maxts];
+ double * ytscor = new double[maxts];
+ double * ztscor = new double[maxts];
+ double * ar = new double[maxts];
+ double * xsphcor = new double[maxts];
+ double * ysphcor = new double[maxts];
+ double * zsphcor = new double[maxts];
+ double * rsph = new double[maxts];
+ double * work = new double[lwork];
+
+ // Clean-up possible heap-crap
+ std::fill_n(xtscor, maxts, 0.0);
+ std::fill_n(ytscor, maxts, 0.0);
+ std::fill_n(ztscor, maxts, 0.0);
+ std::fill_n(ar, maxts, 0.0);
+ std::fill_n(xsphcor, maxts, 0.0);
+ std::fill_n(ysphcor, maxts, 0.0);
+ std::fill_n(zsphcor, maxts, 0.0);
+ std::fill_n(rsph, maxts, 0.0);
+ std::fill_n(work, lwork, 0.0);
+
+ int nts = 0;
+ int ntsirr = 0;
+
+ // If there's an overflow in the number of spheres TsLess will die.
+ // The maximum number of spheres in PEDRA is set to 200 (primitive+additional)
+ // so the integer here declared is just to have enough space C++ side to pass
+ // everything back.
+ int maxAddedSpheres = 200;
+
+ // maximum number of spheres we allow the algorithm to add to our original set.
+ // If this number is exceeded, then the algorithm crashes (should look into
+ // this...)
+ // After the cavity is generated we will update ALL the class data members, both
+ // related
+ // to spheres and finite elements so that the cavity is fully formed.
+
+ Eigen::VectorXd xv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres);
+ Eigen::VectorXd yv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres);
+ Eigen::VectorXd zv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres);
+ Eigen::VectorXd radii_scratch =
+ Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); // Not to be confused with
+ // the data member
+ // inherited from Cavity!!!
+
+ for (int i = 0; i < nSpheres_; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ xv(i) = sphereCenter_(0, i);
+ yv(i) = sphereCenter_(1, i);
+ zv(i) = sphereCenter_(2, i);
}
+ radii_scratch(i) = sphereRadius_(i);
+ }
- double * xe = xv.data();
- double * ye = yv.data();
- double * ze = zv.data();
-
- double *rin = radii_scratch.data();
- double * mass = new double[molecule_.nAtoms()];
- for (size_t i = 0; i < molecule_.nAtoms(); ++i) {
- mass[i] = molecule_.masses(i);
- }
+ double * xe = xv.data();
+ double * ye = yv.data();
+ double * ze = zv.data();
+
+ double * rin = radii_scratch.data();
+ double * mass = new double[molecule_.nAtoms()];
+ for (size_t i = 0; i < molecule_.nAtoms(); ++i) {
+ mass[i] = molecule_.masses(i);
+ }
- addedSpheres = 0;
- // Number of generators and generators of the point group
- int nr_gen = molecule_.pointGroup().nrGenerators();
- int gen1 = molecule_.pointGroup().generators(0);
- int gen2 = molecule_.pointGroup().generators(1);
- int gen3 = molecule_.pointGroup().generators(2);
-
- int weightFunction = 1;
- // Go TsLess, Go!
- TIMER_ON("TsLessCavity::tsless_driver");
- tsless_driver(&maxts, &maxsph, &maxvert, &nSpheres_, &nts, &ntsirr, &addedSpheres,
- xtscor, ytscor, ztscor, ar,
- xsphcor, ysphcor, zsphcor, rsph,
- xe, ye, ze, rin, mass,
- &nr_gen, &gen1, &gen2, &gen3,
- &averageArea_, &minDistance_, &derOrder_, &weightFunction, &probeRadius_, work);
- TIMER_OFF("TsLessCavity::tsless_driver");
-
- // The "intensive" part of updating the spheres related class data members will be of course
- // executed iff addedSpheres != 0
- if (addedSpheres != 0) {
- // Save the number of original spheres
- int orig = nSpheres_;
- // Update the nSpheres
- nSpheres_ += addedSpheres;
- // Resize sphereRadius and sphereCenter...
- sphereRadius_.resize(nSpheres_);
- sphereCenter_.resize(Eigen::NoChange, nSpheres_);
- // Transfer radii and centers.
- // Eigen has no push_back function, so we need to traverse all the spheres...
- sphereRadius_ = radii_scratch.head(nSpheres_);
- for ( int i = 0; i < nSpheres_; ++i ) {
- sphereCenter_(0, i) = xv(i);
- sphereCenter_(1, i) = yv(i);
- sphereCenter_(2, i) = zv(i);
- }
- // Now grow the vector<Sphere> containing the list of spheres
- for ( int i = orig; i < nSpheres_; ++i ) {
- spheres_.push_back(Sphere(sphereCenter_.col(i), sphereRadius_(i)));
- }
+ addedSpheres = 0;
+ // Number of generators and generators of the point group
+ int nr_gen = molecule_.pointGroup().nrGenerators();
+ int gen1 = molecule_.pointGroup().generators(0);
+ int gen2 = molecule_.pointGroup().generators(1);
+ int gen3 = molecule_.pointGroup().generators(2);
+
+ int weightFunction = 1;
+ // Go TsLess, Go!
+ TIMER_ON("TsLessCavity::tsless_driver");
+ tsless_driver(&maxts,
+ &maxsph,
+ &maxvert,
+ &nSpheres_,
+ &nts,
+ &ntsirr,
+ &addedSpheres,
+ xtscor,
+ ytscor,
+ ztscor,
+ ar,
+ xsphcor,
+ ysphcor,
+ zsphcor,
+ rsph,
+ xe,
+ ye,
+ ze,
+ rin,
+ mass,
+ &nr_gen,
+ &gen1,
+ &gen2,
+ &gen3,
+ &averageArea_,
+ &minDistance_,
+ &derOrder_,
+ &weightFunction,
+ &probeRadius_,
+ work);
+ TIMER_OFF("TsLessCavity::tsless_driver");
+
+ // The "intensive" part of updating the spheres related class data members will be
+ // of course
+ // executed iff addedSpheres != 0
+ if (addedSpheres != 0) {
+ // Save the number of original spheres
+ int orig = nSpheres_;
+ // Update the nSpheres
+ nSpheres_ += addedSpheres;
+ // Resize sphereRadius and sphereCenter...
+ sphereRadius_.resize(nSpheres_);
+ sphereCenter_.resize(Eigen::NoChange, nSpheres_);
+ // Transfer radii and centers.
+ // Eigen has no push_back function, so we need to traverse all the spheres...
+ sphereRadius_ = radii_scratch.head(nSpheres_);
+ for (int i = 0; i < nSpheres_; ++i) {
+ sphereCenter_(0, i) = xv(i);
+ sphereCenter_(1, i) = yv(i);
+ sphereCenter_(2, i) = zv(i);
}
-
- nElements_ = static_cast<int>(nts);
- nIrrElements_ = static_cast<int>(ntsirr);
- elementCenter_.resize(Eigen::NoChange, nElements_);
- elementSphereCenter_.resize(Eigen::NoChange, nElements_);
- elementNormal_.resize(Eigen::NoChange, nElements_);
- elementArea_.resize(nElements_);
- elementRadius_.resize(nElements_);
- for( int i = 0; i < nElements_; ++i ) {
- elementCenter_(0,i) = xtscor[i];
- elementCenter_(1,i) = ytscor[i];
- elementCenter_(2,i) = ztscor[i];
- elementArea_(i) = ar[i];
- elementSphereCenter_(0,i) = xsphcor[i];
- elementSphereCenter_(1,i) = ysphcor[i];
- elementSphereCenter_(2,i) = zsphcor[i];
- elementRadius_(i) = rsph[i];
+ // Now grow the vector<Sphere> containing the list of spheres
+ for (int i = orig; i < nSpheres_; ++i) {
+ spheres_.push_back(Sphere(sphereCenter_.col(i), sphereRadius_(i)));
}
+ }
- elementNormal_ = elementCenter_ - elementSphereCenter_;
- for( int i = 0; i < nElements_; ++i) {
- elementNormal_.col(i) /= elementNormal_.col(i).norm();
- }
+ nElements_ = static_cast<int>(nts);
+ nIrrElements_ = static_cast<int>(ntsirr);
+ elementCenter_.resize(Eigen::NoChange, nElements_);
+ elementSphereCenter_.resize(Eigen::NoChange, nElements_);
+ elementNormal_.resize(Eigen::NoChange, nElements_);
+ elementArea_.resize(nElements_);
+ elementRadius_.resize(nElements_);
+ for (int i = 0; i < nElements_; ++i) {
+ elementCenter_(0, i) = xtscor[i];
+ elementCenter_(1, i) = ytscor[i];
+ elementCenter_(2, i) = ztscor[i];
+ elementArea_(i) = ar[i];
+ elementSphereCenter_(0, i) = xsphcor[i];
+ elementSphereCenter_(1, i) = ysphcor[i];
+ elementSphereCenter_(2, i) = zsphcor[i];
+ elementRadius_(i) = rsph[i];
+ }
- // Fill elements_ vector
- for (int i = 0; i < nElements_; ++i) {
- bool irr = false;
- int nv = 1; // TsLess does not generate spherical polygons!!
- // TsLess puts the irreducible tesserae first (? Check with Cris!)
- if (i < nIrrElements_) irr = true;
- Sphere sph(elementSphereCenter_.col(i), elementRadius_(i));
- Eigen::Matrix3Xd vertices, arcs;
- vertices.resize(Eigen::NoChange, nv);
- arcs.resize(Eigen::NoChange, nv);
- // FIXME index of the sphere the element belongs to
- elements_.push_back(Element(nv, 0,
- elementArea_(i),
- elementCenter_.col(i),
- elementNormal_.col(i),
- irr, sph,
- vertices, arcs));
- }
+ elementNormal_ = elementCenter_ - elementSphereCenter_;
+ for (int i = 0; i < nElements_; ++i) {
+ elementNormal_.col(i) /= elementNormal_.col(i).norm();
+ }
- delete[] xtscor;
- delete[] ytscor;
- delete[] ztscor;
- delete[] ar;
- delete[] xsphcor;
- delete[] ysphcor;
- delete[] zsphcor;
- delete[] rsph;
- delete[] work;
- delete[] mass;
+ // Fill elements_ vector
+ for (int i = 0; i < nElements_; ++i) {
+ bool irr = false;
+ int nv = 1; // TsLess does not generate spherical polygons!!
+ // TsLess puts the irreducible tesserae first (? Check with Cris!)
+ if (i < nIrrElements_)
+ irr = true;
+ Sphere sph(elementSphereCenter_.col(i), elementRadius_(i));
+ Eigen::Matrix3Xd vertices, arcs;
+ vertices.resize(Eigen::NoChange, nv);
+ arcs.resize(Eigen::NoChange, nv);
+ // FIXME index of the sphere the element belongs to
+ elements_.push_back(Element(nv,
+ 0,
+ elementArea_(i),
+ elementCenter_.col(i),
+ elementNormal_.col(i),
+ irr,
+ sph,
+ vertices,
+ arcs));
+ }
- built = true;
+ delete[] xtscor;
+ delete[] ytscor;
+ delete[] ztscor;
+ delete[] ar;
+ delete[] xsphcor;
+ delete[] ysphcor;
+ delete[] zsphcor;
+ delete[] rsph;
+ delete[] work;
+ delete[] mass;
+ built = true;
}
-std::ostream & TsLessCavity::printCavity(std::ostream & os)
-{
- os << "Cavity type: TsLess" << std::endl;
- os << "Average point weight = " << averageArea_ << " AU^2" << std::endl;
- os << "Minimal distance between sampling points = " << minDistance_ << " AU" << std::endl;
- os << "Switch function is of class C^" << derOrder_ << std::endl;
- os << "Addition of extra spheres enabled" << std::endl;
- os << "Probe radius = " << probeRadius_ << " AU" << std::endl;
- os << "Number of spheres = " << nSpheres_ << " [initial = " << nSpheres_ -
- addedSpheres << "; added = " << addedSpheres << "]" << std::endl;
- os << "Number of finite elements = " << nElements_;
+std::ostream & TsLessCavity::printCavity(std::ostream & os) {
+ os << "Cavity type: TsLess" << std::endl;
+ os << "Average point weight = " << averageArea_ << " AU^2" << std::endl;
+ os << "Minimal distance between sampling points = " << minDistance_ << " AU"
+ << std::endl;
+ os << "Switch function is of class C^" << derOrder_ << std::endl;
+ os << "Addition of extra spheres enabled" << std::endl;
+ os << "Probe radius = " << probeRadius_ << " AU" << std::endl;
+ os << "Number of spheres = " << nSpheres_
+ << " [initial = " << nSpheres_ - addedSpheres << "; added = " << addedSpheres
+ << "]" << std::endl;
+ os << "Number of finite elements = " << nElements_;
os << "Number of irreducible finite elements = " << nIrrElements_ << std::endl;
os << "============ Spheres list (in Angstrom)" << std::endl;
os << " Sphere on Radius Alpha X Y Z \n";
@@ -291,17 +327,21 @@
os << boost::format(" %s") % "Du";
os << boost::format("%9.4f") % (sphereRadius_(idx) * bohrToAngstrom());
os << boost::format(" %1.2f ") % 1.00;
- os << boost::format("%10.6f ") % (sphereCenter_(0, idx) * bohrToAngstrom());
+ os << boost::format("%10.6f ") % (sphereCenter_(0, idx) * bohrToAngstrom());
os << boost::format(" %10.6f ") % (sphereCenter_(1, idx) * bohrToAngstrom());
os << boost::format(" %10.6f ") % (sphereCenter_(2, idx) * bohrToAngstrom());
os << std::endl;
}
- return os;
+ return os;
}
ICavity * createTsLessCavity(const CavityData & data) {
- return new TsLessCavity(data.molecule, data.area, data.probeRadius, data.minimalRadius,
- data.minDistance, data.derOrder);
+ return new TsLessCavity(data.molecule,
+ data.area,
+ data.probeRadius,
+ data.minimalRadius,
+ data.minDistance,
+ data.derOrder);
}
} // namespace cavity
} // namespace pcm
src/cavity/TsLessCavity.hpp--- src/cavity/TsLessCavity.hpp
+++ src/cavity/TsLessCavity.hpp
@@ -105,13 +105,20 @@
* \param[out] ytscor y-coordinate of tesserae centers (dimension maxts)
* \param[out] ztscor z-coordinate of tesserae centers (dimension maxts)
* \param[out] ar area of the tessera (dimension maxts)
- * \param[out] xsphcor x-coordinate of the sphere center the tessera belongs to (dimension maxts)
- * \param[out] ysphcor y-coordinate of the sphere center the tessera belongs to (dimension maxts)
- * \param[out] zsphcor z-coordinate of the sphere center the tessera belongs to (dimension maxts)
- * \param[out] rsph radii of the sphere the tessera belongs to, i.e. its curvature (dimension maxts)
- * \param[out] xe x-coordinate of the sphere center (dimension nSpheres_ + maxAddedSpheres)
- * \param[out] ye y-coordinate of the sphere center (dimension nSpheres_ + maxAddedSpheres)
- * \param[out] ze z-coordinate of the sphere center (dimension nSpheres_ + maxAddedSpheres)
+ * \param[out] xsphcor x-coordinate of the sphere center the tessera belongs to
+ * (dimension maxts)
+ * \param[out] ysphcor y-coordinate of the sphere center the tessera belongs to
+ * (dimension maxts)
+ * \param[out] zsphcor z-coordinate of the sphere center the tessera belongs to
+ * (dimension maxts)
+ * \param[out] rsph radii of the sphere the tessera belongs to, i.e. its curvature
+ * (dimension maxts)
+ * \param[out] xe x-coordinate of the sphere center (dimension nSpheres_ +
+ * maxAddedSpheres)
+ * \param[out] ye y-coordinate of the sphere center (dimension nSpheres_ +
+ * maxAddedSpheres)
+ * \param[out] ze z-coordinate of the sphere center (dimension nSpheres_ +
+ * maxAddedSpheres)
* \param[out] rin radius of the spheres (dimension nSpheres_ + maxAddedSpheres)
* \param[in] masses atomic masses (for inertia tensor formation in TSLESS)
* \param[in] nr_gen number of symmetry generators
@@ -121,20 +128,44 @@
* \param[in] avgArea average tesserae area
* \param[in] dmin mininal distance between sampling points
* \param[in] nord maximum order of continuous derivative of weight function
- * \param[in] ifun whether to use the normalized or unnormalized form of the weight function
+ * \param[in] ifun whether to use the normalized or unnormalized form of the weight
+ * function
* \param[in] rsolv solvent probe radius
* \param[in] work scratch space
*/
-#define tsless_driver \
- FortranCInterface_MODULE_(tsless_cavity, tsless_driver, TSLESS_CAVITY, TSLESS_DRIVER)
-extern "C" void tsless_driver(int * maxts, int * maxsph, int * maxvert,
- int * nesfp, int * nts, int * ntsirr, int * addsph,
- double * xtscor, double * ytscor, double * ztscor, double * ar,
- double * xsphcor, double * ysphcor, double * zsphcor, double * rsph,
- double * xe, double * ye, double * ze, double * rin, double * masses,
- int * nr_gen, int * gen1, int * gen2, int * gen3,
- double * avgArea, double * dmin, int * nord, int * ifun, double * rsolv,
- double * work);
+#define tsless_driver \
+ FortranCInterface_MODULE_( \
+ tsless_cavity, tsless_driver, TSLESS_CAVITY, TSLESS_DRIVER)
+extern "C" void tsless_driver(int * maxts,
+ int * maxsph,
+ int * maxvert,
+ int * nesfp,
+ int * nts,
+ int * ntsirr,
+ int * addsph,
+ double * xtscor,
+ double * ytscor,
+ double * ztscor,
+ double * ar,
+ double * xsphcor,
+ double * ysphcor,
+ double * zsphcor,
+ double * rsph,
+ double * xe,
+ double * ye,
+ double * ze,
+ double * rin,
+ double * masses,
+ int * nr_gen,
+ int * gen1,
+ int * gen2,
+ int * gen3,
+ double * avgArea,
+ double * dmin,
+ int * nord,
+ int * ifun,
+ double * rsolv,
+ double * work);
ICavity * createTsLessCavity(const CavityData & data);
} // namespace cavity
tests/TestingMolecules.hpp--- tests/TestingMolecules.hpp
+++ tests/TestingMolecules.hpp
@@ -127,34 +127,32 @@
return dummy;
};
-Molecule HF()
-{
- int nAtoms = 2;
-
- Eigen::Vector3d F(0.000000000, 0.00000000, 0.08729478);
- Eigen::Vector3d H(0.000000000, 0.00000000, -1.64558444);
-
- Eigen::MatrixXd geom(3, nAtoms);
- geom.col(0) = F.transpose();
- geom.col(1) = H.transpose();
- Eigen::Vector2d charges, masses;
- charges << 9.0, 1.0;
- masses << 18.9984030, 1.0078250;
- std::vector<Atom> atoms;
- atoms.push_back( Atom("Fluorine", "F", charges(0), masses(0), 2.777897403, F, 1.0) );
- atoms.push_back( Atom("Hydrogen", "H", charges(1), masses(1), 2.267671349, H,
- 1.0) );
-
- std::vector<Sphere> spheres;
- Sphere sph1(F, 2.777897403);
- Sphere sph2(H, 2.267671349);
- spheres.push_back(sph1);
- spheres.push_back(sph2);
+Molecule HF() {
+ int nAtoms = 2;
- // C1
- Symmetry pGroup = buildGroup(0, 0, 0, 0);
+ Eigen::Vector3d F(0.000000000, 0.00000000, 0.08729478);
+ Eigen::Vector3d H(0.000000000, 0.00000000, -1.64558444);
- return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
+ Eigen::MatrixXd geom(3, nAtoms);
+ geom.col(0) = F.transpose();
+ geom.col(1) = H.transpose();
+ Eigen::Vector2d charges, masses;
+ charges << 9.0, 1.0;
+ masses << 18.9984030, 1.0078250;
+ std::vector<Atom> atoms;
+ atoms.push_back(Atom("Fluorine", "F", charges(0), masses(0), 2.777897403, F, 1.0));
+ atoms.push_back(Atom("Hydrogen", "H", charges(1), masses(1), 2.267671349, H, 1.0));
+
+ std::vector<Sphere> spheres;
+ Sphere sph1(F, 2.777897403);
+ Sphere sph2(H, 2.267671349);
+ spheres.push_back(sph1);
+ spheres.push_back(sph2);
+
+ // C1
+ Symmetry pGroup = buildGroup(0, 0, 0, 0);
+
+ return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
};
Molecule NH3() {
tests/cpcm/cpcm_tsless-NH3.cpp--- tests/cpcm/cpcm_tsless-NH3.cpp
+++ tests/cpcm/cpcm_tsless-NH3.cpp
@@ -27,15 +27,15 @@
#include <Eigen/Core>
+#include "TestingMolecules.hpp"
#include "bi_operators/Collocation.hpp"
-#include "green/DerivativeTypes.hpp"
#include "cavity/TsLessCavity.hpp"
-#include "utils/Molecule.hpp"
-#include "green/Vacuum.hpp"
+#include "green/DerivativeTypes.hpp"
#include "green/UniformDielectric.hpp"
+#include "green/Vacuum.hpp"
#include "solver/CPCMSolver.hpp"
+#include "utils/Molecule.hpp"
#include "utils/Symmetry.hpp"
-#include "TestingMolecules.hpp"
using namespace pcm;
using bi_operators::Collocation;
@@ -47,39 +47,40 @@
/*! \class CPCMSolver
* \test \b NH3TsLess tests CPCMSolver using ammonia and a TsLess cavity
*/
-TEST_CASE("Test solver for the C-PCM with NH3 molecule and a TsLess cavity", "[solver][cpcm][cpcm_tsless-NH3]")
-{
- Molecule molec = NH3();
-
- double area = 0.08;
- double minDistance = 0.1;
- double probeRadius = 0.0;
- int derOrder = 8;
- double minRadius = 100.0;
- TsLessCavity cavity = TsLessCavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
-
- double permittivity = 78.39;
- Vacuum<> gfInside;
- UniformDielectric<> gfOutside(permittivity);
- bool symm = true;
- double correction = 0.0;
-
- Collocation S;
-
- CPCMSolver solver(symm, correction);
- solver.buildSystemMatrix(cavity, gfInside, gfOutside, S);
-
- double Ncharge = 7.0;
- double Hcharge = 1.0;
- size_t size = cavity.size();
- Eigen::VectorXd fake_mep = computeMEP(molec, cavity.elements());
- // The total ASC for a conductor is -Q
- // for CPCM it will be -Q*[(epsilon-1)/epsilon + correction]
- Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
- fake_asc = solver.computeCharge(fake_mep);
- double totalASC = - (Ncharge + 3.0 * Hcharge) * (permittivity - 1) /
- (permittivity + correction);
- double totalFakeASC = fake_asc.sum();
- CAPTURE(totalASC - totalFakeASC);
- REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
+TEST_CASE("Test solver for the C-PCM with NH3 molecule and a TsLess cavity",
+ "[solver][cpcm][cpcm_tsless-NH3]") {
+ Molecule molec = NH3();
+
+ double area = 0.08;
+ double minDistance = 0.1;
+ double probeRadius = 0.0;
+ int derOrder = 8;
+ double minRadius = 100.0;
+ TsLessCavity cavity =
+ TsLessCavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
+
+ double permittivity = 78.39;
+ Vacuum<> gfInside;
+ UniformDielectric<> gfOutside(permittivity);
+ bool symm = true;
+ double correction = 0.0;
+
+ Collocation S;
+
+ CPCMSolver solver(symm, correction);
+ solver.buildSystemMatrix(cavity, gfInside, gfOutside, S);
+
+ double Ncharge = 7.0;
+ double Hcharge = 1.0;
+ size_t size = cavity.size();
+ Eigen::VectorXd fake_mep = computeMEP(molec, cavity.elements());
+ // The total ASC for a conductor is -Q
+ // for CPCM it will be -Q*[(epsilon-1)/epsilon + correction]
+ Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
+ fake_asc = solver.computeCharge(fake_mep);
+ double totalASC =
+ -(Ncharge + 3.0 * Hcharge) * (permittivity - 1) / (permittivity + correction);
+ double totalFakeASC = fake_asc.sum();
+ CAPTURE(totalASC - totalFakeASC);
+ REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
}
tests/cpcm/cpcm_tsless-point.cpp--- tests/cpcm/cpcm_tsless-point.cpp
+++ tests/cpcm/cpcm_tsless-point.cpp
@@ -27,13 +27,13 @@
#include <Eigen/Core>
+#include "TestingMolecules.hpp"
#include "bi_operators/Collocation.hpp"
-#include "green/DerivativeTypes.hpp"
#include "cavity/TsLessCavity.hpp"
-#include "green/Vacuum.hpp"
+#include "green/DerivativeTypes.hpp"
#include "green/UniformDielectric.hpp"
+#include "green/Vacuum.hpp"
#include "solver/CPCMSolver.hpp"
-#include "TestingMolecules.hpp"
using namespace pcm;
using bi_operators::Collocation;
@@ -42,101 +42,101 @@
using green::UniformDielectric;
using solver::CPCMSolver;
-SCENARIO("Test solver for the C-PCM for a point charge and a TsLess cavity", "[solver][cpcm][cpcm_tsless-point]")
-{
- GIVEN("An isotropic environment modelled as a perfect conductor and a point charge")
- {
- double permittivity = 78.39;
- Vacuum<> gfInside;
- UniformDielectric<> gfOutside(permittivity);
- bool symm = true;
- double correction = 0.0;
-
- Collocation S;
-
- double charge = 8.0;
- double totalASC = - charge * (permittivity - 1) / (permittivity + correction);
-
- /*! \class CPCMSolver
- * \test \b pointChargeTsLess tests CPCMSolver using a point charge with a TsLess cavity
- * The point charge is at the origin.
- */
- WHEN("the point charge is located at the origin")
- {
- Molecule point = dummy<0>(2.929075493);
- double area = 0.4;
- double probeRadius = 0.0;
- double minRadius = 100.0;
- double minDistance = 0.1;
- int derOrder = 4;
- TsLessCavity cavity = TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
- cavity.saveCavity("tsless_point.npz");
-
- CPCMSolver solver(symm, correction);
- solver.buildSystemMatrix(cavity, gfInside, gfOutside, S);
-
- size_t size = cavity.size();
- Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge);
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_mep(" << i << ") = " << fake_mep(i));
- }
-
- THEN("the apparent surface charge is")
- {
- Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
- fake_asc = solver.computeCharge(fake_mep);
- double totalFakeASC = fake_asc.sum();
-
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_asc(" << i << ") = " << fake_asc(i));
- }
-
- CAPTURE(totalASC);
- CAPTURE(totalFakeASC);
- CAPTURE(totalASC - totalFakeASC);
- REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
- }
+SCENARIO("Test solver for the C-PCM for a point charge and a TsLess cavity",
+ "[solver][cpcm][cpcm_tsless-point]") {
+ GIVEN("An isotropic environment modelled as a perfect conductor and a point "
+ "charge") {
+ double permittivity = 78.39;
+ Vacuum<> gfInside;
+ UniformDielectric<> gfOutside(permittivity);
+ bool symm = true;
+ double correction = 0.0;
+
+ Collocation S;
+
+ double charge = 8.0;
+ double totalASC = -charge * (permittivity - 1) / (permittivity + correction);
+
+ /*! \class CPCMSolver
+ * \test \b pointChargeTsLess tests CPCMSolver using a point charge with a
+ * TsLess cavity
+ * The point charge is at the origin.
+ */
+ WHEN("the point charge is located at the origin") {
+ Molecule point = dummy<0>(2.929075493);
+ double area = 0.4;
+ double probeRadius = 0.0;
+ double minRadius = 100.0;
+ double minDistance = 0.1;
+ int derOrder = 4;
+ TsLessCavity cavity =
+ TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
+ cavity.saveCavity("tsless_point.npz");
+
+ CPCMSolver solver(symm, correction);
+ solver.buildSystemMatrix(cavity, gfInside, gfOutside, S);
+
+ size_t size = cavity.size();
+ Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge);
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_mep(" << i << ") = " << fake_mep(i));
+ }
+
+ THEN("the apparent surface charge is") {
+ Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
+ fake_asc = solver.computeCharge(fake_mep);
+ double totalFakeASC = fake_asc.sum();
+
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_asc(" << i << ") = " << fake_asc(i));
}
- /*! \class CPCMSolver
- * \test \b pointChargeShiftedTsLess tests CPCMSolver using a point charge with a TsLess cavity
- * The point charge is away from the origin.
- */
- AND_WHEN("the point charge is located away from the origin")
- {
- Eigen::Vector3d origin = 100 * Eigen::Vector3d::Random();
- Molecule point = dummy<0>(2.929075493, origin);
- double area = 0.4;
- double probeRadius = 0.0;
- double minRadius = 100.0;
- double minDistance = 0.1;
- int derOrder = 4;
- TsLessCavity cavity = TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
-
- CPCMSolver solver(symm, correction);
- solver.buildSystemMatrix(cavity, gfInside, gfOutside, S);
-
- size_t size = cavity.size();
- Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge, origin);
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_mep(" << i << ") = " << fake_mep(i));
- }
-
- THEN("the surface charge is")
- {
- Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
- fake_asc = solver.computeCharge(fake_mep);
- double totalFakeASC = fake_asc.sum();
-
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_asc(" << i << ") = " << fake_asc(i));
- }
-
- CAPTURE(totalASC);
- CAPTURE(totalFakeASC);
- CAPTURE(totalASC - totalFakeASC);
- REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
- }
+ CAPTURE(totalASC);
+ CAPTURE(totalFakeASC);
+ CAPTURE(totalASC - totalFakeASC);
+ REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
+ }
+ }
+
+ /*! \class CPCMSolver
+ * \test \b pointChargeShiftedTsLess tests CPCMSolver using a point charge with
+ * a TsLess cavity
+ * The point charge is away from the origin.
+ */
+ AND_WHEN("the point charge is located away from the origin") {
+ Eigen::Vector3d origin = 100 * Eigen::Vector3d::Random();
+ Molecule point = dummy<0>(2.929075493, origin);
+ double area = 0.4;
+ double probeRadius = 0.0;
+ double minRadius = 100.0;
+ double minDistance = 0.1;
+ int derOrder = 4;
+ TsLessCavity cavity =
+ TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
+
+ CPCMSolver solver(symm, correction);
+ solver.buildSystemMatrix(cavity, gfInside, gfOutside, S);
+
+ size_t size = cavity.size();
+ Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge, origin);
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_mep(" << i << ") = " << fake_mep(i));
+ }
+
+ THEN("the surface charge is") {
+ Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
+ fake_asc = solver.computeCharge(fake_mep);
+ double totalFakeASC = fake_asc.sum();
+
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_asc(" << i << ") = " << fake_asc(i));
}
+
+ CAPTURE(totalASC);
+ CAPTURE(totalFakeASC);
+ CAPTURE(totalASC - totalFakeASC);
+ REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
+ }
}
+ }
}
tests/iefpcm/iefpcm_tsless-NH3.cpp--- tests/iefpcm/iefpcm_tsless-NH3.cpp
+++ tests/iefpcm/iefpcm_tsless-NH3.cpp
@@ -27,14 +27,14 @@
#include <Eigen/Core>
+#include "TestingMolecules.hpp"
#include "bi_operators/Collocation.hpp"
-#include "green/DerivativeTypes.hpp"
#include "cavity/TsLessCavity.hpp"
-#include "utils/Molecule.hpp"
-#include "green/Vacuum.hpp"
-#include "TestingMolecules.hpp"
+#include "green/DerivativeTypes.hpp"
#include "green/UniformDielectric.hpp"
+#include "green/Vacuum.hpp"
#include "solver/IEFSolver.hpp"
+#include "utils/Molecule.hpp"
#include "utils/Symmetry.hpp"
using namespace pcm;
@@ -47,36 +47,37 @@
/*! \class IEFSolver
* \test \b NH3TsLess tests IEFSolver using ammonia and a TsLess cavity
*/
-TEST_CASE("Test solver for the IEFPCM with NH3 molecule and a TsLess cavity", "[solver][iefpcm][iefpcm_tsless-NH3]")
-{
- Molecule molec = NH3();
-
- double area = 0.08;
- double minDistance = 0.1;
- double probeRadius = 0.0;
- int derOrder = 8;
- double minRadius = 100.0;
- TsLessCavity cavity = TsLessCavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
-
- double permittivity = 78.39;
- Vacuum<> gfInside;
- UniformDielectric<> gfOutside(permittivity);
-
- Collocation op;
-
- bool symm = true;
- IEFSolver solver(symm);
- solver.buildSystemMatrix(cavity, gfInside, gfOutside, op);
-
- double Ncharge = 7.0;
- double Hcharge = 1.0;
- size_t size = cavity.size();
- Eigen::VectorXd fake_mep = computeMEP(molec, cavity.elements());
- // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon]
- Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
- fake_asc = solver.computeCharge(fake_mep);
- double totalASC = - (Ncharge + 3.0 * Hcharge) * (permittivity - 1) / permittivity;
- double totalFakeASC = fake_asc.sum();
- CAPTURE(totalASC - totalFakeASC);
- REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
+TEST_CASE("Test solver for the IEFPCM with NH3 molecule and a TsLess cavity",
+ "[solver][iefpcm][iefpcm_tsless-NH3]") {
+ Molecule molec = NH3();
+
+ double area = 0.08;
+ double minDistance = 0.1;
+ double probeRadius = 0.0;
+ int derOrder = 8;
+ double minRadius = 100.0;
+ TsLessCavity cavity =
+ TsLessCavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
+
+ double permittivity = 78.39;
+ Vacuum<> gfInside;
+ UniformDielectric<> gfOutside(permittivity);
+
+ Collocation op;
+
+ bool symm = true;
+ IEFSolver solver(symm);
+ solver.buildSystemMatrix(cavity, gfInside, gfOutside, op);
+
+ double Ncharge = 7.0;
+ double Hcharge = 1.0;
+ size_t size = cavity.size();
+ Eigen::VectorXd fake_mep = computeMEP(molec, cavity.elements());
+ // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon]
+ Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
+ fake_asc = solver.computeCharge(fake_mep);
+ double totalASC = -(Ncharge + 3.0 * Hcharge) * (permittivity - 1) / permittivity;
+ double totalFakeASC = fake_asc.sum();
+ CAPTURE(totalASC - totalFakeASC);
+ REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
}
tests/iefpcm/iefpcm_tsless-point.cpp--- tests/iefpcm/iefpcm_tsless-point.cpp
+++ tests/iefpcm/iefpcm_tsless-point.cpp
@@ -27,13 +27,13 @@
#include <Eigen/Core>
+#include "TestingMolecules.hpp"
#include "bi_operators/Collocation.hpp"
-#include "green/DerivativeTypes.hpp"
#include "cavity/TsLessCavity.hpp"
-#include "green/Vacuum.hpp"
+#include "green/DerivativeTypes.hpp"
#include "green/UniformDielectric.hpp"
+#include "green/Vacuum.hpp"
#include "solver/IEFSolver.hpp"
-#include "TestingMolecules.hpp"
using namespace pcm;
using bi_operators::Collocation;
@@ -42,100 +42,99 @@
using green::UniformDielectric;
using solver::IEFSolver;
-SCENARIO("Test solver for the IEFPCM for a point charge and a TsLess cavity", "[solver][iefpcm][iefpcm_tsless-point]")
-{
- GIVEN("An isotropic environment and a point charge")
- {
- double permittivity = 78.39;
- Vacuum<> gfInside;
- UniformDielectric<> gfOutside(permittivity);
- bool symm = true;
-
- Collocation op;
-
- double charge = 8.0;
- double totalASC = - charge * (permittivity - 1) / permittivity;
-
- /*! \class IEFSolver
- * \test \b pointChargeTsLess tests IEFSolver using a point charge with a TsLess cavity
- * The point charge is at the origin.
- */
- WHEN("the point charge is located at the origin")
- {
- Molecule point = dummy<0>(2.929075493);
- double area = 0.4;
- double probeRadius = 0.0;
- double minRadius = 100.0;
- double minDistance = 0.1;
- int derOrder = 4;
- TsLessCavity cavity = TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
- cavity.saveCavity("tsless_point.npz");
-
- IEFSolver solver(symm);
- solver.buildSystemMatrix(cavity, gfInside, gfOutside, op);
-
- size_t size = cavity.size();
- Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge);
- // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon]
- Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
- fake_asc = solver.computeCharge(fake_mep);
-
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_mep(" << i << ") = " << fake_mep(i));
- }
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_asc(" << i << ") = " << fake_asc(i));
- }
-
- double totalFakeASC = fake_asc.sum();
- THEN("the apparent surface charge is")
- {
- CAPTURE(totalASC);
- CAPTURE(totalFakeASC);
- CAPTURE(totalASC - totalFakeASC);
- REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
- }
- }
-
- /*! \class IEFSolver
- * \test \b pointChargeShiftedTsLess tests IEFSolver using a point charge with a TsLess cavity
- * The point charge is away from the origin.
- */
- AND_WHEN("the point charge is located away from the origin")
- {
- Eigen::Vector3d origin = 100 * Eigen::Vector3d::Random();
- Molecule point = dummy<0>(2.929075493, origin);
- double area = 0.4;
- double probeRadius = 0.0;
- double minRadius = 100.0;
- double minDistance = 0.1;
- int derOrder = 4;
- TsLessCavity cavity = TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
-
- IEFSolver solver(symm);
- solver.buildSystemMatrix(cavity, gfInside, gfOutside, op);
-
- size_t size = cavity.size();
- Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge, origin);
- // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon]
- Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
- fake_asc = solver.computeCharge(fake_mep);
-
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_mep(" << i << ") = " << fake_mep(i));
- }
- for (size_t i = 0; i < size; ++i) {
- INFO("fake_asc(" << i << ") = " << fake_asc(i));
- }
-
- double totalFakeASC = fake_asc.sum();
- THEN("the surface charge is")
- {
- CAPTURE(totalASC);
- CAPTURE(totalFakeASC);
- CAPTURE(totalASC - totalFakeASC);
- REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
- }
- }
+SCENARIO("Test solver for the IEFPCM for a point charge and a TsLess cavity",
+ "[solver][iefpcm][iefpcm_tsless-point]") {
+ GIVEN("An isotropic environment and a point charge") {
+ double permittivity = 78.39;
+ Vacuum<> gfInside;
+ UniformDielectric<> gfOutside(permittivity);
+ bool symm = true;
+
+ Collocation op;
+
+ double charge = 8.0;
+ double totalASC = -charge * (permittivity - 1) / permittivity;
+
+ /*! \class IEFSolver
+ * \test \b pointChargeTsLess tests IEFSolver using a point charge with a TsLess
+ * cavity
+ * The point charge is at the origin.
+ */
+ WHEN("the point charge is located at the origin") {
+ Molecule point = dummy<0>(2.929075493);
+ double area = 0.4;
+ double probeRadius = 0.0;
+ double minRadius = 100.0;
+ double minDistance = 0.1;
+ int derOrder = 4;
+ TsLessCavity cavity =
+ TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
+ cavity.saveCavity("tsless_point.npz");
+
+ IEFSolver solver(symm);
+ solver.buildSystemMatrix(cavity, gfInside, gfOutside, op);
+
+ size_t size = cavity.size();
+ Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge);
+ // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon]
+ Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
+ fake_asc = solver.computeCharge(fake_mep);
+
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_mep(" << i << ") = " << fake_mep(i));
+ }
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_asc(" << i << ") = " << fake_asc(i));
+ }
+
+ double totalFakeASC = fake_asc.sum();
+ THEN("the apparent surface charge is") {
+ CAPTURE(totalASC);
+ CAPTURE(totalFakeASC);
+ CAPTURE(totalASC - totalFakeASC);
+ REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
+ }
+ }
+
+ /*! \class IEFSolver
+ * \test \b pointChargeShiftedTsLess tests IEFSolver using a point charge with a
+ * TsLess cavity
+ * The point charge is away from the origin.
+ */
+ AND_WHEN("the point charge is located away from the origin") {
+ Eigen::Vector3d origin = 100 * Eigen::Vector3d::Random();
+ Molecule point = dummy<0>(2.929075493, origin);
+ double area = 0.4;
+ double probeRadius = 0.0;
+ double minRadius = 100.0;
+ double minDistance = 0.1;
+ int derOrder = 4;
+ TsLessCavity cavity =
+ TsLessCavity(point, area, probeRadius, minRadius, minDistance, derOrder);
+
+ IEFSolver solver(symm);
+ solver.buildSystemMatrix(cavity, gfInside, gfOutside, op);
+
+ size_t size = cavity.size();
+ Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge, origin);
+ // The total ASC for a dielectric is -Q*[(epsilon-1)/epsilon]
+ Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
+ fake_asc = solver.computeCharge(fake_mep);
+
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_mep(" << i << ") = " << fake_mep(i));
+ }
+ for (size_t i = 0; i < size; ++i) {
+ INFO("fake_asc(" << i << ") = " << fake_asc(i));
+ }
+
+ double totalFakeASC = fake_asc.sum();
+ THEN("the surface charge is") {
+ CAPTURE(totalASC);
+ CAPTURE(totalFakeASC);
+ CAPTURE(totalASC - totalFakeASC);
+ REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
+ }
}
+ }
}
tests/tsless/tsless_HF.cpp--- tests/tsless/tsless_HF.cpp
+++ tests/tsless/tsless_HF.cpp
@@ -2,55 +2,52 @@
#include <Eigen/Core>
-#include "cavity/TsLessCavity.hpp"
#include "TestingMolecules.hpp"
+#include "cavity/TsLessCavity.hpp"
using namespace pcm;
using cavity::TsLessCavity;
-TEST_CASE("TsLess cavity for an hydrogen fluoride molecule", "[tsless][tsless_HF]")
-{
- Molecule molec = HF();
- double area = 0.02 / bohr2ToAngstrom2();
- double probeRadius = 0.0 / bohrToAngstrom();
- double minRadius = 0.2 / bohrToAngstrom();
- double minDistance = 0.01;
- int derOrder = 4;
- TsLessCavity cavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityHFTest_size tests TsLess cavity size for ammonia
- */
- SECTION("Test size")
- {
- size_t ref_size = 1498;
- size_t size = cavity.size();
- REQUIRE(size == ref_size);
- }
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityHFTest_area tests TsLess cavity surface area for ammonia
- */
- SECTION("Test surface area")
- {
- double ref_area = 111.07794350824591;
- double area = cavity.elementArea().sum();
- REQUIRE(area == Approx(ref_area));
- }
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityHFTest_volume tests TsLess cavity volume for ammonia
- */
- SECTION("Test volume")
- {
- double ref_volume = 106.59560252994531;
- Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
- Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
- double volume = 0;
- for ( size_t i = 0; i < cavity.size(); ++i ) {
- volume += cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
- }
- volume /= 3;
- REQUIRE(volume == Approx(ref_volume));
+TEST_CASE("TsLess cavity for an hydrogen fluoride molecule", "[tsless][tsless_HF]") {
+ Molecule molec = HF();
+ double area = 0.02 / bohr2ToAngstrom2();
+ double probeRadius = 0.0 / bohrToAngstrom();
+ double minRadius = 0.2 / bohrToAngstrom();
+ double minDistance = 0.01;
+ int derOrder = 4;
+ TsLessCavity cavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityHFTest_size tests TsLess cavity size for ammonia
+ */
+ SECTION("Test size") {
+ size_t ref_size = 1498;
+ size_t size = cavity.size();
+ REQUIRE(size == ref_size);
+ }
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityHFTest_area tests TsLess cavity surface area for ammonia
+ */
+ SECTION("Test surface area") {
+ double ref_area = 111.07794350824591;
+ double area = cavity.elementArea().sum();
+ REQUIRE(area == Approx(ref_area));
+ }
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityHFTest_volume tests TsLess cavity volume for ammonia
+ */
+ SECTION("Test volume") {
+ double ref_volume = 106.59560252994531;
+ Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
+ Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
+ double volume = 0;
+ for (size_t i = 0; i < cavity.size(); ++i) {
+ volume +=
+ cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
}
+ volume /= 3;
+ REQUIRE(volume == Approx(ref_volume));
+ }
}
tests/tsless/tsless_NH3.cpp--- tests/tsless/tsless_NH3.cpp
+++ tests/tsless/tsless_NH3.cpp
@@ -2,56 +2,53 @@
#include <Eigen/Core>
-#include "cavity/TsLessCavity.hpp"
#include "TestingMolecules.hpp"
+#include "cavity/TsLessCavity.hpp"
using namespace pcm;
using cavity::TsLessCavity;
-TEST_CASE("TsLess cavity for an ammonia molecule", "[tsless][tsless_NH3]")
-{
- Molecule molec = NH3();
- double area = 0.02 / bohr2ToAngstrom2();
- double probeRadius = 0.0 / bohrToAngstrom();
- double minRadius = 0.2 / bohrToAngstrom();
- double minDistance = 0.01;
- int derOrder = 4;
- TsLessCavity cavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
- cavity.saveCavity("nh3.npz");
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityNH3Test_size tests TsLess cavity size for ammonia
- */
- SECTION("Test size")
- {
- size_t ref_size = 2042;
- size_t size = cavity.size();
- REQUIRE(size == ref_size);
- }
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityNH3Test_area tests TsLess cavity surface area for ammonia
- */
- SECTION("Test surface area")
- {
- double ref_area = 147.59900736316496;
- double area = cavity.elementArea().sum();
- REQUIRE(area == Approx(ref_area));
- }
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityNH3Test_volume tests TsLess cavity volume for ammonia
- */
- SECTION("Test volume")
- {
- double ref_volume = 153.7079082474709;
- Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
- Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
- double volume = 0;
- for ( size_t i = 0; i < cavity.size(); ++i ) {
- volume += cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
- }
- volume /= 3;
- REQUIRE(volume == Approx(ref_volume));
+TEST_CASE("TsLess cavity for an ammonia molecule", "[tsless][tsless_NH3]") {
+ Molecule molec = NH3();
+ double area = 0.02 / bohr2ToAngstrom2();
+ double probeRadius = 0.0 / bohrToAngstrom();
+ double minRadius = 0.2 / bohrToAngstrom();
+ double minDistance = 0.01;
+ int derOrder = 4;
+ TsLessCavity cavity(molec, area, probeRadius, minRadius, minDistance, derOrder);
+ cavity.saveCavity("nh3.npz");
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityNH3Test_size tests TsLess cavity size for ammonia
+ */
+ SECTION("Test size") {
+ size_t ref_size = 2042;
+ size_t size = cavity.size();
+ REQUIRE(size == ref_size);
+ }
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityNH3Test_area tests TsLess cavity surface area for ammonia
+ */
+ SECTION("Test surface area") {
+ double ref_area = 147.59900736316496;
+ double area = cavity.elementArea().sum();
+ REQUIRE(area == Approx(ref_area));
+ }
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityNH3Test_volume tests TsLess cavity volume for ammonia
+ */
+ SECTION("Test volume") {
+ double ref_volume = 153.7079082474709;
+ Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
+ Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
+ double volume = 0;
+ for (size_t i = 0; i < cavity.size(); ++i) {
+ volume +=
+ cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
}
+ volume /= 3;
+ REQUIRE(volume == Approx(ref_volume));
+ }
}
tests/tsless/tsless_point.cpp--- tests/tsless/tsless_point.cpp
+++ tests/tsless/tsless_point.cpp
@@ -2,54 +2,52 @@
#include <Eigen/Core>
-#include "cavity/TsLessCavity.hpp"
#include "TestingMolecules.hpp"
+#include "cavity/TsLessCavity.hpp"
using namespace pcm;
using cavity::TsLessCavity;
-TEST_CASE("TsLess cavity for a single sphere", "[tsless][tsless_point]")
-{
- Molecule point = dummy<0>();
- double area = 0.4;
- double minDistance = 0.1;
- double probeRadius = 0.0;
- int derOrder = 4;
- TsLessCavity cavity(point, area, probeRadius, 100.0, minDistance, derOrder);
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityTest_size tests TsLess cavity size for a point charge
- */
- SECTION("Test size")
- {
- size_t ref_size = 32;
- size_t size = cavity.size();
- REQUIRE(ref_size == size);
- }
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityTest_area tests TsLess cavity surface area for a point charge
- */
- SECTION("Test surface area")
- {
- double ref_area = 4.0 * M_PI * pow(1.0, 2);
- double area = cavity.elementArea().sum();
- REQUIRE(ref_area == Approx(area));
- }
-
- /*! \class TsLessCavity
- * \test \b TsLessCavityTest_volume tests TsLess cavity volume for a point charge
- */
- SECTION("Test volume")
- {
- double ref_volume = 4.0 * M_PI * pow(1.0, 3) / 3.0;
- Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
- Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
- double volume = 0;
- for ( size_t i = 0; i < cavity.size(); ++i ) {
- volume += cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
- }
- volume /= 3;
- REQUIRE(ref_volume == Approx(volume));
+TEST_CASE("TsLess cavity for a single sphere", "[tsless][tsless_point]") {
+ Molecule point = dummy<0>();
+ double area = 0.4;
+ double minDistance = 0.1;
+ double probeRadius = 0.0;
+ int derOrder = 4;
+ TsLessCavity cavity(point, area, probeRadius, 100.0, minDistance, derOrder);
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityTest_size tests TsLess cavity size for a point charge
+ */
+ SECTION("Test size") {
+ size_t ref_size = 32;
+ size_t size = cavity.size();
+ REQUIRE(ref_size == size);
+ }
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityTest_area tests TsLess cavity surface area for a point
+ * charge
+ */
+ SECTION("Test surface area") {
+ double ref_area = 4.0 * M_PI * pow(1.0, 2);
+ double area = cavity.elementArea().sum();
+ REQUIRE(ref_area == Approx(area));
+ }
+
+ /*! \class TsLessCavity
+ * \test \b TsLessCavityTest_volume tests TsLess cavity volume for a point charge
+ */
+ SECTION("Test volume") {
+ double ref_volume = 4.0 * M_PI * pow(1.0, 3) / 3.0;
+ Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
+ Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
+ double volume = 0;
+ for (size_t i = 0; i < cavity.size(); ++i) {
+ volume +=
+ cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
}
+ volume /= 3;
+ REQUIRE(ref_volume == Approx(volume));
+ }
}
Generated by 🚫 Danger |
|
Superseded by #140 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Add the tessellationless integration grid for the PCM as described in Pomelli 2004
Description
The GePol algorithm is based on a partition of the spheres that involves the use of spherical polygons. The TsLess algorithm uses a different strategy for the partition of the spheres that makes away with the polygons altogether in favor of just points and weights.
How Has This Been Tested?
Todos
FCMangle.hppheader is used insideGePolCavityandTsLessCavityto handle name-mangling of Fortran backends.Types of changes
Questions
Status