diff --git a/CHANGELOG.md b/CHANGELOG.md index efbeae3..b4366b3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,13 +2,15 @@ ## [Unreleased] -## [Version 2.0.0-alpha1] - 2019-04-XY - We have introduced a number of breaking changes, motivated by the need to modernize the library. - The version 1 release series is available in the archived repository [PCMSolver/pcmsolver_v1] +### Added + +- A new cavity generator using C. S. Pomelli's [TsLess algorithm](http://dx.doi.org/10.1002/jcc.20076) + ### Changed - **BREAKING** We require a fully compliant C++14 compiler. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bd305a9..8d0a18e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -65,6 +65,7 @@ add_subdirectory(interface) add_subdirectory(pedra) add_subdirectory(solver) add_subdirectory(utils) +add_subdirectory(tsless) add_library(pcm $) diff --git a/src/cavity/CMakeLists.txt b/src/cavity/CMakeLists.txt index 28fe1ac..3cf989b 100644 --- a/src/cavity/CMakeLists.txt +++ b/src/cavity/CMakeLists.txt @@ -4,6 +4,7 @@ target_sources(pcm-object ${CMAKE_CURRENT_LIST_DIR}/GePolCavity.cpp ${CMAKE_CURRENT_LIST_DIR}/ICavity.cpp ${CMAKE_CURRENT_LIST_DIR}/RestartCavity.cpp + ${CMAKE_CURRENT_LIST_DIR}/TsLessCavity.cpp ) # List of headers @@ -14,6 +15,7 @@ list(APPEND headers_list GePolCavity.hpp ICavity.hpp RestartCavity.hpp + TsLessCavity.hpp ) # Sets install directory for all the headers in the list foreach(_header IN LISTS headers_list) diff --git a/src/cavity/Cavity.hpp b/src/cavity/Cavity.hpp index 9b0a668..3f24883 100644 --- a/src/cavity/Cavity.hpp +++ b/src/cavity/Cavity.hpp @@ -26,6 +26,7 @@ #include "GePolCavity.hpp" #include "ICavity.hpp" #include "RestartCavity.hpp" +#include "TsLessCavity.hpp" #include "utils/Factory.hpp" /*! @@ -49,6 +50,7 @@ inline Factory bootstrapFactory() { factory_.subscribe("GEPOL", createGePolCavity); factory_.subscribe("RESTART", createRestartCavity); + factory_.subscribe("TSLESS", createTsLessCavity); return factory_; } diff --git a/src/cavity/CavityData.hpp b/src/cavity/CavityData.hpp index ced375a..03d1883 100644 --- a/src/cavity/CavityData.hpp +++ b/src/cavity/CavityData.hpp @@ -40,6 +40,12 @@ struct CavityData { double area; /*! The radius of the spherical probe representing the solvent */ double probeRadius; + /*! The minimal distance between two sampling points on different spheres. Relevant + * for TsLessCavity */ + double minDistance; + /*! The maximum derivative order to be used in the definition of the smoothing + * function. Relevant for TsLessCavity */ + int derOrder; /*! Triggers the addition of spheres not centered on atoms, relevant for * GePolCavity */ double minimalRadius; @@ -50,12 +56,16 @@ struct CavityData { const Molecule & _molec, double _area, double _probeRadius, + double _minDistance, + int _derOrder, double _minRadius, const std::string & _fname) : cavityType(type), molecule(_molec), area(_area), probeRadius(_probeRadius), + minDistance(_minDistance), + derOrder(_derOrder), minimalRadius(_minRadius), filename(_fname) {} }; diff --git a/src/cavity/GePolCavity.cpp b/src/cavity/GePolCavity.cpp index bef6f22..1444f8e 100644 --- a/src/cavity/GePolCavity.cpp +++ b/src/cavity/GePolCavity.cpp @@ -81,7 +81,7 @@ void GePolCavity::build(const std::string & suffix, int maxsph, int maxvert) { - // This is a wrapper for the generatecavity_cpp function defined in the Fortran + // This is a wrapper for the pedra_driver function defined in the Fortran // code PEDRA. // Here we allocate the necessary arrays to be passed to PEDRA, in particular we // allow @@ -172,41 +172,41 @@ void GePolCavity::build(const std::string & suffix, pedra << "PEDRA.OUT_" << suffix << "_" << getpid(); int len_f_pedra = std::strlen(pedra.str().c_str()); // Go PEDRA, Go! - timer::timerON("GePolCavity::generatecavity_cpp"); - generatecavity_cpp(&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); - timer::timerOFF("GePolCavity::generatecavity_cpp"); + timer::timerON("GePolCavity::pedra_driver"); + 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); + timer::timerOFF("GePolCavity::pedra_driver"); // The "intensive" part of updating the spheres related class data members will be // of course diff --git a/src/cavity/GePolCavity.hpp b/src/cavity/GePolCavity.hpp index 9a40898..eb6430a 100644 --- a/src/cavity/GePolCavity.hpp +++ b/src/cavity/GePolCavity.hpp @@ -89,18 +89,7 @@ class GePolCavity final : public ICavity { void writeOFF(const std::string & suffix); }; -/*! \fn extern "C" void generatecavity_cpp(int * maxts, int * maxsph, int * maxvert, - * double * xtscor, double * ytscor, double * ztscor, - * double * ar, - * double * xsphcor, double * ysphcor, double * - * zsphcor, double * rsph, - * int * nts, int * ntsirr, int * nesfp, int * - * addsph, - * double * xe, double * ye, double * ze, double * - * rin, - * double * avgArea, double * rsolv, double * ret, - * int * nr_gen, int * gen1, int * gen2, int * gen3, - * int * nvert, double * vert, double * centr) +/*! \brief Interface to the Fortran PEDRA code * \param[in] maxts maximum number of tesserae allowed * \param[in] maxsph maximum number of spheres allowed * \param[in] maxvert maximum number of vertices allowed @@ -139,39 +128,39 @@ class GePolCavity final : public ICavity { * \param[out] vert coordinates of tesserae vertices * \param[out] centr centers of arcs defining the edges of the tesserae */ -extern "C" void generatecavity_cpp(int * maxts, - int * maxsph, - int * maxvert, - double * xtscor, - double * ytscor, - double * ztscor, - double * ar, - double * xsphcor, - double * ysphcor, - double * zsphcor, - double * rsph, - int * nts, - int * ntsirr, - int * nesfp, - int * addsph, - double * xe, - double * ye, - double * ze, - double * rin, - double * masses, - double * avgArea, - double * rsolv, - double * ret, - int * nr_gen, - int * gen1, - int * gen2, - int * gen3, - int * nvert, - double * vert, - double * centr, - int * isphe, - const char * pedra, - int * len_f_pedra); +extern "C" void pedra_driver(int * maxts, + int * maxsph, + int * maxvert, + double * xtscor, + double * ytscor, + double * ztscor, + double * ar, + double * xsphcor, + double * ysphcor, + double * zsphcor, + double * rsph, + int * nts, + int * ntsirr, + int * nesfp, + int * addsph, + double * xe, + double * ye, + double * ze, + double * rin, + double * masses, + double * avgArea, + double * rsolv, + double * ret, + int * nr_gen, + int * gen1, + int * gen2, + int * gen3, + int * nvert, + double * vert, + double * centr, + int * isphe, + const char * pedra, + int * len_f_pedra); ICavity * createGePolCavity(const CavityData & data); } // namespace cavity diff --git a/src/cavity/TsLessCavity.cpp b/src/cavity/TsLessCavity.cpp new file mode 100644 index 0000000..bb619ce --- /dev/null +++ b/src/cavity/TsLessCavity.cpp @@ -0,0 +1,341 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "TsLessCavity.hpp" + +#include +#include +#include +#include +#include + +#include + +#include "CavityData.hpp" +#include "utils/PhysicalConstants.hpp" +#include "utils/Sphere.hpp" +#include "utils/Symmetry.hpp" +#include "utils/Timer.hpp" + +namespace pcm { +namespace cavity { +TsLessCavity::TsLessCavity(const Molecule & molec, + double a, + double pr, + double minR, + double minD, + int der) + : ICavity(molec), + averageArea_(a), + probeRadius_(pr), + minimalRadius_(minR), + minDistance_(minD), + derOrder_(der) { + timer::timerON("TsLessCavity::build from Molecule"); + build(10000, 200, 25000); + timer::timerOFF("TsLessCavity::build from Molecule"); +} + +TsLessCavity::TsLessCavity(const Sphere & sph, + double a, + double pr, + double minR, + double minD, + int der) + : ICavity(sph), + averageArea_(a), + probeRadius_(pr), + minimalRadius_(minR), + minDistance_(minD), + derOrder_(der) { + timer::timerON("TsLessCavity::build from a single sphere"); + build(10000, 200, 25000); + timer::timerOFF("TsLessCavity::build from a single sphere"); +} + +TsLessCavity::TsLessCavity(const std::vector & sph, + double a, + double pr, + double minR, + double minD, + int der) + : ICavity(sph), + averageArea_(a), + probeRadius_(pr), + minimalRadius_(minR), + minDistance_(minD), + derOrder_(der) { + timer::timerON("TsLessCavity::build from list of spheres"); + build(10000, 200, 25000); + timer::timerOFF("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; + auto xtscor = new double[maxts]; + auto ytscor = new double[maxts]; + auto ztscor = new double[maxts]; + auto ar = new double[maxts]; + auto xsphcor = new double[maxts]; + auto ysphcor = new double[maxts]; + auto zsphcor = new double[maxts]; + auto rsph = new double[maxts]; + auto 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); + } + + auto xe = xv.data(); + auto ye = yv.data(); + auto ze = zv.data(); + + auto rin = radii_scratch.data(); + auto 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::timerON("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::timerOFF("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 containing the list of spheres + for (int i = orig; i < nSpheres_; ++i) { + spheres_.push_back(Sphere(sphereCenter_.col(i), sphereRadius_(i))); + } + } + + nElements_ = static_cast(nts); + nIrrElements_ = static_cast(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]; + } + + elementNormal_ = elementCenter_ - elementSphereCenter_; + for (int i = 0; i < nElements_; ++i) { + elementNormal_.col(i) /= elementNormal_.col(i).norm(); + } + + // 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)); + } + + 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_; + 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"; + os << "-------- ---- -------- ------- ----------- ----------- -----------\n"; + // Print original set of spheres + int original = nSpheres_ - addedSpheres; + Eigen::IOFormat CleanFmt(6, Eigen::DontAlignCols, " ", "\n", "", ""); + for (int i = 0; i < original; ++i) { + os << std::setw(4) << i + 1; + os << " " << molecule_.atoms()[i].symbol << " "; + os << std::fixed << std::setprecision(4) << molecule_.atoms()[i].radius; + os << std::fixed << std::setprecision(2) << " " + << molecule_.atoms()[i].radiusScaling << " "; + os << (molecule_.geometry().col(i).transpose() * bohrToAngstrom()) + .format(CleanFmt); + os << std::endl; + } + // Print added spheres + for (int j = 0; j < addedSpheres; ++j) { + int idx = original + j; + os << std::setw(4) << idx + 1; + os << " Du "; + os << std::fixed << std::setprecision(4) + << sphereRadius_(idx) * bohrToAngstrom(); + os << std::fixed << std::setprecision(2) << " 1.00"; + os << (sphereCenter_.col(idx).transpose() * bohrToAngstrom()).format(CleanFmt); + os << std::endl; + } + return os; +} + +ICavity * createTsLessCavity(const CavityData & data) { + return new TsLessCavity(data.molecule, + data.area, + data.probeRadius, + data.minimalRadius, + data.minDistance, + data.derOrder); +} +} // namespace cavity +} // namespace pcm diff --git a/src/cavity/TsLessCavity.hpp b/src/cavity/TsLessCavity.hpp new file mode 100644 index 0000000..8c12164 --- /dev/null +++ b/src/cavity/TsLessCavity.hpp @@ -0,0 +1,156 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#pragma once + +#include +#include +#include + +namespace pcm { +struct CavityData; +} // namespace pcm + +#include "ICavity.hpp" + +/*! \file TsLessCavity.hpp */ + +namespace pcm { +namespace cavity { +/*! \class TsLessCavity + * \brief A class for TsLess cavity. + * \author Roberto Di Remigio + * \date 2013 + * + * This class is a wrapper for the Fortran routines generating + * a tessellationless grid on the cavity surface. The original + * algoritm is described in \cite Pomelli2004 + */ +class TsLessCavity final : public ICavity { +public: + TsLessCavity() {} + TsLessCavity(const Molecule & molec, + double a, + double pr, + double minR, + double minD, + int der); + TsLessCavity(const Sphere & sph, + double a, + double pr, + double minR, + double minD, + int der); + TsLessCavity(const std::vector & sph, + double a, + double pr, + double minR, + double minD, + int der); + virtual ~TsLessCavity() {} + friend std::ostream & operator<<(std::ostream & os, TsLessCavity & cavity) { + return cavity.printCavity(os); + } + +private: + double averageArea_; + double probeRadius_; + double minimalRadius_; + double minDistance_; + int derOrder_; + int addedSpheres; + virtual std::ostream & printCavity(std::ostream & os) override; + virtual void makeCavity() override { build(10000, 200, 25000); } + /*! \brief Driver for TsLess Fortran module. + * \param[in] maxts maximum number of tesserae + * \param[in] maxsp maximum number of spheres (original + added) + * \param[in] maxvert maximum number of vertices + */ + void build(int maxts, int maxsp, int maxvert); +}; + +/*! \brief Fortran interface function to TsLess cavity generation + * \param[in] maxts maximum number of tesserae allowed + * \param[in] maxsph maximum number of spheres allowed + * \param[in] maxvert maximum number of vertices allowed + * \param[out] nesfp number of spheres (original + added) + * \param[out] nts number of generated tesserae + * \param[out] ntsirr number of generated irreducible tesserae + * \param[out] addsph number of added spheres + * \param[out] xtscor x-coordinate of tesserae centers (dimension maxts) + * \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] 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 \param[in] gen1 first generator \param[in] gen2 second generator + * \param[in] gen3 third generator + * \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] rsolv solvent probe radius \param[in] work scratch space + */ +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 +} // namespace pcm diff --git a/src/interface/Input.cpp b/src/interface/Input.cpp index 162e3b4..2f620ac 100644 --- a/src/interface/Input.cpp +++ b/src/interface/Input.cpp @@ -72,6 +72,8 @@ void Input::reader(const std::string & filename) { cavityType_ = cavity.getStr("TYPE"); area_ = cavity.getDbl("AREA"); + minDistance_ = cavity.getDbl("MINDISTANCE"); + derOrder_ = cavity.getInt("DERORDER"); if (cavityType_ == "RESTART") { cavFilename_ = cavity.getStr("NPZFILE"); } @@ -196,6 +198,8 @@ void Input::reader(const PCMInput & host_input) { cavityType_ = detail::trim_and_upper(host_input.cavity_type); area_ = host_input.area * angstrom2ToBohr2(); + minDistance_ = host_input.min_distance * angstromToBohr(); + derOrder_ = host_input.der_order; if (cavityType_ == "RESTART") { cavFilename_ = detail::trim(host_input.restart_name); // No case conversion here! } @@ -333,8 +337,14 @@ void Input::initMolecule() { } CavityData Input::cavityParams() const { - return CavityData( - cavityType_, molecule_, area_, probeRadius_, minimalRadius_, cavFilename_); + return CavityData(cavityType_, + molecule_, + area_, + probeRadius_, + minDistance_, + derOrder_, + minimalRadius_, + cavFilename_); } GreenData Input::insideGreenParams() const { diff --git a/src/interface/Input.hpp b/src/interface/Input.hpp index 9015daa..db2ef43 100644 --- a/src/interface/Input.hpp +++ b/src/interface/Input.hpp @@ -169,8 +169,12 @@ class PCMSolver_EXPORT Input { std::string cavityType_; /// Filename for the .npz cavity restart file std::string cavFilename_; - /// GePol cavity average element area + /// GePol and TsLess cavity average element area double area_; + /// TsLess cavity minimal distance between sampling points + double minDistance_; + /// TsLess cavity maximum derivative order of switch function + int derOrder_; /// Whether the radii should be scaled by 1.2 bool scaling_; /// The set of radii to be used diff --git a/src/pedra/pedra_cavity_interface.F90 b/src/pedra/pedra_cavity_interface.F90 index 6974750..6e46c13 100644 --- a/src/pedra/pedra_cavity_interface.F90 +++ b/src/pedra/pedra_cavity_interface.F90 @@ -27,14 +27,14 @@ ! ! RDR, 280114. Put things in makecav.F inside here directly. ! -subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & +subroutine pedra_driver(maxts_, maxsph_, maxvert_, & xtscor_, ytscor_, ztscor_, ar_, & xsphcor_, ysphcor_, zsphcor_, rsph_, & nts_, ntsirr_, nesfp_, addsph_, & xe_, ye_, ze_, rin_, masses_, avgArea_, rsolv_, ret_, & nr_gen_, gen1_, gen2_, gen3_, & nvert_, vert_, centr_, isphe_, pedra_, len_pedra_) & - bind(C, name='generatecavity_cpp') + bind(C, name='pedra_driver') use, intrinsic :: iso_c_binding use pedra_precision @@ -70,7 +70,6 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & character(len=len_pedra_) :: pedra type(point_group) :: pgroup - pedra = carray_to_fstring(pedra_) !lvpri = 121201_regint_k ! The following INQUIRE statement returns whether the file named cavity.off is @@ -172,4 +171,4 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & close(pedra_unit) -end subroutine generatecavity_cpp +end subroutine diff --git a/src/tools/pcmparser.py b/src/tools/pcmparser.py index 4e2fec1..90a713b 100644 --- a/src/tools/pcmparser.py +++ b/src/tools/pcmparser.py @@ -170,47 +170,57 @@ def setup_keywords(): # Cavity section cavity = Section("CAVITY", callback=verify_cavity) # Type of the cavity - # Valid values: GEPOL and RESTART + # Valid values: GEPOL, TSLESS and RESTART cavity.add_kw("TYPE", "STR") # Name of the file containing data for restarting a cavity # Valid for: Restart cavity # Default: empty string cavity.add_kw("NPZFILE", "STR", "") # Average area (in au^2) - # Valid for: GePol + # Valid for: GePol and TsLess cavities # Valid values: double strictly greater than 0.01 au^2 # Default: 0.3 au^2 cavity.add_kw("AREA", "DBL", 0.3) + # Minimal distance between sampling points (in au) + # Valid for: TsLess cavity + # Valid values: double + # Default: 0.1 au + cavity.add_kw("MINDISTANCE", "DBL", 0.1) + # Derivative order for the weight function + # Valid for: TsLess cavity + # Valid values: integer + # Default: 4 + cavity.add_kw("DERORDER", "INT", 4) # Scaling of the atomic radii - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: boolean # Default: True cavity.add_kw("SCALING", "BOOL", True) # Built-in radii set - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: BONDI, UFF or ALLINGER # Default: BONDI cavity.add_kw("RADIISET", "STR", "BONDI") # Minimal radius of added spheres (in au) - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: double greater than 0.4 au # Default: 100.0 au (no added spheres) cavity.add_kw("MINRADIUS", "DBL", 100.0) # Spheres geometry creation mode - # Valid for: GePol + # Valid for: GePol and TsLess # Valid values: EXPLICIT, ATOMS or IMPLICIT # Default: IMPLICIT cavity.add_kw("MODE", "STR", "IMPLICIT") # List of atoms with custom radius - # Valid for: GePol in MODE=ATOMS + # Valid for: GePol and TsLess in MODE=ATOMS # Valid values: array of integers cavity.add_kw("ATOMS", "INT_ARRAY") # List of custom radii - # Valid for: GePol in MODE=ATOMS + # Valid for: GePol and TsLess in MODE=ATOMS # Valid values: array of doubles cavity.add_kw("RADII", "DBL_ARRAY") # List of spheres - # Valid for: GePol in MODE=EXPLICIT + # Valid for: GePol and TsLess in MODE=EXPLICIT # Valid values: array of doubles in format [x, y, z, R] cavity.add_kw("SPHERES", "DBL_ARRAY", callback=verify_spheres) top.add_sect(cavity) @@ -378,7 +388,7 @@ def verify_top(section): def verify_cavity(section): - allowed = ("GEPOL", "RESTART") + allowed = ("GEPOL", "TSLESS", "RESTART") type = section.get("TYPE") if type.get() not in allowed: print( @@ -396,6 +406,8 @@ def verify_cavity(section): convert_area_scalar(section["AREA"]) if section["MINRADIUS"].is_set(): convert_length_scalar(section["MINRADIUS"]) + if section["MINDISTANCE"].is_set(): + convert_length_scalar(section["MINDISTANCE"]) if type.get() == "GEPOL": area = section.get("AREA") @@ -420,6 +432,18 @@ def verify_cavity(section): ) ) sys.exit(1) + elif type.get() == "TSLESS": + area = section.get("AREA") + a = area.get() + if a < 0.01: + print( + ( + "Requested area value too small: {}. Minimal value is: 0.01 au^2".format( + a + ) + ) + ) + sys.exit(1) elif type.get() == "RESTART": npzfile = section.get("NPZFILE") # Restart string is the filename, with extension, where the cavity specs were saved. diff --git a/src/tsless/CMakeLists.txt b/src/tsless/CMakeLists.txt new file mode 100644 index 0000000..984ebad --- /dev/null +++ b/src/tsless/CMakeLists.txt @@ -0,0 +1,9 @@ +target_sources(pcm-object + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_cavity.f90 + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_lapack.f90 + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_precision.f90 + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_symmetry.f90 + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_weight_function.f90 + ) + diff --git a/src/tsless/tsless_cavity.f90 b/src/tsless/tsless_cavity.f90 new file mode 100644 index 0000000..f773ca5 --- /dev/null +++ b/src/tsless/tsless_cavity.f90 @@ -0,0 +1,390 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can 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 3 of the License, or +! (at your option) any later version. +! +! PCMSolver 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +!> Originally written by Christian Silvio Pomelli (ca. 2003-2004) +!> Introduced in PCMSolver by Christian Silvio Pomelli (ca. 2014) +!> Wrapped in a F90 module by Roberto Di Remigio (2015) +module tsless_cavity + +use, intrinsic :: iso_c_binding +use tsless_precision +use tsless_symmetry +use tsless_weight_function + +implicit none + +public tsless_driver + +private + +contains + + !> \brief Interface function to TsLess cavity generation + !> \author Roberto Di Remigio, Christian S. Pomelli + !> \year 2014, 2015 + !> \param[in] maxts maximum number of tesserae allowed + !> \param[in] maxsph maximum number of spheres allowed + !> \param[in] maxvert maximum number of vertices allowed + !> \param[out] nesfp number of spheres (original + added) + !> \param[out] nts number of generated tesserae + !> \param[out] ntsirr number of generated irreducible tesserae + !> \param[out] addsph number of added spheres + !> \param[out] xtscor x-coordinate of tesserae centers + !> \param[out] ytscor y-coordinate of tesserae centers + !> \param[out] ztscor z-coordinate of tesserae centers + !> \param[out] ar area of the tessera + !> \param[out] xsphcor x-coordinate of the sphere center the tessera belongs to + !> \param[out] ysphcor y-coordinate of the sphere center the tessera belongs to + !> \param[out] zsphcor z-coordinate of the sphere center the tessera belongs to + !> \param[out] rsph radii of the sphere the tessera belongs to, i.e. its curvature + !> \param[out] xe x-coordinate of the sphere center + !> \param[out] ye y-coordinate of the sphere center + !> \param[out] ze z-coordinate of the sphere center + !> \param[out] rin radius of the spheres + !> \param[in] masses atomic masses (for inertia tensor formation in TSLESS) + !> \param[in] nr_gen number of symmetry generators + !> \param[in] gen1 first generator + !> \param[in] gen2 second generator + !> \param[in] gen3 third generator + !> \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] rsolv solvent probe radius + !> \param[in] work scratch space + subroutine tsless_driver(maxts, maxsph, maxvert, nesfp, nts, ntsirr, addsph, & + xtscor, ytscor, ztscor, ar, & + xsphcor, ysphcor, zsphcor, rsph, & + xe, ye, ze, rin, & + masses, & + nr_gen, gen1, gen2, gen3, & + avgArea, dmin, nord, ifun, rsolv, work) & + bind(C, name='tsless_driver') + + !> Passed variables + integer(c_size_t), intent(in) :: maxts, maxsph, maxvert + integer(c_int), intent(in) :: nesfp + real(c_double), intent(out) :: xtscor(maxts), ytscor(maxts), ztscor(maxts), ar(maxts) + real(c_double), intent(out) :: xsphcor(maxts), ysphcor(maxts), zsphcor(maxts), rsph(maxts) + real(c_double), intent(in) :: xe(nesfp), ye(nesfp), ze(nesfp), rin(nesfp) + real(c_double), intent(in) :: masses(nesfp) + real(c_double), intent(in) :: work(maxts*maxsph) + real(c_double), intent(in) :: avgArea, rsolv, dmin + integer(c_int), intent(in) :: nord, ifun + integer(c_int), intent(out) :: nts, ntsirr, addsph + integer(c_int), intent(in) :: nr_gen, gen1, gen2, gen3 + !> Parameters + real(kind=dp), parameter :: pi = acos(-1.0_dp) + real(kind=dp), parameter :: tpi = 2.0_dp * pi + real(kind=dp), parameter :: fpi = 4.0_dp * pi + !> Local variables + integer :: print_unit + integer :: i, j, k, ntssp + logical :: tsless_file_exists, tsless_file_opened + integer(kind=regint_k) :: nesf_orig + type(weight_function) :: weights + real(kind=dp) :: x, y, z, r, dscale(3) + real(kind=dp) :: dist, dist2, rho, t, w0 + type(point_group) :: pgroup + + print_unit = 121221_regint_k + tsless_file_exists = .false. + inquire(file = 'TSLESS.OUT', exist = tsless_file_exists, opened = tsless_file_opened) + if (.not. tsless_file_opened) then + open(print_unit, & + file = 'TSLESS.OUT', & + status = 'replace', & + form = 'formatted', & + access = 'sequential') + end if + + !> Save original number of spheres + nesf_orig = nesfp + !> Create point group + pgroup = build_point_group(nr_gen, gen1, gen2, gen3, print_unit) + write(print_unit, *) 'Passed from host: nord, dmin, avgarea', nord, dmin, avgarea + + ! version 0.1 10/1/2014: plain implementation + ! added sphere positions, radii and derivatives (TBI) + ! generation and weighting of sampling points + + ! Generate weight function + weights = create_weight_function(nderiv=nord, normalization=ifun, dmin=dmin, printer=print_unit) + nts = 0_regint_k + + ! main loop on spheres (in future will implement linear scaling) + LoopOnSpheres: do i = 1, nesfp + ! point generation with weighting according to Leopardi scheme [ref] + ntssp = int(fpi * rin(i)**2 / avgArea) + 1_regint_k + !> Original weight factor + w0 = fpi * rin(i)**2 / ntssp + !> Interaction radius + rho = sqrt(w0 * pi) + write(print_unit, '(a, f12.5, a, f12.5)') 'Original weight factor ', w0, ' and interaction radius ', rho + x = xe(i) + y = ye(i) + z = ze(i) + r = rin(i) + + ! generate sphere sampling points + write(print_unit, '(a, i4, a, i4, a)') 'Generating ', ntssp, ' sampling points on sphere', i, ' according to Leopardi' + write(print_unit, '(a, f12.5, a, f12.5, a, f12.5, a, f12.5)') 'Center = (', x, ', ', y, ', ', z, ') Radius = ', r + call leopardi_grid(ntssp, nts, x, y, z, r, xtscor, ytscor, ztscor, print_unit) + write(print_unit, '(a)') 'Sampling points generated, now refining for intersections between spheres' + + ! loop on sphere sampling points + LoopOnSamplingPoints: do j = nts + 1, nts + ntssp + ar(j) = w0 + ! loop on intersecting spheres + LoopOnIntersectingSpheres: do k = 1, nesfp + if (k .eq. i) then + goto 10 + endif + dist2 = (xe(k) - xtscor(j))**2 + (ye(k) - ytscor(j))**2 + (ze(k) - ztscor(j))**2 + dist = sqrt(dist2) + t = (dist - rin(k)) / rho + write(print_unit, *) 't = ', t + if (t .le. dmin) then + ar(j) = -1.0_dp + write(print_unit, *) 'point ', j, ' was deleted' + goto 20 + else if(t .le. 1.0_dp) then + write(print_unit, *) 'point ', j, ' was scaled' + write(print_unit, *) 'Original value ', ar(j) + dscale = scaling_factor(weights, t, dmin) + ar(j) = ar(j) * dscale(1) + write(print_unit, *) 'Current value ', ar(j), ' Scaling factor ', dscale(1) + endif + 10 continue + end do LoopOnIntersectingSpheres + 20 continue + + ! save points with non zero weights + SavePoints: if (ar(j) .gt. 0.0_dp) then + nts = nts + 1_regint_k + xtscor(nts) = xtscor(j) + ytscor(nts) = ytscor(j) + ztscor(nts) = ztscor(j) + ar(nts) = ar(j) + xsphcor(nts) = x + ysphcor(nts) = y + zsphcor(nts) = z + rsph(nts) = r + end if SavePoints + end do LoopOnSamplingPoints + end do LoopOnSpheres + + !> Number of irreducible sampling points, hardcoded to nts for the moment + ntsirr = nts + !> Number of added spheres + addsph = nesfp - nesf_orig + + call destroy_weight_function(weights) + + write(print_unit, *) ' Surface = ', cavity_surface(nts, ar), ' volume = ', cavity_volume(nts, rsph, ar) + close(print_unit) + + call write_visualization_file(nts, nesfp, xtscor, ytscor, ztscor, ar, xe, ye, ze, rin) + + end subroutine tsless_driver + + !> \brief Create and translate Leopardi points + !> \author Christian S. Pomelli + !> \param[in] ntssp number of points on the sphere + !> \param[in] nts a counter for the total number of points generated + !> \param[in] x x-coordinate of sphere center + !> \param[in] y y-coordinate of sphere center + !> \param[in] z z-coordinate of sphere center + !> \param[in] r radius of sphere + !> \param[in] xtscor x-coordinates of points on the sphere + !> \param[in] ytscor y-coordinates of points on the sphere + !> \param[in] ztscor z-coordinates of points on the sphere + !> \param[in] print_unit + subroutine leopardi_grid(ntssp, nts, x, y, z, r, xtscor, ytscor, ztscor, print_unit) + + use, intrinsic :: iso_fortran_env, only: output_unit + + !> Passed variables + integer(kind=regint_k), intent(in) :: ntssp, nts + integer, optional, intent(in) :: print_unit + real(kind=dp), intent(in) :: x, y, z, r + real(kind=dp), intent(inout) :: xtscor(*), ytscor(*), ztscor(*) + !> Parameters + real(kind=dp), parameter :: pi = acos(-1.0_dp) + real(kind=dp), parameter :: tpi = 2.0_dp * pi + real(kind=dp), parameter :: fpi = 4.0_dp * pi + !> Local variables + integer :: i, j, k, n, print_out + real(kind=dp) :: thetaf(ntssp), ypsi(ntssp), theta(ntssp), apsi(ntssp), m(ntssp) + real(kind=dp) :: vr, thetac, delta1, deltaf + real(kind=dp) :: thetam, costm, sintm, omegaj + + if (present(print_unit)) then + print_out = print_unit + else + print_out = output_unit + end if + + thetaf = 0.0_dp + ypsi = 0.0_dp + theta = 0.0_dp + apsi = 0.0_dp + m = 0.0_dp + ! Notation consistent with [ref] + vr = fpi / ntssp + thetac = 2.0_dp * (asin(sqrt(vr / fpi))) + delta1 = sqrt(vr) + + ! Polar points + ! North pole + xtscor(nts+1) = x + ytscor(nts+1) = y + ztscor(nts+1) = z + r + ! South pole + xtscor(nts+ntssp) = x + ytscor(nts+ntssp) = y + ztscor(nts+ntssp) = z - r + + ! two points case (rare) + if (ntssp == 2) return + + ! Number of collars + n = int((pi - 2.0_dp * thetac) / delta1 + 0.5_dp) + if (n == 0) n = 1 + deltaf = (pi - 2.0_dp * thetac) / (1.0_dp * n) + + ! Colatitudes of collars + CollarColatitudes: do i = 1, n + 1 + thetaF(i) = thetac + (i - 1) * deltaf + end do CollarColatitudes + + ! Collars building + ypsi(1) = fpi*((sin(thetaf(2) / 2.0_dp))**2 - & + (sin(Thetaf(1) / 2.0_dp))**2) / vr + m(1) = int(ypsi(1) + 0.5_dp) + apsi(1) = ypsi(1) - m(1) + Collars: do i = 2, n + ypsi(i) = fpi * ((sin(thetaf(i+1) / 2.0_dp))**2 - & + (sin(thetaf(i ) / 2.0_dp))**2) / vr + m(i) = int(ypsi(i) + apsi(i-1) + 0.5_dp) + apsi(i) = ypsi(i) - m(i) + apsi(i-1) + end do Collars + + j = 1 + do i = 1, n + 1 + theta(i) = 2.0_dp * asin(sqrt((vr * j) / fpi)) + j = j + int(m(i) + 0.5_dp) + end do + + k = 2 + do i = 1, n + thetam = (theta(i) + theta(i+1)) / 2.0_dp + costm = cos(thetam) + sintm = sin(thetam) + TranslatePoints: do j = 1, int(m(i) + 0.5_dp) + omegaj = tpi * j / (1.0_dp * m(i)) + xtscor(nts+k) = x + r * sintm * cos(omegaj) + ytscor(nts+k) = y + r * sintm * sin(omegaj) + ztscor(nts+k) = z + r * costm + k=k+1 + end do TranslatePoints + end do + + end subroutine leopardi_grid + + subroutine write_visualization_file(nts, nesfp, xtscor, ytscor, ztscor, ar, xe, ye, ze, rin) + + !> writes a script for VMD + + !> Passed variables + integer(kind=regint_k), intent(in) :: nts, nesfp + real(kind=dp), intent(in) :: xtscor(nts), ytscor(nts), ztscor(nts), ar(nts), xe(nesfp), ye(nesfp), ze(nesfp), rin(nesfp) + !> Local variables + integer(kind=regint_k) :: file_unit + logical :: file_exists, file_opened + integer :: k + + file_unit = 221221_regint_k + file_exists = .false. + inquire(file = 'tsless_visual.vmd', exist = file_exists, opened = file_opened) + if (.not. file_opened) then + open(file_unit, & + file = 'tsless_visual.vmd', & + status = 'replace', & + form = 'formatted', & + access = 'sequential') + end if + + write(file_unit, *) 'mol load graphics name' + + do k = 1, nesfp + write(file_unit,1000) xe(k), ye(k), ze(k), rin(k) + enddo + + write(file_unit, *) 'graphics top color 3' + + do k = 1, nts + write(file_unit,1000) xtscor(k), ytscor(k), ztscor(k), 0.1_dp + end do + write(file_unit, *) + + close(file_unit) + +1000 format ('graphics top sphere {', 3F10.4, '} radius', F10.4, ' resolution 120') + + end subroutine write_visualization_file + + pure function cavity_surface(nts, ar) result(surface) + + !> Passed variables + integer(kind=regint_k), intent(in) :: nts + real(kind=dp), intent(in) :: ar(nts) + real(kind=dp) :: surface + !> Local variables + integer :: i + + surface = 0.0_dp + do i = 1, nts + surface = surface + ar(i) + end do + + end function cavity_surface + + pure function cavity_volume(nts, rsph, ar) result(volume) + + !> Passed variables + integer(kind=regint_k), intent(in) :: nts + real(kind=dp), intent(in) :: rsph(nts), ar(nts) + real(kind=dp) :: volume + !> Local variables + integer :: i + + volume = 0.0_dp + do i = 1, nts + volume = volume + ar(i) * rsph(i) + end do + volume = volume / 3.0_dp + + end function cavity_volume + +end module tsless_cavity diff --git a/src/tsless/tsless_lapack.f90 b/src/tsless/tsless_lapack.f90 new file mode 100644 index 0000000..aed4731 --- /dev/null +++ b/src/tsless/tsless_lapack.f90 @@ -0,0 +1,158 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can 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 3 of the License, or +! (at your option) any later version. +! +! PCMSolver 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_lapack + +use tsless_precision + +implicit none + +public solve_linear_system + +private + +contains + + !> \brief Solve linear system of equations + !> \author Roberto Di Remigio + !> \year 2015 + !> \param[in] ndim dimension of the problem + !> \param[in] A system matrix + !> \param[in] b right-hand side + !> \param[out] x solution vector + !> Computes the explicit inverse by LU decomposition and then + !> performs matrix-vector multiplication A^-1 b = x + pure subroutine solve_linear_system(ndim, A, b, x) + + !> Passed variables + integer(kind=regint_k), intent(in) :: ndim + real(kind=dp), intent(inout) :: A(ndim, ndim) + real(kind=dp), intent(in) :: b(ndim) + real(kind=dp), intent(out) :: x(ndim) + !> Local variables + real(kind=dp), allocatable :: Ainverse(:, :) + integer(kind=regint_k) :: i, j + + !> Compute inverse by LU decomposition + allocate(Ainverse(ndim, ndim)); Ainverse = 0.0_dp + call inverse_LU(ndim, A, Ainverse) + + !> Solve linear system + x = 0.0_dp + do i = 1_regint_k, ndim + do j = 1_regint_k, ndim + x(i) = x(i) + Ainverse(i, j) * b(j) + end do + end do + + deallocate(Ainverse) + + end subroutine solve_linear_system + + !> \brief Computes inverse matrix + !> \author Alexander L. Godunov + !> \year 2009 + !> \param[in, out] A matrix to be inverted + !> \param[in, out] C inverse matrix + !> \param[in] ndim dimension of the matrix + !> The method is based on Doolittle LU factorization for Ax = b + !> notice that the contents of the original matrix will be destroyed + pure subroutine inverse_LU(ndim, A, C) + + !> Passed variables + integer(kind=regint_k), intent(in) :: ndim + real(kind=dp), intent(inout) :: A(ndim, ndim) + real(kind=dp), intent(inout) :: C(ndim, ndim) + !> Local variables + real(kind=dp), allocatable :: L(:, :), U(:, :) + real(kind=dp), allocatable :: b(:), d(:), x(:) + real(kind=dp) :: coeff + integer(kind=regint_k) :: i, j, k + + ! step 0: initialization for matrices L and U and b + allocate(L(ndim, ndim)); L = 0.0_dp + allocate(U(ndim, ndim)); U = 0.0_dp + allocate(b(ndim)); b = 0.0_dp + allocate(d(ndim)); d = 0.0_dp + allocate(x(ndim)); x = 0.0_dp + + ! step 1: forward elimination + do k = 1, ndim-1 + do i = k + 1, ndim + coeff = A(i, k) / A(k, k) + L(i, k) = coeff + do j = k + 1, ndim + A(i, j) = A(i, j) - coeff * A(k, j) + end do + end do + end do + + ! Step 2: prepare L and U matrices + ! L matrix is a matrix of the elimination coefficient + ! + the diagonal elements are 1.0 + do i = 1, ndim + L(i, i) = 1.0_dp + end do + ! U matrix is the upper triangular part of A + do j = 1, ndim + do i = 1, j + U(i, j) = A(i, j) + end do + end do + + ! Step 3: compute columns of the inverse matrix C + do k = 1, ndim + b(k) = 1.0_dp + d(1) = b(1) + ! Step 3a: Solve Ld=b using the forward substitution + do i = 2, ndim + d(i) = b(i) + do j = 1, i - 1 + d(i) = d(i) - L(i, j) * d(j) + end do + end do + ! Step 3b: Solve Ux=d using the back substitution + x(ndim) = d(ndim) / U(ndim, ndim) + do i = ndim - 1, 1, -1 + x(i) = d(i) + do j = ndim, i + 1, -1 + x(i) = x(i) - U(i, j) * x(j) + end do + x(i) = x(i) / U(i, i) + end do + ! Step 3c: fill the solutions x(ndim) into column k of C + do i = 1, ndim + C(i, k) = x(i) + end do + b(k) = 0.0_dp + end do + + deallocate(L) + deallocate(U) + deallocate(b) + deallocate(d) + deallocate(x) + + end subroutine inverse_LU + +end module tsless_lapack diff --git a/src/tsless/tsless_precision.f90 b/src/tsless/tsless_precision.f90 new file mode 100644 index 0000000..627605a --- /dev/null +++ b/src/tsless/tsless_precision.f90 @@ -0,0 +1,42 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can 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 3 of the License, or +! (at your option) any later version. +! +! PCMSolver 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_precision +! Read this: http://stackoverflow.com/a/3170438/2528668 +! and this: http://stackoverflow.com/a/3204981/2528668 + +implicit none + +! Integer types +! 32-bit integers +integer, parameter :: regint_k = selected_int_kind(8) +! 64-bit integers +integer, parameter :: largeint_k = selected_int_kind(18) + +! Real types +! Single-precision real +integer, parameter :: sp = kind(1.0) +! Double-precision real +integer, parameter :: dp = selected_real_kind(2*precision(1.0_sp)) + +end module tsless_precision diff --git a/src/tsless/tsless_symmetry.f90 b/src/tsless/tsless_symmetry.f90 new file mode 100644 index 0000000..92e1c2c --- /dev/null +++ b/src/tsless/tsless_symmetry.f90 @@ -0,0 +1,354 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can 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 3 of the License, or +! (at your option) any later version. +! +! PCMSolver 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_symmetry +! NOTE: the ibtfun module is gone, as we can safely use Fortran +! standard intrinsic functions. +! Mapping between intrinsics and ibtfun: +! ibtand(i, j) <-> iand(i, j) +! ibtor(i, j) <-> ior(i, j) +! ibtshl(i, j) <-> ishft(i, j) +! ibtshr(i, j) <-> ishft(i, -j) !WARNING! +! ibtxor(i, j) <-> ieor(i, j) + +use, intrinsic :: iso_c_binding +use tsless_precision + +implicit none + +public get_pt +public build_point_group +! MINI-MANUAL +! A. Indexing of symmetry operations and their mapping to a bitstring: +! zyx Parity +! 0 000 E 1.0 +! 1 001 Oyz -1.0 +! 2 010 Oxz -1.0 +! 3 011 C2z 1.0 +! 4 100 Oxy -1.0 +! 5 101 C2y 1.0 +! 6 110 C2x 1.0 +! 7 111 i -1.0 +! B. Indexing of isymax matrix +! The isymax array contains the irrep to which the linear functions +! (first column) and the rotations (second column) belong. The indexing +! is given above. +! C. Indexing of the jsop array +! The jsop array contains the position at which the operation +! i (index given in point A) appears +! + +type, public :: point_group +! A type containing all you need to know about the point group + + ! String with the group name: + ! C1, C2, Cs, Ci, D2, C2v, C2h, D2h + character(len=3) :: group_name + ! Integer representing the group + ! 0, 1, 2, 3, 4, 5, 6, 7 + integer(kind=regint_k) :: group_int + ! Number of generators + integer(kind=regint_k) :: nr_generators + ! Number of not-trivial symmetry operations (2**nr_generators - 1) + integer(kind=regint_k) :: maxrep + ! group%isymax(i, 1): behaviour of principal axes under basic operations + ! (x-y-z) + ! group%isymax(i, 2): behaviour of principal rotations under basic operations + ! (Rx-Ry-Rz) + integer(kind=regint_k) :: isymax(3, 2) + ! Symmetry operations in the Abelian groups. + ! Bitstring: 1 coordinate changes sign under operation; + ! 0 coordinate does not change sign. + ! Of course, that's also the binary representation of + ! numbers from 0 to 7! + integer(kind=regint_k) :: jsop(0:7) + ! + integer(kind=regint_k) :: nr_rotations + ! + integer(kind=regint_k) :: nr_reflections + ! + integer(kind=regint_k) :: nr_inversion +end type point_group + +private +contains + + !> \brief returns parity of a bitstring + !> \param[in] bit_rep bitstring represenation of symmetry operator + !> PT is the parity of a bitstring: + !> 1 for an even number of ones: 000,011,110,101 + !> -1 for an odd number of ones: 001,010,100,111 + real(kind=dp) function get_pt(bit_rep) + + integer(kind=regint_k), intent(in) :: bit_rep + + real(kind=dp) :: pt(0:7) + + pt(0) = 1.0d0 + pt(1) = -1.0d0 + pt(2) = -1.0d0 + pt(3) = 1.0d0 + pt(4) = -1.0d0 + pt(5) = 1.0d0 + pt(6) = 1.0d0 + pt(7) = -1.0d0 + + get_pt = pt(bit_rep) + + end function get_pt + + !> \brief Builds point group given the generators + !> \param[in] nr_gen number of generators + !> \param[in] gen1 first generator + !> \param[in] gen2 second generator + !> \param[in] gen3 third generator + !> \param[in] print_unit logical unit for printing + !> Originally written by Trond Saue for DALTON/DIRAC + !> Copied and adapted by Roberto Di Remigio + !> The generators are encoded as bitstrings + function build_point_group(nr_gen, gen1, gen2, gen3, printer) result(pg) + + use, intrinsic :: iso_fortran_env, only: output_unit + + !> Passed variables + integer(kind=regint_k), intent(in) :: nr_gen + integer(kind=regint_k), intent(in) :: gen1, gen2, gen3 + integer, optional, intent(in) :: printer + !> Output variables + type(point_group) :: pg + !> Local variables + integer :: print_out + integer(kind=regint_k) :: maxrep + integer(kind=regint_k) :: isymax(3, 2) + integer(kind=regint_k) :: igen(3) + !> Integer representation of the rotations bitmaps + integer(kind=regint_k), parameter :: irots(3) = [3, 5, 6] + integer(kind=regint_k), parameter :: rots(3) = [6, 5, 3] + !> Integer representation of the reflections bitmaps + integer(kind=regint_k), parameter :: irefl(3) = [4, 2, 1] + !> Parity of the symmetry operations bitmaps + integer(kind=regint_k), parameter :: jpar(0:7) = [1, -1, -1, 1, -1, 1, 1, -1] + integer(kind=regint_k) :: i, j, k, l, i0, i1, i2, ind, ipos, bitmap + integer(kind=regint_k) :: nrots, nrefl, ninvc, igroup + integer(kind=regint_k) :: char_tab(0:7, 0:7) + logical :: lsymop(0:7) + integer(kind=regint_k) :: jsop(0:7), ipar(0:7) + character(3) :: group, rep(0:7) + character(3), parameter :: groups(0:7) = ['C1 ', 'C2 ', 'Cs ', & + 'Ci ', 'D2 ', 'C2v', & + 'C2h', 'D2h'] + character(3), parameter :: symop(0:7) = [' E ', 'Oyz', 'Oxz', & + 'C2z', 'Oxy', 'C2y', & + 'C2x', ' i '] + + if (present(printer)) then + print_out = printer + else + print_out = output_unit + end if + + isymax = 0 + igen = 0 + maxrep = 2**nr_gen - 1 + ! igen contains the bitmap for the generators + igen = [gen1, gen2, gen3] + ! Build isymax(:, 1) + ! determine to which irrep the translations belong to + ! Loop over Cartesian axes + do i = 1, 3 + bitmap = 0 + ! Loop over generators + do j = 1, nr_gen + ! Apply generators on Cartesian axes rots(i) and check the character + if (nint(get_pt(ior(igen(j), rots(i)))) == -1) then + ! Set the bitmap + bitmap = ibset(bitmap, j) + end if + end do + ! Right-shift the bitmap and assign to isymax + isymax(i, 1) = ishft(bitmap, -1) + end do + + ! Build isymax(:, 2) + ! determine to which irrep the rotations belong to + ! R_x = (y XOR z) and cyclic permutations + isymax(1, 2) = ieor(isymax(2, 1), isymax(3, 1)) + isymax(2, 2) = ieor(isymax(3, 1), isymax(1, 1)) + isymax(3, 2) = ieor(isymax(1, 1), isymax(2, 1)) + + ! Build the character table + lsymop = .false. + ! Activate all symmetry operations of the group + lsymop(0) = .true. + jsop(0) = 0 + ipar(0) = 1 + do i = 1, maxrep + i0 = iand(1_regint_k, i) * igen(1) + i1 = iand(1_regint_k, ishft(i, -1)) * igen(2) + i2 = iand(1_regint_k, ishft(i, -2)) * igen(3) + ind = ieor(ieor(i0, i1),i2) + lsymop(ind) = .true. + ipar(i) = jpar(ind) + end do + ! List group operations in preferred order + ! Identity, E + ind = 0 + jsop(ind) = 0 + ! Rotations + nrots = 0 + do i = 1, 3 + if (lsymop(irots(i))) then + ind = ind + 1 + jsop(ind) = irots(i) + nrots = nrots + 1 + end if + end do + ! Inversion + ninvc = 0 + if (lsymop(7)) then + ind = ind + 1 + jsop(ind) = 7 + ninvc = 1 + end if + ! Reflections + nrefl = 0 + do i = 1, 3 + if (lsymop(irefl(i))) then + ind = ind + 1 + jsop(ind) = irefl(i) + nrefl = nrefl + 1 + end if + end do + ! Classify group + ! ============== + ! tsaue - Here I have devised a highly empirical formula, but it works !!! + igroup = min(7, nint((4 * nrots + 8 * ninvc + 6 * nrefl) / 3.0)) + group = groups(igroup) + char_tab = 0 + ! Now generate the character table + do i = 0, maxrep + ! The character of the identity is always +1 + char_tab(0, i) = 1 + do j = 1, nr_gen + char_tab(igen(j), i) = nint(get_pt(iand(ishft(i,-(j-1)), 1_regint_k))) + do k = 1, (j-1) + ind = ieor(igen(j), igen(k)) + char_tab(ind, i) = char_tab(igen(j), i) * char_tab(igen(k), i) + do l = 1, (k-1) + char_tab(ieor(ind, igen(l)), i) = char_tab(ind, i) * char_tab(igen(l), i) + end do + end do + end do + end do + ! Classify irrep + do i = 0, maxrep + rep(i) = 'A ' + ipos = 2 + ! 1. Rotational symmetry + if (nrots == 3) then + ind = (1 - char_tab(jsop(1), i)) + (1 - char_tab(jsop(2), i))/2 + if (ind /= 0) then + rep(i)(1:1) = 'B' + rep(i)(2:2) = char(ichar('0') + ind) + ipos = 3 + end if + else if (nrots == 1) then + if (char_tab(jsop(1), i) == -1) then + rep(i)(1:1) = 'B' + end if + if (nrefl == 2) then + if (iand(ishft(jsop(1), -1), 1_regint_k) == 1) then + ind = 2 + else + ind = 3 + end if + if (char_tab(jsop(ind), i) == 1) then + rep(i)(2:2) = '1' + else + rep(i)(2:2) = '2' + end if + end if + else if (nrefl == 1) then + ! 2. Mirror symmetry + if (char_tab(jsop(1), i) == 1) then + rep(i)(2:2) = '''' + else if (char_tab(jsop(1), i) == -1) then + rep(i)(2:2) = '"' + end if + end if + ! 3. Inversion symmetry + if (ninvc == 1) then + ind = nrots + 1 + if (char_tab(jsop(ind), i) == 1) then + rep(i)(ipos:ipos) = 'g' + else + rep(i)(ipos:ipos) = 'u' + end if + end if + end do + ! Output + ! 1. Group name and generators + write(print_out, '(a, a3)') 'Point group: ', group + if (nr_gen > 0) then + write(print_out, '(/3x, a/)') '* The point group was generated by:' + do i = 1, nr_gen + if (symop(igen(i))(1:1) == 'C') then + write(print_out, '(6x, 3a)') 'Rotation about the ', symop(igen(i))(3:3),'-axis' + else if (symop(igen(i))(1:1) == 'O') then + write(print_out, '(6x, 3a)') 'Reflection in the ', symop(igen(i))(2:3),'-plane' + else + write(print_out, '(6x, a)') 'Inversion center' + end if + end do + ! 2. Group multiplication table + write(print_out,'(/3x, a/)') '* Group multiplication table' + write(print_out,'(8x, a1, 8(1x, a3, 1x))') '|', (symop(jsop(i)), i = 0, maxrep) + write(print_out,'(3x,a6,8a5)') '-----+', ('-----', i = 0, maxrep) + do i = 0, maxrep + write(print_out,'(4x, a3, 1x, a1, 8(1x, a3, 1x))') symop(jsop(i)), '|', & + (symop(ieor(jsop(i), jsop(j))), j = 0, maxrep) + end do + ! 3. Character table + write(print_out,'(/3x, a/)') '* Character table' + write(print_out,'(8x, a1, 8(1x, a3, 1x))') '|', (symop(jsop(j)), j = 0, maxrep) + write(print_out,'(3x, a6, 8a5)') '-----+', ('-----', i = 0, maxrep) + do i = 0, maxrep + write(print_out,'(4x, a3, 1x, a1, 8(1x, i3, 1x))') rep(i), '|', (char_tab(jsop(j), i), j=0, maxrep) + end do + ! 4. Direct product table + write(print_out,'(/3x, a/)') '* Direct product table' + write(print_out,'(8x, a1, 8(1x, a3, 1x))') '|', (rep(i), i = 0, maxrep) + write(print_out,'(3x, a6, 8a5)') '-----+', ('-----', i = 0, maxrep) + do i = 0, maxrep + write(print_out,'(3x, 1x, a3, 1x, a1, 8(1x, a3, 1x))') rep(i), '|', (rep(ieor(i, j)), j = 0, maxrep) + end do + end if + ! Fields: group name, group integer(kind=regint_k), number of generators, + ! number of nontrivial operations, isymax, jsop, + ! number of rotations, number of reflections, + ! number of inversions. + pg = point_group(group, igroup, nr_gen, maxrep, isymax, jsop, nrots, nrefl, ninvc) + + end function build_point_group + + end module tsless_symmetry diff --git a/src/tsless/tsless_weight_function.f90 b/src/tsless/tsless_weight_function.f90 new file mode 100644 index 0000000..a84663e --- /dev/null +++ b/src/tsless/tsless_weight_function.f90 @@ -0,0 +1,248 @@ +! +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. +! +! This file is part of PCMSolver. +! +! PCMSolver is free software: you can 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 3 of the License, or +! (at your option) any later version. +! +! PCMSolver 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with PCMSolver. If not, see . +! +! For information on the complete list of contributors to the +! PCMSolver API, see: +! + +module tsless_weight_function + +use, intrinsic :: iso_c_binding +use tsless_precision + +implicit none + +!> \brief Describes a weight function +!> Data fiels have dimensions 2*nr_derivative+3 +type, public :: weight_function + !> Number of continuous derivatives + integer(kind=regint_k) :: nr_derivative + !> + real(kind=dp), pointer :: weight_1(:) + !> + real(kind=dp), pointer :: weight_2(:) + !> + real(kind=dp), pointer :: weight_3(:) +end type weight_function + +public create_weight_function +public destroy_weight_function +public scaling_factor + +private + +contains + + function create_weight_function(nderiv, normalization, dmin, printer) result(wfun) + + use, intrinsic :: iso_fortran_env, only: output_unit + + !> Passed variables + integer(kind=regint_k), intent(in) :: nderiv + integer(kind=regint_k), intent(in) :: normalization + real(kind=dp), intent(in) :: dmin + integer, optional, intent(in) :: printer + !> Output variables + type(weight_function) :: wfun + !> Local variables + integer :: print_out + + if (present(printer)) then + print_out = printer + else + print_out = output_unit + end if + + call allocate_weight_function(wfun, nderiv) + call compute(wfun, normalization, dmin, print_out) + + end function create_weight_function + + pure subroutine allocate_weight_function(wfun, nderiv) + + !> Passed variables + type(weight_function), intent(inout) :: wfun + integer(kind=regint_k), intent(in) :: nderiv + !> Local variables + integer(kind=regint_k) :: nparms + + wfun%nr_derivative = nderiv + nparms = 2_regint_k * nderiv + 3_regint_k + allocate(wfun%weight_1(nparms)); wfun%weight_1 = 0.0_dp + allocate(wfun%weight_2(nparms)); wfun%weight_2 = 0.0_dp + allocate(wfun%weight_3(nparms)); wfun%weight_3 = 0.0_dp + + end subroutine allocate_weight_function + + pure subroutine destroy_weight_function(wfun) + + !> Passed variables + type(weight_function), intent(inout) :: wfun + + wfun%nr_derivative = 0_regint_k + deallocate(wfun%weight_1) + deallocate(wfun%weight_2) + deallocate(wfun%weight_3) + + end subroutine destroy_weight_function + + !> \brief + !> \author Christian S. Pomelli + !> \param[inout] wfun weight function + !> \param[in] nord maximum number of continuous derivatives + !> \param[in] ifun + !> \param[in] dmin minimal distance between sampling points + !> \param[in] printer logical unit for printing + !> ifun = 0 non-normalized function (gamma) + !> ifun = 1 normalized function (omega) + subroutine compute(wfun, ifun, dmin, printer) + + use, intrinsic :: iso_fortran_env, only: output_unit + + use tsless_lapack, only: solve_linear_system + + !> Passed variables + type(weight_function), intent(inout) :: wfun + integer(kind=regint_k), intent(in) :: ifun + real(kind=dp), intent(in) :: dmin + integer, optional, intent(in) :: printer + !> Local variables + real(kind=dp), allocatable :: C(:, :) + real(kind=dp), allocatable :: b(:) + integer :: i, j, k, n, m + integer :: print_out + + if (present(printer)) then + print_out = printer + else + print_out = output_unit + end if + + ! check parameters + if ((ifun .lt. 0) .or. (ifun .gt. 1)) Then + write(print_out, *) 'Unknown IFUN' + stop + end if + + ! m = number of parameters + n = wfun%nr_derivative + if (ifun .eq. 0) then + m = 2*n + 2 + else + m = 2*n + 3 + end if + !> Allocate, clean up and fill coefficient matrix C + allocate(C(m, m)); C = 0.0_dp + ! l = 1 boundary conditions + do j = 0_regint_k, n + do k = j, m - 1_regint_k + C(j+1, k+1) = coefficients(j, k) + end do + end do + ! l = d boundary conditions + do j = 0_regint_k, n + do k = j, m - 1_regint_k + C(j+2+n, k+1) = coefficients(j, k) + C(j+2+n, k+1) = C(j+2+n, k+1) * dmin**(k-j) + end do + end do + ! normalization condition + if (ifun .eq. 1_regint_k) then + do k = 0_regint_k, m - 1_regint_k + C(m, k+1) = (1.0_dp - dmin**(k+1)) / (1.0_dp * (k+1)) + end do + end if + ! Allocate, clean up and fill b vector (RHS) + allocate(b(m)); b = 0.0_dp + b(1) = 1.0_dp + if (ifun .eq. 1) b(m) = 1.0_dp - dmin + ! 2. Solve linear system of equations + call solve_linear_system(ndim=m, A=C, b=b, x=wfun%weight_1) + ! write function and derivatives + do i = 1_regint_k, m - 1_regint_k + wfun%weight_2(i) = i * wfun%weight_1(i+1) + end do + + do i = 1_regint_k, m - 2_regint_k + wfun%weight_3(i) = i * wfun%weight_2(i+1) + end do + + write(print_out, *) 'ifun = ', ifun + do i = 1_regint_k, m + write(print_out, '(a, i2, a, 3F15.4)') 'c(', i, ') = ', wfun%weight_1(i), wfun%weight_2(i), wfun%weight_3(i) + end do + write(print_out, *) '-------------------' + + deallocate(b) + deallocate(C) + + end subroutine compute + + !> \brief Calculates scaling factor for weight of point + !> \author Christian S. Pomelli + !> \param[in] wfun weight function + !> \param[in] arg argument of the regularized step function + !> \param[in] dmin minimal distance between sampling points + pure function scaling_factor(wfun, arg, dmin) result(scaling) + + !> Input variables + type(weight_function), intent(in) :: wfun + real(kind=dp), intent(in) :: arg + real(kind=dp), intent(in) :: dmin + !> Output variables + real(kind=dp) :: scaling(3) + !> Local variables + integer :: i, nr_params + real(kind=dp) :: tn + + ! Simply a polynomial. later add derivatives + scaling = 0.0_dp + tn = 1.0_dp + ! number of polynomial parameters + nr_params = 2_regint_k * wfun%nr_derivative + 3_regint_k + do i = 1_regint_k, nr_params + scaling(1) = scaling(1) + wfun%weight_1(i) * tn + scaling(2) = scaling(2) + wfun%weight_2(i) * tn + scaling(3) = scaling(3) + wfun%weight_3(i) * tn + tn = tn * arg + end do + + end function scaling_factor + + pure function coefficients(j, k) result(f) + + !> Passed variables + integer(kind=regint_k), intent(in) :: j, k + real(kind=dp) :: f + !> Local variables + integer :: i + real(kind=dp) :: dkk + + f = 1.0_dp + dkk = dble(k) + if (j .ne. 0) then + do i = 1, j + f = f * dkk + dkk = dkk - 1.0_dp + end do + end if + + end function coefficients + +end module tsless_weight_function diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index cf0f920..7182305 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -81,6 +81,7 @@ add_subdirectory(bi_operators ${PARENT_BINARY_DIR}/tests/bi_operator) add_subdirectory(cpcm ${PARENT_BINARY_DIR}/tests/cpcm) add_subdirectory(iefpcm ${PARENT_BINARY_DIR}/tests/iefpcm) add_subdirectory(utils ${PARENT_BINARY_DIR}/tests/utils) +add_subdirectory(tsless ${PARENT_BINARY_DIR}/tests/tsless) # Little hack to avoid prefixing every include in every source file with "PCMSolver" target_include_directories(unit_tests diff --git a/tests/TestingMolecules.hpp b/tests/TestingMolecules.hpp index 38e4b3d..6db2d68 100644 --- a/tests/TestingMolecules.hpp +++ b/tests/TestingMolecules.hpp @@ -35,6 +35,10 @@ #include "utils/Symmetry.hpp" namespace pcm { +/*! Returns the hydrogen fluoride molecule + */ +inline Molecule HF(); + /*! Returns the ammonia molecule */ inline Molecule NH3(); @@ -123,6 +127,34 @@ template Molecule dummy(double radius, const Eigen::Vector3d & cente 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 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 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() { int nAtoms = 4; diff --git a/tests/cpcm/CMakeLists.txt b/tests/cpcm/CMakeLists.txt index 526bfc7..feae426 100644 --- a/tests/cpcm/CMakeLists.txt +++ b/tests/cpcm/CMakeLists.txt @@ -7,6 +7,8 @@ target_sources(unit_tests ${CMAKE_CURRENT_SOURCE_DIR}/cpcm_gepol-point.cpp ${CMAKE_CURRENT_SOURCE_DIR}/cpcm_gepol-point_from-file.cpp ${CMAKE_CURRENT_SOURCE_DIR}/cpcm_symmetry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/cpcm_tsless-point.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/cpcm_tsless-NH3.cpp ) # cpcm_gepol-NH3_from-file.cpp test @@ -94,3 +96,22 @@ add_Catch_test( REFERENCE_FILES ASC-cpcm_gepol-C2H4_D2h.npy ) + +# FIXME Reference files for these two tests +# cpcm_tsless-point.cpp test +add_Catch_test( + NAME + cpcm_tsless-point + LABELS + cpcm + cpcm_tsless-point + ) + +# cpcm_tsless-NH3.cpp test +add_Catch_test( + NAME + cpcm_tsless-NH3 + LABELS + cpcm + cpcm_tsless-NH3 + ) diff --git a/tests/cpcm/cpcm_tsless-NH3.cpp b/tests/cpcm/cpcm_tsless-NH3.cpp new file mode 100644 index 0000000..8c05cf7 --- /dev/null +++ b/tests/cpcm/cpcm_tsless-NH3.cpp @@ -0,0 +1,84 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.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" + +using namespace pcm; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +using solver::CPCMSolver; + +/*! \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)); +} diff --git a/tests/cpcm/cpcm_tsless-point.cpp b/tests/cpcm/cpcm_tsless-point.cpp new file mode 100644 index 0000000..5fe1a7d --- /dev/null +++ b/tests/cpcm/cpcm_tsless-point.cpp @@ -0,0 +1,138 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.hpp" +#include "green/DerivativeTypes.hpp" +#include "green/UniformDielectric.hpp" +#include "green/Vacuum.hpp" +#include "solver/CPCMSolver.hpp" + +using namespace pcm; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +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)); + } + } + + /*! \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)); + } + } + } +} diff --git a/tests/gepol/CMakeLists.txt b/tests/gepol/CMakeLists.txt index 4751574..3957562 100644 --- a/tests/gepol/CMakeLists.txt +++ b/tests/gepol/CMakeLists.txt @@ -12,6 +12,7 @@ target_sources(unit_tests ${CMAKE_CURRENT_SOURCE_DIR}/gepol_point.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gepol_point_from-file.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gepol_point_symmetry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/gepol_HF.cpp ) # gepol_point_symmetry.cpp test @@ -127,3 +128,12 @@ add_Catch_test( DEPENDS gepol_NH3 ) + +# gepol_HF.cpp test +add_Catch_test( + NAME + gepol_HF + LABELS + gepol + gepol_HF + ) diff --git a/tests/gepol/gepol_HF.cpp b/tests/gepol/gepol_HF.cpp new file mode 100644 index 0000000..5919f12 --- /dev/null +++ b/tests/gepol/gepol_HF.cpp @@ -0,0 +1,76 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "cavity/GePolCavity.hpp" +#include "utils/Molecule.hpp" +#include "utils/Symmetry.hpp" + +using namespace pcm; +using cavity::GePolCavity; + +TEST_CASE("GePol cavity for an hydrogen fluoride molecule", "[gepol][gepol_HF]") { + Molecule molec = HF(); + double area = 0.02 / bohr2ToAngstrom2(); + double probeRadius = 0.0 / bohrToAngstrom(); + double minRadius = 0.2 / bohrToAngstrom(); + GePolCavity cavity(molec, area, probeRadius, minRadius); + + /*! \class GePolCavity + * \test \b GePolCavityHFTest_size tests GePol cavity size for ammonia + */ + SECTION("Test size") { + size_t ref_size = 1688; + size_t size = cavity.size(); + REQUIRE(size == ref_size); + } + + /*! \class GePolCavity + * \test \b GePolCavityHFTest_area tests GePol cavity surface area for ammonia + */ + SECTION("Test surface area") { + double ref_area = 110.64517236179323; + double area = cavity.elementArea().sum(); + REQUIRE(area == Approx(ref_area)); + } + + /*! \class GePolCavity + * \test \b GePolCavityHFTest_volume tests GePol cavity volume for ammonia + */ + SECTION("Test volume") { + double ref_volume = 105.98002389543836; + 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)); + } +} diff --git a/tests/iefpcm/CMakeLists.txt b/tests/iefpcm/CMakeLists.txt index e619222..55fb1f5 100644 --- a/tests/iefpcm/CMakeLists.txt +++ b/tests/iefpcm/CMakeLists.txt @@ -9,6 +9,8 @@ target_sources(unit_tests ${CMAKE_CURRENT_SOURCE_DIR}/iefpcm_gepol-point.cpp ${CMAKE_CURRENT_SOURCE_DIR}/iefpcm_gepol-point_from-file.cpp ${CMAKE_CURRENT_SOURCE_DIR}/iefpcm_symmetry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/iefpcm_tsless-point.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/iefpcm_tsless-NH3.cpp ) # iefpcm_gepol-NH3.cpp test @@ -107,3 +109,23 @@ add_Catch_test( DEPENDS iefpcm_gepol-point ) + +# iefpcm_tsless-point.cpp test +add_Catch_test( + NAME + iefpcm_tsless-point + LABELS + solver + iefpcm + iefpcm_tsless-point + ) + +# iefpcm_tsless-NH3.cpp test +add_Catch_test( + NAME + iefpcm_tsless-NH3 + LABELS + solver + iefpcm + iefpcm_tsless-NH3 + ) diff --git a/tests/iefpcm/iefpcm_tsless-NH3.cpp b/tests/iefpcm/iefpcm_tsless-NH3.cpp new file mode 100644 index 0000000..9a06607 --- /dev/null +++ b/tests/iefpcm/iefpcm_tsless-NH3.cpp @@ -0,0 +1,81 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.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; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +using solver::IEFSolver; + +/*! \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)); +} diff --git a/tests/iefpcm/iefpcm_tsless-point.cpp b/tests/iefpcm/iefpcm_tsless-point.cpp new file mode 100644 index 0000000..82f8261 --- /dev/null +++ b/tests/iefpcm/iefpcm_tsless-point.cpp @@ -0,0 +1,136 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#include "TestingMolecules.hpp" +#include "bi_operators/Collocation.hpp" +#include "cavity/TsLessCavity.hpp" +#include "green/DerivativeTypes.hpp" +#include "green/UniformDielectric.hpp" +#include "green/Vacuum.hpp" +#include "solver/IEFSolver.hpp" + +using namespace pcm; +using bi_operators::Collocation; +using cavity::TsLessCavity; +using green::UniformDielectric; +using green::Vacuum; +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)); + } + } + } +} diff --git a/tests/tsless/CMakeLists.txt b/tests/tsless/CMakeLists.txt new file mode 100644 index 0000000..b56ae67 --- /dev/null +++ b/tests/tsless/CMakeLists.txt @@ -0,0 +1,34 @@ +target_sources(unit_tests + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_HF.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_NH3.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/tsless_point.cpp + ) + +# tsless_HF.cpp test +add_Catch_test( + NAME + tsless_HF + LABELS + tsless + tsless_HF + ) + +# tsless_point.cpp test +add_Catch_test( + NAME + tsless_point + LABELS + tsless + tsless_point + ) + +# tsless_NH3.cpp test +add_Catch_test( + NAME + tsless_NH3 + LABELS + tsless + tsless_NH3 + ) + diff --git a/tests/tsless/tsless_HF.cpp b/tests/tsless/tsless_HF.cpp new file mode 100644 index 0000000..c0245f3 --- /dev/null +++ b/tests/tsless/tsless_HF.cpp @@ -0,0 +1,76 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#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)); + } +} diff --git a/tests/tsless/tsless_NH3.cpp b/tests/tsless/tsless_NH3.cpp new file mode 100644 index 0000000..ec611bb --- /dev/null +++ b/tests/tsless/tsless_NH3.cpp @@ -0,0 +1,77 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#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)); + } +} diff --git a/tests/tsless/tsless_point.cpp b/tests/tsless/tsless_point.cpp new file mode 100644 index 0000000..6d032e1 --- /dev/null +++ b/tests/tsless/tsless_point.cpp @@ -0,0 +1,76 @@ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors. + * + * This file is part of PCMSolver. + * + * PCMSolver is free software: you can 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 3 of the License, or + * (at your option) any later version. + * + * PCMSolver 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with PCMSolver. If not, see . + * + * For information on the complete list of contributors to the + * PCMSolver API, see: + */ + +#include "catch.hpp" + +#include + +#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)); + } +}