From 630c65786f104bacbf64bbfd5cc2efc5dabac355 Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Mon, 5 Oct 2015 14:41:52 +0200 Subject: [PATCH 1/6] IMPORTED. This is an artificial re-creation of the Rob-tsless branch. --- src/CMakeLists.txt | 1 + src/bin/run_pcm.cpp | 2 +- src/cavity/CMakeLists.txt | 5 +- src/cavity/GePolCavity.cpp | 100 +++---- src/cavity/GePolCavity.hpp | 4 +- src/cavity/RegisterCavityToFactory.hpp | 23 +- src/cavity/TsLessCavity.cpp | 275 +++++++++++++++++ src/cavity/TsLessCavity.hpp | 89 ++++++ src/pedra/pedra_cavity_interface.F90 | 31 +- src/tsless/CMakeLists.txt | 11 + src/tsless/tsless_cavity.F90 | 390 +++++++++++++++++++++++++ src/tsless/tsless_lapack.F90 | 158 ++++++++++ src/tsless/tsless_precision.F90 | 42 +++ src/tsless/tsless_symmetry.F90 | 354 ++++++++++++++++++++++ src/tsless/tsless_weight_function.F90 | 251 ++++++++++++++++ src/utils/Citation.hpp | 4 +- src/utils/Molecule.cpp | 20 +- src/utils/TestingMolecules.hpp | 44 ++- tests/CMakeLists.txt | 1 + tests/cpcm/CMakeLists.txt | 4 + tests/cpcm/cpcm_gepol-NH3.cpp | 12 +- tests/cpcm/cpcm_gepol-point.cpp | 10 +- tests/cpcm/cpcm_tsless-NH3.cpp | 76 ++--- tests/cpcm/cpcm_tsless-point.cpp | 155 ++++++---- tests/gepol/gepol_HF.cpp | 80 +++++ tests/iefpcm/CMakeLists.txt | 4 + tests/iefpcm/iefpcm_gepol-point.cpp | 12 +- tests/iefpcm/iefpcm_tsless-NH3.cpp | 61 ++-- tests/iefpcm/iefpcm_tsless-point.cpp | 147 +++++++--- tests/tsless/CMakeLists.txt | 6 + tests/tsless/tsless_HF.cpp | 82 ++++++ tests/tsless/tsless_NH3.cpp | 83 ++++++ tests/tsless/tsless_point.cpp | 79 +++++ tools/pcmsolver.py.in | 13 +- 34 files changed, 2321 insertions(+), 308 deletions(-) create mode 100644 src/cavity/TsLessCavity.cpp create mode 100644 src/cavity/TsLessCavity.hpp create mode 100644 src/tsless/CMakeLists.txt create mode 100644 src/tsless/tsless_cavity.F90 create mode 100644 src/tsless/tsless_lapack.F90 create mode 100644 src/tsless/tsless_precision.F90 create mode 100644 src/tsless/tsless_symmetry.F90 create mode 100644 src/tsless/tsless_weight_function.F90 create mode 100644 tests/gepol/gepol_HF.cpp create mode 100644 tests/tsless/CMakeLists.txt create mode 100644 tests/tsless/tsless_HF.cpp create mode 100644 tests/tsless/tsless_NH3.cpp create mode 100644 tests/tsless/tsless_point.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bda52c167..006da8cbb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,6 +14,7 @@ add_subdirectory(metal) add_subdirectory(bi_operators) add_subdirectory(pedra) add_subdirectory(solver) +add_subdirectory(tsless) add_subdirectory(utils) get_property(_pcmsolver_cxx_sources GLOBAL PROPERTY PCMSolver_CXX_SOURCES) diff --git a/src/bin/run_pcm.cpp b/src/bin/run_pcm.cpp index 4e8ea911d..6ebf56cf0 100644 --- a/src/bin/run_pcm.cpp +++ b/src/bin/run_pcm.cpp @@ -111,7 +111,7 @@ int main(int argc, char * argv[]) // Compute apparent surface charge Eigen::VectorXd asc = solver->computeCharge(mep); // Compute energy and print it out - out_stream << "Solvation energy = " << 0.5 * (asc * mep) << std::endl; + out_stream << "Solvation energy = " << std::setprecision(14) << 0.5 * asc.dot(mep) << std::endl; out_stream << "DONE!" << std::endl; out_stream.close(); diff --git a/src/cavity/CMakeLists.txt b/src/cavity/CMakeLists.txt index 379876635..bd407771a 100644 --- a/src/cavity/CMakeLists.txt +++ b/src/cavity/CMakeLists.txt @@ -1,8 +1,7 @@ -# List of headers -list(APPEND headers_list Cavity.hpp Element.hpp GePolCavity.hpp RegisterCavityToFactory.hpp RestartCavity.hpp) +list(APPEND headers_list Cavity.hpp Element.hpp GePolCavity.hpp RegisterCavityToFactory.hpp RestartCavity.hpp TsLessCavity.hpp) # List of sources -list(APPEND sources_list Cavity.cpp Element.cpp GePolCavity.cpp RestartCavity.cpp) +list(APPEND sources_list Cavity.cpp Element.cpp GePolCavity.cpp RestartCavity.cpp TsLessCavity.cpp) # Update bar chart if(BUILD_CHARTS) diff --git a/src/cavity/GePolCavity.cpp b/src/cavity/GePolCavity.cpp index 7ae0fb9ed..c032383b8 100644 --- a/src/cavity/GePolCavity.cpp +++ b/src/cavity/GePolCavity.cpp @@ -31,6 +31,7 @@ #include #include "Config.hpp" +#include "FCMangle.hpp" #include #include @@ -38,14 +39,7 @@ #include "Sphere.hpp" #include "Symmetry.hpp" -/*! \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 @@ -77,20 +71,22 @@ * \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); +#define pedra_driver\ + FortranCInterface_GLOBAL_(pedra_driver, PEDRA_DRIVER) +extern "C" void pedra_driver(size_t * maxts, size_t * maxsph, size_t * 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); -void GePolCavity::build(int maxts, int maxsph, int maxvert) +void GePolCavity::build(size_t maxts, size_t maxsph, size_t maxvert) { - // This is a wrapper for the generatecavity_cpp function defined in the Fortran code PEDRA. + // 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 // for the insertion of additional spheres as in the most general formulation of the // GePol algorithm. @@ -168,7 +164,7 @@ void GePolCavity::build(int maxts, int maxsph, int maxvert) // Go PEDRA, Go! TIMER_ON("GePolCavity::generatecavity_cpp"); - generatecavity_cpp(&maxts, &maxsph, &maxvert, + pedra_driver(&maxts, &maxsph, &maxvert, xtscor, ytscor, ztscor, ar, xsphcor, ysphcor, zsphcor, rsph, &nts, &ntsirr, &nSpheres_, &addedSpheres, xe, ye, ze, rin, mass, @@ -254,31 +250,31 @@ void GePolCavity::build(int maxts, int maxsph, int maxvert) // Fill elements_ vector for (int i = 0; i < nElements_; ++i) { - int i_off = i + 1; - bool irr = false; - // PEDRA puts the irreducible tesserae first - if (i < nIrrElements_) irr = true; - Sphere sph(elementSphereCenter_.col(i), elementRadius_(i)); - int nv = nvert[i]; - Eigen::Matrix3Xd vertices, arcs; - vertices.resize(Eigen::NoChange, nv); - arcs.resize(Eigen::NoChange, nv); - // Populate vertices and arcs - for (int j = 0; j < nv; ++j) { - int j_off = (j + 1) * nElements_ - 1; - for (int k = 0; k < 3; ++k) { - int k_off = (k + 1) * nElements_ * nv; - int offset = i_off + j_off + k_off; - vertices(k, j) = vert[offset]; - arcs(k, j) = centr[offset]; - } - } - elements_.push_back(Element(nv, - elementArea_(i), - elementCenter_.col(i), - elementNormal_.col(i), - irr, sph, - vertices, arcs)); + int i_off = i + 1; + bool irr = false; + // PEDRA puts the irreducible tesserae first + if (i < nIrrElements_) irr = true; + Sphere sph(elementSphereCenter_.col(i), elementRadius_(i)); + int nv = nvert[i]; + Eigen::Matrix3Xd vertices, arcs; + vertices.resize(Eigen::NoChange, nv); + arcs.resize(Eigen::NoChange, nv); + // Populate vertices and arcs + for (int j = 0; j < nv; ++j) { + int j_off = (j + 1) * nElements_ - 1; + for (int k = 0; k < 3; ++k) { + int k_off = (k + 1) * nElements_ * nv; + int offset = i_off + j_off + k_off; + vertices(k, j) = vert[offset]; + arcs(k, j) = centr[offset]; + } + } + elements_.push_back(Element(nv, + elementArea_(i), + elementCenter_.col(i), + elementNormal_.col(i), + irr, sph, + vertices, arcs)); } // Clean-up @@ -313,15 +309,15 @@ std::ostream & GePolCavity::printCavity(std::ostream & os) os << "\nNumber of irreducible finite elements = " << nIrrElements_; } /* - for(int i = 0; i < nElements_; i++) + for (int i = 0; i < nElements_; i++) { - os << std::endl; - os << i+1 << " "; - os << elementCenter_(0,i) << " "; - os << elementCenter_(1,i) << " "; - os << elementCenter_(2,i) << " "; - os << elementArea_(i) << " "; - } + os << std::endl; + os << i+1 << " "; + os << elementCenter_(0,i) << " "; + os << elementCenter_(1,i) << " "; + os << elementCenter_(2,i) << " "; + os << elementArea_(i) << " "; + } */ return os; } diff --git a/src/cavity/GePolCavity.hpp b/src/cavity/GePolCavity.hpp index 45610d3c8..ff73aa4ec 100644 --- a/src/cavity/GePolCavity.hpp +++ b/src/cavity/GePolCavity.hpp @@ -32,8 +32,6 @@ #include "Config.hpp" -#include - #include "Cavity.hpp" /*! \file GePolCavity.hpp @@ -82,7 +80,7 @@ class GePolCavity : public Cavity * \param[in] maxsp maximum number of spheres (original + added) * \param[in] maxvert maximum number of vertices */ - void build(int maxts, int maxsp, int maxvert); + void build(size_t maxts, size_t maxsp, size_t maxvert); }; #endif // GEPOLCAVITY_HPP diff --git a/src/cavity/RegisterCavityToFactory.hpp b/src/cavity/RegisterCavityToFactory.hpp index a32f548f5..271c61f3f 100644 --- a/src/cavity/RegisterCavityToFactory.hpp +++ b/src/cavity/RegisterCavityToFactory.hpp @@ -2,22 +2,22 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ @@ -32,6 +32,7 @@ #include "Factory.hpp" #include "GePolCavity.hpp" #include "RestartCavity.hpp" +#include "TsLessCavity.hpp" /*! \file RegisterCavityToFactory.hpp * \brief Register each cavity to the factory. @@ -67,4 +68,16 @@ namespace RESTART, createRestartCavity); } +namespace +{ + Cavity * createTsLessCavity(const cavityData & data) + { + return new TsLessCavity(data.molecule, data.area, data.probeRadius, data.minimalRadius, + data.minDistance, data.derOrder); + } + const std::string TSLESS("TSLESS"); + const bool registeredTsLess = Factory::TheFactory().registerObject(TSLESS, + createTsLessCavity); +} + #endif // REGISTERCAVITYTOFACTORY_HPP diff --git a/src/cavity/TsLessCavity.cpp b/src/cavity/TsLessCavity.cpp new file mode 100644 index 000000000..34afbd3c4 --- /dev/null +++ b/src/cavity/TsLessCavity.cpp @@ -0,0 +1,275 @@ +/* pcmsolver_copyright_start */ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2013 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: + */ +/* pcmsolver_copyright_end */ + +#include "TsLessCavity.hpp" + +#include +#include +#include +#include + +#include "Config.hpp" +#include "FCMangle.hpp" + +#include + +#include "Exception.hpp" +#include "Sphere.hpp" +#include "Symmetry.hpp" +#include "TimerInterface.hpp" + +/*! \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 + */ +#define tsless_driver \ + FortranCInterface_MODULE_(tsless_cavity, tsless_driver, TSLESS_CAVITY, TSLESS_DRIVER) +extern "C" void tsless_driver(size_t * maxts, size_t * maxsph, size_t * 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); + +void TsLessCavity::build(size_t maxts, size_t maxsph, size_t maxvert) +{ + // This is a wrapper for the generatecavity_cpp_ function defined in the Fortran code TsLess. + // Here we allocate the necessary arrays to be passed to TsLess, in particular we allow + // for the insertion of additional spheres as in the most general formulation of the + // GePol algorithm. + + int lwork = maxts*maxsph; + double * xtscor = new double[maxts]; + double * ytscor = new double[maxts]; + double * ztscor = new double[maxts]; + double * ar = new double[maxts]; + double * xsphcor = new double[maxts]; + double * ysphcor = new double[maxts]; + double * zsphcor = new double[maxts]; + double * rsph = new double[maxts]; + double * work = new double[lwork]; + + // Clean-up possible heap-crap + std::fill_n(xtscor, maxts, 0.0); + std::fill_n(ytscor, maxts, 0.0); + std::fill_n(ztscor, maxts, 0.0); + std::fill_n(ar, maxts, 0.0); + std::fill_n(xsphcor, maxts, 0.0); + std::fill_n(ysphcor, maxts, 0.0); + std::fill_n(zsphcor, maxts, 0.0); + std::fill_n(rsph, maxts, 0.0); + std::fill_n(work, lwork, 0.0); + + int nts = 0; + int ntsirr = 0; + + // If there's an overflow in the number of spheres TsLess will die. + // The maximum number of spheres in PEDRA is set to 200 (primitive+additional) + // so the integer here declared is just to have enough space C++ side to pass everything back. + int maxAddedSpheres = 200; + + // Allocate vectors of size equal to nSpheres_ + maxAddedSpheres where maxAddedSpheres is the + // maximum number of spheres we allow the algorithm to add to our original set. + // If this number is exceeded, then the algorithm crashes (should look into this...) + // After the cavity is generated we will update ALL the class data members, both related + // to spheres and finite elements so that the cavity is fully formed. + + Eigen::VectorXd xv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); + Eigen::VectorXd yv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); + Eigen::VectorXd zv = Eigen::VectorXd::Zero(nSpheres_ + maxAddedSpheres); + Eigen::VectorXd radii_scratch = Eigen::VectorXd::Zero(nSpheres_ + + maxAddedSpheres); // Not to be confused with the data member inherited from Cavity!!! + + for ( int i = 0; i < nSpheres_; ++i ) { + for ( int j = 0; j < 3; ++j ) { + xv(i) = sphereCenter_(0, i); + yv(i) = sphereCenter_(1, i); + zv(i) = sphereCenter_(2, i); + } + radii_scratch(i) = sphereRadius_(i); + } + + double * xe = xv.data(); + double * ye = yv.data(); + double * ze = zv.data(); + + double *rin = radii_scratch.data(); + double * mass = new double[molecule_.nAtoms()]; + for (size_t i = 0; i < molecule_.nAtoms(); ++i) { + mass[i] = molecule_.masses(i); + } + + addedSpheres = 0; + // Number of generators and generators of the point group + int nr_gen = molecule_.pointGroup().nrGenerators(); + int gen1 = molecule_.pointGroup().generators(0); + int gen2 = molecule_.pointGroup().generators(1); + int gen3 = molecule_.pointGroup().generators(2); + + int weightFunction = 1; + // Go TsLess, Go! + TIMER_ON("TsLessCavity::tsless_driver"); + tsless_driver(&maxts, &maxsph, &maxvert, &nSpheres_, &nts, &ntsirr, &addedSpheres, + xtscor, ytscor, ztscor, ar, + xsphcor, ysphcor, zsphcor, rsph, + xe, ye, ze, rin, mass, + &nr_gen, &gen1, &gen2, &gen3, + &averageArea_, &minDistance_, &derOrder_, &weightFunction, &probeRadius_, work); + TIMER_OFF("TsLessCavity::tsless_driver"); + + // The "intensive" part of updating the spheres related class data members will be of course + // executed iff addedSpheres != 0 + if (addedSpheres != 0) { + // Save the number of original spheres + int orig = nSpheres_; + // Update the nSpheres + nSpheres_ += addedSpheres; + // Resize sphereRadius and sphereCenter... + sphereRadius_.resize(nSpheres_); + sphereCenter_.resize(Eigen::NoChange, nSpheres_); + // Transfer radii and centers. + // Eigen has no push_back function, so we need to traverse all the spheres... + sphereRadius_ = radii_scratch.head(nSpheres_); + for ( int i = 0; i < nSpheres_; ++i ) { + sphereCenter_(0, i) = xv(i); + sphereCenter_(1, i) = yv(i); + sphereCenter_(2, i) = zv(i); + } + // Now grow the vector 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); + elements_.push_back(Element(nv, + 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_; + /* + for (int i = 0; i < nElements_; i++) + { + os << std::endl; + os << i+1 << " "; + os << elementCenter_(0,i) << " "; + os << elementCenter_(1,i) << " "; + os << elementCenter_(2,i) << " "; + os << elementArea_(i) << " "; + } + */ + return os; +} diff --git a/src/cavity/TsLessCavity.hpp b/src/cavity/TsLessCavity.hpp new file mode 100644 index 000000000..537db5bbc --- /dev/null +++ b/src/cavity/TsLessCavity.hpp @@ -0,0 +1,89 @@ +/* pcmsolver_copyright_start */ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2013 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: + */ +/* pcmsolver_copyright_end */ + +#ifndef TSLESSCAVITY_HPP +#define TSLESSCAVITY_HPP + +#include +#include +#include + +#include "Config.hpp" + +#include "Cavity.hpp" + +/*! \file TsLessCavity.hpp + * \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 : public Cavity +{ +public: + TsLessCavity() {} + TsLessCavity(const Molecule & molec, double a, double pr, double minR, double minD, int der) : + Cavity(molec), averageArea_(a), probeRadius_(pr), minimalRadius_(minR), minDistance_(minD), derOrder_(der) + { + std::string checkpointName = "TsLessCavity::build"; + TIMER_ON(checkpointName); + build(10000, 200, 25000); + TIMER_OFF(checkpointName); + } + TsLessCavity(const std::vector & sph, double a, double pr, double minR, double minD, int der) : + Cavity(sph), averageArea_(a), probeRadius_(pr), minimalRadius_(minR), minDistance_(minD), derOrder_(der) + { + std::string checkpointName = "TsLessCavity::build"; + TIMER_ON(checkpointName); + build(10000, 200, 25000); + TIMER_OFF(checkpointName); + } + 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); + virtual void makeCavity() { 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(size_t maxts, size_t maxsp, size_t maxvert); +}; + +#endif // TSLESSCAVITY_HPP diff --git a/src/pedra/pedra_cavity_interface.F90 b/src/pedra/pedra_cavity_interface.F90 index 9d4b26660..31e25753a 100644 --- a/src/pedra/pedra_cavity_interface.F90 +++ b/src/pedra/pedra_cavity_interface.F90 @@ -24,30 +24,29 @@ ! ! simple input reader for cavity generator ! written by Krzysztof Mozgawa, 2010 -! +! ! 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_) & - bind(c, name='generatecavity_cpp') + nvert_, vert_, centr_) -use, intrinsic :: iso_c_binding +use, intrinsic :: iso_c_binding use pedra_precision use pedra_symmetry, only: get_point_group, point_group use pedra_cavity, only: polyhedra_driver implicit none - + #include "pcm_pcmdef.h" #include "pcm_mxcent.h" #include "pcm_pcm.h" -integer(c_int) :: maxts_, maxsph_, maxvert_ +integer(c_size_t) :: maxts_, maxsph_, maxvert_ real(c_double) :: xtscor_(maxts_), ytscor_(maxts_), ztscor_(maxts_) real(c_double) :: xsphcor_(maxts_), ysphcor_(maxts_), zsphcor_(maxts_), rsph_(maxts_) real(c_double) :: ar_(maxts_), xe_(maxts_), ye_(maxts_), ze_(maxts_), rin_(maxts_) @@ -56,9 +55,9 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & integer(c_int) :: nts_, ntsirr_, nesfp_, addsph_ integer(c_int) :: nr_gen_, gen1_, gen2_, gen3_ integer(c_int) :: nvert_(maxts_) -real(c_double) :: vert_(maxts_ * 30), centr_(maxts_ * 30) +real(c_double) :: vert_(maxts_ * 30), centr_(maxts_ * 30) -integer(c_int) :: i, j, k, offset +integer(c_int) :: i, j, k, offset integer(c_int) :: error_code integer(kind=regint_k) :: lvpri logical :: pedra_file_exists @@ -69,10 +68,10 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & pedra_file_exists = .false. inquire(file = 'PEDRA.OUT', exist = pedra_file_exists) if (pedra_file_exists) then - open(lvpri, & - file = 'PEDRA.OUT', & - status = 'unknown', & - form = 'formatted', & + open(lvpri, & + file = 'PEDRA.OUT', & + status = 'unknown', & + form = 'formatted', & access = 'sequential') close(lvpri, status = 'delete') end if @@ -87,11 +86,11 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & iprpcm = 3 call get_point_group(lvpri, pgroup, nr_gen_, gen1_, gen2_, gen3_) rsolv = rsolv_ -! These parameters are fixed see one of the original GePol papers +! These parameters are fixed see one of the original GePol papers omega = 40.0d+00 fro = 0.7d+00 ! ret is the minimum radius of added spheres -ret = ret_ +ret = ret_ ! nesfp is the number of initial spheres. ! nesf is the total number of spheres: initial + added nesfp = nesfp_ @@ -149,4 +148,4 @@ subroutine generatecavity_cpp(maxts_, maxsph_, maxvert_, & close(lvpri) -end subroutine generatecavity_cpp +end subroutine pedra_driver diff --git a/src/tsless/CMakeLists.txt b/src/tsless/CMakeLists.txt new file mode 100644 index 000000000..2ddb4f428 --- /dev/null +++ b/src/tsless/CMakeLists.txt @@ -0,0 +1,11 @@ +# List of sources +list(APPEND sources_list tsless_cavity.F90 tsless_lapack.F90 tsless_precision.F90 tsless_symmetry.F90 tsless_weight_function.F90) + +# Update bar chart +if(BUILD_CHARTS) + update_bar_chart(${CMAKE_CURRENT_LIST_DIR} Fortran) +endif() + +foreach(_source ${sources_list}) + set_property(GLOBAL APPEND PROPERTY PCMSolver_Fortran_SOURCES ${CMAKE_CURRENT_LIST_DIR}/${_source}) +endforeach() diff --git a/src/tsless/tsless_cavity.F90 b/src/tsless/tsless_cavity.F90 new file mode 100644 index 000000000..4718c005b --- /dev/null +++ b/src/tsless/tsless_cavity.F90 @@ -0,0 +1,390 @@ +!pcmsolver_copyright_start +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2013 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: +!pcmsolver_copyright_end + +!> 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 tsless_precision + +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) + + use, intrinsic :: iso_c_binding + use tsless_symmetry + use tsless_weight_function + + !> 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, 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 + !> \param[in] nts + !> \param[in] x + !> \param[in] y + !> \param[in] z + !> \param[in] r + !> \param[in] xtscor + !> \param[in] ytscor + !> \param[in] ztscor + !> \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 000000000..fec813f01 --- /dev/null +++ b/src/tsless/tsless_lapack.F90 @@ -0,0 +1,158 @@ +!pcmsolver_copyright_start +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2013 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: +!pcmsolver_copyright_end + +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, k + + !> 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 000000000..f3fffd20a --- /dev/null +++ b/src/tsless/tsless_precision.F90 @@ -0,0 +1,42 @@ +!pcmsolver_copyright_start +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2013 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: +!pcmsolver_copyright_end + +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 000000000..bab45e202 --- /dev/null +++ b/src/tsless/tsless_symmetry.F90 @@ -0,0 +1,354 @@ +!pcmsolver_copyright_start +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2013 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: +!pcmsolver_copyright_end + +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 000000000..9d42f10d5 --- /dev/null +++ b/src/tsless/tsless_weight_function.F90 @@ -0,0 +1,251 @@ +!pcmsolver_copyright_start +! PCMSolver, an API for the Polarizable Continuum Model +! Copyright (C) 2013 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: +!pcmsolver_copyright_end + +module tsless_weight_function + +use, intrinsic :: iso_c_binding +use tsless_precision + +implicit none + +!> \brief Describes a weight function +!> Data fiels have dimensiont 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, 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 + 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, 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] printer logical unit for printing + !> ifun = 0 non-normalized function (gamma) + !> ifun = 1 normalized function (omega) + !> d0: threshold value + subroutine compute(wfun, ifun, 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 + 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 + !> Parameters + real(kind=dp), parameter :: d0 = 0.0_dp + + 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 + !> Allocations and clean-up + allocate(C(m, m)); C = 0.0_dp + allocate(b(m)); b = 0.0_dp + + ! Solve C * wfun%weight_0 = b + ! 0. Fill coefficient matrix + ! 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) * d0**(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 - d0**(k+1)) / (1.0_dp * (k+1)) + end do + end if + ! 1. Fill b vector (RHS) + b(1) = 1.0_dp + if (ifun .eq. 1) b(m) = 1.0_dp + ! 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/src/utils/Citation.hpp b/src/utils/Citation.hpp index 03f5bbd83..d515f57b4 100644 --- a/src/utils/Citation.hpp +++ b/src/utils/Citation.hpp @@ -52,8 +52,8 @@ inline std::string citation_message() rest << " With contributions from:" << std::endl; rest << " R. Bast (CMake framework)" << std::endl; rest << " U. Ekstroem (automatic differentiation library)" << std::endl; - rest << " J. Juselius (input parsing library and CMake framework)" << - std::endl; + rest << " J. Juselius (input parsing library and CMake framework)" << std::endl; + rest << " C. S. Pomelli (TsLess cavity library)" << std::endl; rest << " Theory: - J. Tomasi, B. Mennucci and R. Cammi:" << std::endl; rest << " \"Quantum Mechanical Continuum Solvation Models\", Chem. Rev., 105 (2005) 2999" << std::endl; diff --git a/src/utils/Molecule.cpp b/src/utils/Molecule.cpp index f854e367d..d3bf8e294 100644 --- a/src/utils/Molecule.cpp +++ b/src/utils/Molecule.cpp @@ -36,6 +36,8 @@ #include #include +#include + #include "Atom.hpp" #include "Element.hpp" #include "MathUtils.hpp" @@ -237,20 +239,16 @@ Molecule& Molecule::operator=(const Molecule& other) std::ostream & operator<<(std::ostream &os, const Molecule &m) { - // Declare formatting of Eigen output. - std::string sep = " "; - Eigen::IOFormat CleanFmt(Eigen::FullPrecision, Eigen::DontAlignCols, sep, "\n", "", - ""); - os << "Rotor type: " << rotorTypeList[m.rotor_] << std::endl; if (m.nAtoms_ != 0) { - os << " Center X Y Z " << - std::endl; - os << " ------------ ----------------- ----------------- -----------------" << - std::endl; + os << " Center X Y Z " << std::endl; + os << " ------------ ----------------- ----------------- -----------------" << std::endl; for (size_t i = 0; i < m.nAtoms_; ++i) { - os << std::setw(10) << m.atoms_[i].atomSymbol() << std::setw(15) <. - * + * * For information on the complete list of contributors to the * PCMSolver API, see: */ @@ -39,6 +39,10 @@ #include "Sphere.hpp" #include "Symmetry.hpp" +/*! Returns the hydrogen fluoride molecule + */ +inline Molecule HF(); + /*! Returns the ammonia molecule */ inline Molecule NH3(); @@ -126,6 +130,36 @@ Molecule dummy(double radius, const Eigen::Vector3d & center) 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/CMakeLists.txt b/tests/CMakeLists.txt index ef25f673a..8c8d851a3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -4,6 +4,7 @@ add_subdirectory(surface_function) add_subdirectory(input) add_subdirectory(numerical_quadrature) add_subdirectory(gepol) +add_subdirectory(tsless) add_subdirectory(dielectric_profile) add_subdirectory(green) add_subdirectory(bi_operators) diff --git a/tests/cpcm/CMakeLists.txt b/tests/cpcm/CMakeLists.txt index 334d3a85b..e8c5627c5 100644 --- a/tests/cpcm/CMakeLists.txt +++ b/tests/cpcm/CMakeLists.txt @@ -10,3 +10,7 @@ add_test(NAME cpcm_gepol-NH3_from-file COMMAND unit_tests.x [cpcm_gepol-NH3_from set_tests_properties(cpcm_gepol-NH3_from-file PROPERTIES DEPENDS cpcm_gepol-NH3) set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/cpcm_symmetry.cpp) add_test(NAME cpcm_symmetry COMMAND unit_tests.x [cpcm_symmetry]) +set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/cpcm_tsless-point.cpp) +add_test(NAME cpcm_tsless-point COMMAND unit_tests.x [cpcm_tsless-point]) +set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/cpcm_tsless-NH3.cpp) +add_test(NAME cpcm_tsless-NH3 COMMAND unit_tests.x [cpcm_tsless-NH3]) diff --git a/tests/cpcm/cpcm_gepol-NH3.cpp b/tests/cpcm/cpcm_gepol-NH3.cpp index 4c2c390e3..69a6ceef5 100644 --- a/tests/cpcm/cpcm_gepol-NH3.cpp +++ b/tests/cpcm/cpcm_gepol-NH3.cpp @@ -2,22 +2,22 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ @@ -25,8 +25,6 @@ #include "catch.hpp" -#include - #include "Config.hpp" #include diff --git a/tests/cpcm/cpcm_gepol-point.cpp b/tests/cpcm/cpcm_gepol-point.cpp index c0d09ccef..2e39c0c66 100644 --- a/tests/cpcm/cpcm_gepol-point.cpp +++ b/tests/cpcm/cpcm_gepol-point.cpp @@ -2,22 +2,22 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ diff --git a/tests/cpcm/cpcm_tsless-NH3.cpp b/tests/cpcm/cpcm_tsless-NH3.cpp index 5e234b95e..3ed5a2d3a 100644 --- a/tests/cpcm/cpcm_tsless-NH3.cpp +++ b/tests/cpcm/cpcm_tsless-NH3.cpp @@ -2,72 +2,60 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ /* pcmsolver_copyright_end */ -#define BOOST_TEST_MODULE CPCMSolverNH3TsLess - -#include -#include - -#include +#include "catch.hpp" #include "Config.hpp" #include #include "CollocationIntegrator.hpp" +#include "CPCMSolver.hpp" #include "DerivativeTypes.hpp" +#include "Molecule.hpp" +#include "TestingMolecules.hpp" #include "TsLessCavity.hpp" -#include "Vacuum.hpp" #include "UniformDielectric.hpp" -#include "CPCMSolver.hpp" +#include "Vacuum.hpp" /*! \class CPCMSolver * \test \b NH3TsLess tests CPCMSolver using ammonia and a TsLess cavity */ -BOOST_AUTO_TEST_CASE(NH3TsLess) +TEST_CASE("Test solver for the C-PCM with NH3 molecule and a TsLess cavity", "[solver][cpcm][cpcm_tsless-NH3]") { - // Set up cavity - Eigen::Vector3d N( -0.000000000, -0.104038047, 0.000000000); - Eigen::Vector3d H1(-0.901584415, 0.481847022, -1.561590016); - Eigen::Vector3d H2(-0.901584415, 0.481847022, 1.561590016); - Eigen::Vector3d H3( 1.803168833, 0.481847022, 0.000000000); - std::vector spheres; - Sphere sph1(N, 2.929075493); - Sphere sph2(H1, 2.267671349); - Sphere sph3(H2, 2.267671349); - Sphere sph4(H3, 2.267671349); - spheres.push_back(sph1); - spheres.push_back(sph2); - spheres.push_back(sph3); - spheres.push_back(sph4); - double area = 0.4; - TsLessCavity cavity(spheres, area); + 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); - CollocationIntegrator * diag = new CollocationIntegrator(); double permittivity = 78.39; - Vacuum gfInside = Vacuum(diag); - UniformDielectric gfOutside = - UniformDielectric(permittivity, diag); + Vacuum gfInside = Vacuum(); + UniformDielectric gfOutside = + UniformDielectric(permittivity); bool symm = true; double correction = 0.0; CPCMSolver solver(symm, correction); @@ -76,22 +64,14 @@ BOOST_AUTO_TEST_CASE(NH3TsLess) double Ncharge = 7.0; double Hcharge = 1.0; size_t size = cavity.size(); - Eigen::VectorXd fake_mep = computeMEP(molecule, cavity.elements()); - for (size_t i = 0; i < size; ++i) { - Eigen::Vector3d center = cavity.elementCenter(i); - double Ndistance = (center - N).norm(); - double H1distance = (center - H1).norm(); - double H2distance = (center - H2).norm(); - double H3distance = (center - H3).norm(); - fake_mep(i) = Ncharge / Ndistance + Hcharge / H1distance + Hcharge / H2distance + - Hcharge / H3distance; - } + 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} + // 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; + double totalASC = - (Ncharge + 3.0 * Hcharge) * (permittivity - 1) / + (permittivity + correction); double totalFakeASC = fake_asc.sum(); - std::cout << "totalASC - totalFakeASC = " << totalASC - totalFakeASC << std::endl; - BOOST_REQUIRE_CLOSE(totalASC, totalFakeASC, 4e-02); + 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 index 8b66d3a75..28e2e2760 100644 --- a/tests/cpcm/cpcm_tsless-point.cpp +++ b/tests/cpcm/cpcm_tsless-point.cpp @@ -2,82 +2,135 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ /* pcmsolver_copyright_end */ -#define BOOST_TEST_MODULE CPCMSolverpointChargeTsLess - -#include -#include - -#include +#include "catch.hpp" #include "Config.hpp" #include #include "CollocationIntegrator.hpp" +#include "CPCMSolver.hpp" #include "DerivativeTypes.hpp" -#include "TsLessCavity.hpp" -#include "Vacuum.hpp" #include "UniformDielectric.hpp" -#include "CPCMSolver.hpp" +#include "Vacuum.hpp" +#include "TestingMolecules.hpp" +#include "TsLessCavity.hpp" -/*! \class CPCMSolver - * \test \b pointChargeTsLess tests CPCMSolver using a point charge with a TsLess cavity - */ -BOOST_AUTO_TEST_CASE(pointChargeTsLess) +SCENARIO("Test solver for the C-PCM for a point charge and a TsLess cavity", "[solver][cpcm][cpcm_tsless-point]") { - // Set up cavity - Eigen::Vector3d N(0.0, 0.0, 0.0); - std::vector spheres; - Sphere sph1(N, 2.929075493); - spheres.push_back(sph1); - double area = 0.4; - TsLessCavity cavity(spheres, area); - - CollocationIntegrator * diag = new CollocationIntegrator(); - double permittivity = 78.39; - Vacuum * gfInside = new Vacuum(diag); - UniformDielectric * gfOutside = new - UniformDielectric(permittivity, diag); - bool symm = true; - double correction = 0.0; - CPCMSolver solver(symm, correction); - solver.buildSystemMatrix(cavity, gfInside, gfOutside); - - double charge = 8.0; - size_t size = cavity.size(); - Eigen::VectorXd fake_mep = computeMEP(molecule, cavity.elements()); - for (size_t i = 0; i < size; ++i) { - Eigen::Vector3d center = cavity.elementCenter(i); - double distance = center.norm(); - fake_mep(i) = charge / distance; + GIVEN("An isotropic environment modelled as a perfect conductor and a point charge") + { + double permittivity = 78.39; + Vacuum gfInside = Vacuum(); + UniformDielectric gfOutside = + UniformDielectric(permittivity); + bool symm = true; + double correction = 0.0; + + 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); + + 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); + + 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)); + } + } } - // The total ASC for a conductor is -Q - // for CPCM it will be -Q*[(epsilon-1)/epsilon} - Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); - fake_asc = solver.computeCharge(fake_mep); - double totalASC = - charge * (permittivity - 1) / permittivity; - double totalFakeASC = fake_asc.sum(); - std::cout << "totalASC - totalFakeASC = " << totalASC - totalFakeASC << std::endl; - BOOST_REQUIRE_CLOSE(totalASC, totalFakeASC, 4e-02); } diff --git a/tests/gepol/gepol_HF.cpp b/tests/gepol/gepol_HF.cpp new file mode 100644 index 000000000..fbae2c309 --- /dev/null +++ b/tests/gepol/gepol_HF.cpp @@ -0,0 +1,80 @@ +/* pcmsolver_copyright_start */ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2013 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: + */ +/* pcmsolver_copyright_end */ + +#define BOOST_TEST_MODULE GePolCavityNH3 + +#include +#include + +#include "Config.hpp" + +#include + +#include "GePolCavity.hpp" +#include "LoggerInterface.hpp" +#include "Molecule.hpp" +#include "PhysicalConstants.hpp" +#include "TestingMolecules.hpp" +#include "Symmetry.hpp" + +struct GePolCavityNH3Test { +protected: + GePolCavity cavity; + GePolCavityNH3Test() { SetUp(); } + void SetUp() { + Molecule molec = HF(); + double area = 0.02 / convertBohr2ToAngstrom2; + double probeRadius = 0.0 / convertBohrToAngstrom; + double minRadius = 0.2 / convertBohrToAngstrom; + cavity = GePolCavity(molec, area, probeRadius, minRadius); + } +}; + +BOOST_FIXTURE_TEST_CASE(size, GePolCavityNH3Test) +{ + int ref_size = 1688; + int size = cavity.size(); + BOOST_REQUIRE_EQUAL(size, ref_size); +} + +BOOST_FIXTURE_TEST_CASE(area, GePolCavityNH3Test) +{ + double ref_area = 110.64517236179323; + double area = cavity.elementArea().sum(); + BOOST_REQUIRE_CLOSE(area, ref_area, 1.0e-09); +} + +BOOST_FIXTURE_TEST_CASE(volume, GePolCavityNH3Test) +{ + double ref_volume = 105.98002389543836; + Eigen::Matrix3Xd elementCenter = cavity.elementCenter(); + Eigen::Matrix3Xd elementNormal = cavity.elementNormal(); + double volume = 0; + for ( int i = 0; i < cavity.size(); ++i ) { + volume += cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i)); + } + volume /= 3; + BOOST_REQUIRE_CLOSE(volume, ref_volume, 1.0e-09); +} diff --git a/tests/iefpcm/CMakeLists.txt b/tests/iefpcm/CMakeLists.txt index 88c1a3589..e0a9fd490 100644 --- a/tests/iefpcm/CMakeLists.txt +++ b/tests/iefpcm/CMakeLists.txt @@ -16,3 +16,7 @@ set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/iefpcm add_test(NAME iefpcm_symmetry COMMAND unit_tests.x [iefpcm_symmetry]) set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/iefpcm_anisotropic-symmetry.cpp) add_test(NAME iefpcm_anisotropic-symmetry COMMAND unit_tests.x [iefpcm_anisotropic-symmetry]) +set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/iefpcm_tsless-point.cpp) +add_test(NAME iefpcm_tsless-point COMMAND unit_tests.x [iefpcm_tsless-point]) +set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/iefpcm_tsless-NH3.cpp) +add_test(NAME iefpcm_tsless-NH3 COMMAND unit_tests.x [iefpcm_tsless-NH3]) diff --git a/tests/iefpcm/iefpcm_gepol-point.cpp b/tests/iefpcm/iefpcm_gepol-point.cpp index e98bb28dc..af56cbe7f 100644 --- a/tests/iefpcm/iefpcm_gepol-point.cpp +++ b/tests/iefpcm/iefpcm_gepol-point.cpp @@ -2,22 +2,22 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ @@ -25,8 +25,6 @@ #include "catch.hpp" -#include - #include "Config.hpp" #include diff --git a/tests/iefpcm/iefpcm_tsless-NH3.cpp b/tests/iefpcm/iefpcm_tsless-NH3.cpp index 6840610ed..166a0941b 100644 --- a/tests/iefpcm/iefpcm_tsless-NH3.cpp +++ b/tests/iefpcm/iefpcm_tsless-NH3.cpp @@ -2,33 +2,28 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ /* pcmsolver_copyright_end */ -#define BOOST_TEST_MODULE IEFSolverNH3TsLess - -#include -#include - -#include +#include "catch.hpp" #include "Config.hpp" @@ -37,31 +32,26 @@ #include "CollocationIntegrator.hpp" #include "DerivativeTypes.hpp" #include "TsLessCavity.hpp" +#include "Molecule.hpp" #include "Vacuum.hpp" +#include "TestingMolecules.hpp" #include "UniformDielectric.hpp" #include "IEFSolver.hpp" +#include "Symmetry.hpp" /*! \class IEFSolver * \test \b NH3TsLess tests IEFSolver using ammonia and a TsLess cavity */ -BOOST_AUTO_TEST_CASE(NH3TsLess) +TEST_CASE("Test solver for the IEFPCM with NH3 molecule and a TsLess cavity", "[solver][iefpcm][iefpcm_tsless-NH3]") { - // Set up cavity - Eigen::Vector3d N( -0.000000000, -0.104038047, 0.000000000); - Eigen::Vector3d H1(-0.901584415, 0.481847022, -1.561590016); - Eigen::Vector3d H2(-0.901584415, 0.481847022, 1.561590016); - Eigen::Vector3d H3( 1.803168833, 0.481847022, 0.000000000); - std::vector spheres; - Sphere sph1(N, 2.929075493); - Sphere sph2(H1, 2.267671349); - Sphere sph3(H2, 2.267671349); - Sphere sph4(H3, 2.267671349); - spheres.push_back(sph1); - spheres.push_back(sph2); - spheres.push_back(sph3); - spheres.push_back(sph4); - double area = 0.4; - TsLessCavity cavity(spheres, area); + 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 = Vacuum(); @@ -74,21 +64,12 @@ BOOST_AUTO_TEST_CASE(NH3TsLess) double Ncharge = 7.0; double Hcharge = 1.0; size_t size = cavity.size(); - Eigen::VectorXd fake_mep = Eigen::VectorXd::Zero(size); - for (size_t i = 0; i < size; ++i) { - Eigen::Vector3d center = cavity.elementCenter(i); - double Ndistance = (center - N).norm(); - double H1distance = (center - H1).norm(); - double H2distance = (center - H2).norm(); - double H3distance = (center - H3).norm(); - fake_mep(i) = Ncharge / Ndistance + Hcharge / H1distance + Hcharge / H2distance + - Hcharge / H3distance; - } + 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(); - std::cout << "totalASC - totalFakeASC = " << totalASC - totalFakeASC << std::endl; - BOOST_REQUIRE_CLOSE(totalASC, totalFakeASC, 4e-02); + 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 index 6da1c970b..d909672c7 100644 --- a/tests/iefpcm/iefpcm_tsless-point.cpp +++ b/tests/iefpcm/iefpcm_tsless-point.cpp @@ -2,33 +2,28 @@ /* * PCMSolver, an API for the Polarizable Continuum Model * Copyright (C) 2013-2015 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: */ /* pcmsolver_copyright_end */ -#define BOOST_TEST_MODULE IEFSolverpointChargeTsLess - -#include -#include - -#include +#include "catch.hpp" #include "Config.hpp" @@ -36,45 +31,105 @@ #include "CollocationIntegrator.hpp" #include "DerivativeTypes.hpp" -#include "TsLessCavity.hpp" #include "Vacuum.hpp" #include "UniformDielectric.hpp" #include "IEFSolver.hpp" +#include "TestingMolecules.hpp" +#include "TsLessCavity.hpp" -/*! \class IEFSolver - * \test \b pointChargeTsLess tests IEFSolver using a point charge with a TsLess cavity - */ -BOOST_AUTO_TEST_CASE(pointChargeTsLess) +SCENARIO("Test solver for the IEFPCM for a point charge and a TsLess cavity", "[solver][iefpcm][iefpcm_tsless-point]") { - // Set up cavity - Eigen::Vector3d N(0.0, 0.0, 0.0); - std::vector spheres; - Sphere sph1(N, 2.929075493); - spheres.push_back(sph1); - double area = 0.4; - TsLessCavity cavity(spheres, area); - - double permittivity = 78.39; - Vacuum gfInside = Vacuum(); - UniformDielectric gfOutside = - UniformDielectric(permittivity); - bool symm = true; - IEFSolver solver(symm); - solver.buildSystemMatrix(cavity, gfInside, gfOutside); - - double charge = 8.0; - size_t size = cavity.size(); - Eigen::VectorXd fake_mep = Eigen::VectorXd::Zero(size); - for (size_t i = 0; i < size; ++i) { - Eigen::Vector3d center = cavity.elementCenter(i); - double distance = center.norm(); - fake_mep(i) = charge / distance; + GIVEN("An isotropic environment and a point charge") + { + double permittivity = 78.39; + Vacuum gfInside = Vacuum(); + UniformDielectric gfOutside = + UniformDielectric(permittivity); + bool symm = true; + + 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); + + 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); + + 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)); + } + } } - // 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 = - charge * (permittivity - 1) / permittivity; - double totalFakeASC = fake_asc.sum(); - std::cout << "totalASC - totalFakeASC = " << totalASC - totalFakeASC << std::endl; - BOOST_REQUIRE_CLOSE(totalASC, totalFakeASC, 4e-02); } diff --git a/tests/tsless/CMakeLists.txt b/tests/tsless/CMakeLists.txt new file mode 100644 index 000000000..3d5b4d61a --- /dev/null +++ b/tests/tsless/CMakeLists.txt @@ -0,0 +1,6 @@ +set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/tsless_point.cpp) +add_test(NAME tsless_point COMMAND unit_tests.x [tsless_point]) +set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/tsless_NH3.cpp) +add_test(NAME tsless_NH3 COMMAND unit_tests.x [tsless_NH3]) +set_property(GLOBAL APPEND PROPERTY TestSources ${CMAKE_CURRENT_LIST_DIR}/tsless_HF.cpp) +add_test(NAME tsless_HF COMMAND unit_tests.x [tsless_HF]) diff --git a/tests/tsless/tsless_HF.cpp b/tests/tsless/tsless_HF.cpp new file mode 100644 index 000000000..11842997e --- /dev/null +++ b/tests/tsless/tsless_HF.cpp @@ -0,0 +1,82 @@ +/* pcmsolver_copyright_start */ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2013-2015 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: + */ +/* pcmsolver_copyright_end */ + +#include "catch.hpp" + +#include "Config.hpp" + +#include + +#include "Molecule.hpp" +#include "PhysicalConstants.hpp" +#include "TestingMolecules.hpp" +#include "TsLessCavity.hpp" + +TEST_CASE("TsLess cavity for an hydrogen fluoride molecule", "[tsless][tsless_HF]") +{ + Molecule molec = HF(); + double area = 0.02 / convertBohr2ToAngstrom2; + double probeRadius = 0.0 / convertBohrToAngstrom; + double minRadius = 0.2 / convertBohrToAngstrom; + 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 000000000..5c51cf2f4 --- /dev/null +++ b/tests/tsless/tsless_NH3.cpp @@ -0,0 +1,83 @@ +/* pcmsolver_copyright_start */ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2013-2015 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: + */ +/* pcmsolver_copyright_end */ + +#include "catch.hpp" + +#include "Config.hpp" + +#include + +#include "TsLessCavity.hpp" +#include "Molecule.hpp" +#include "PhysicalConstants.hpp" +#include "TestingMolecules.hpp" + +TEST_CASE("TsLess cavity for an ammonia molecule", "[tsless][tsless_NH3]") +{ + Molecule molec = NH3(); + double area = 0.02 / convertBohr2ToAngstrom2; + double probeRadius = 0.0 / convertBohrToAngstrom; + double minRadius = 0.2 / convertBohrToAngstrom; + 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 000000000..fe3e87963 --- /dev/null +++ b/tests/tsless/tsless_point.cpp @@ -0,0 +1,79 @@ +/* pcmsolver_copyright_start */ +/* + * PCMSolver, an API for the Polarizable Continuum Model + * Copyright (C) 2013-2015 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: + */ +/* pcmsolver_copyright_end */ + +#include "catch.hpp" + +#include "Config.hpp" + +#include + +#include "TsLessCavity.hpp" +#include "TestingMolecules.hpp" + +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)); + } +} diff --git a/tools/pcmsolver.py.in b/tools/pcmsolver.py.in index 1c6c57e33..8a238416e 100644 --- a/tools/pcmsolver.py.in +++ b/tools/pcmsolver.py.in @@ -665,12 +665,13 @@ def verify_geometry(keyword): print('Empty or incoherent geometry') sys.exit(1) # Convert only geometry to AU. Leave charges untouched - j = 0 - for i in range(length/4): - data[j] /= CODATAdict[CODATAyear] - data[j+1] /= CODATAdict[CODATAyear] - data[j+2] /= CODATAdict[CODATAyear] - j += 4 + if (isAngstrom): + j = 0 + for i in range(length/4): + data[j] /= CODATAdict[CODATAyear] + data[j+1] /= CODATAdict[CODATAyear] + data[j+2] /= CODATAdict[CODATAyear] + j += 4 def verify_spheres(keyword): length=len(keyword.get()) From c824fb0e13310f1d667294ddaff6b79ae64d7736 Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Tue, 19 Jul 2016 11:35:34 +0200 Subject: [PATCH 2/6] Remove spurious files --- src/utils/Molecule.cpp.orig | 287 ----------------------- tests/cpcm/cpcm_gepol-NH3.cpp.orig | 81 ------- tests/iefpcm/iefpcm_gepol-point.cpp.orig | 139 ----------- 3 files changed, 507 deletions(-) delete mode 100644 src/utils/Molecule.cpp.orig delete mode 100644 tests/cpcm/cpcm_gepol-NH3.cpp.orig delete mode 100644 tests/iefpcm/iefpcm_gepol-point.cpp.orig diff --git a/src/utils/Molecule.cpp.orig b/src/utils/Molecule.cpp.orig deleted file mode 100644 index 26a0dcbf0..000000000 --- a/src/utils/Molecule.cpp.orig +++ /dev/null @@ -1,287 +0,0 @@ -/* pcmsolver_copyright_start */ -/* - * PCMSolver, an API for the Polarizable Continuum Model -<<<<<<< HEAD - * Copyright (C) 2013-2015 Roberto Di Remigio, Luca Frediani and contributors -======= - * Copyright (C) 2013-2016 Roberto Di Remigio, Luca Frediani and contributors ->>>>>>> master - * - * 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: - */ -/* pcmsolver_copyright_end */ - -#include "Molecule.hpp" - -#include -#include -#include -#include -#include -#include - -#include "Config.hpp" - -#include - -#include -#include - -#include "Atom.hpp" -#include "cavity/Element.hpp" -#include "MathUtils.hpp" -#include "Symmetry.hpp" - -Molecule::Molecule(int nat, const Eigen::VectorXd & chg, const Eigen::VectorXd & m, - const Eigen::Matrix3Xd & geo, const std::vector & at, - const std::vector & sph) - : nAtoms_(nat), charges_(chg), masses_(m), geometry_(geo), atoms_(at), spheres_(sph) -{ - rotor_ = findRotorType(); - pointGroup_ = buildGroup(0, 0, 0, 0); -} - -Molecule::Molecule(int nat, const Eigen::VectorXd & chg, const Eigen::VectorXd & m, - const Eigen::Matrix3Xd & geo, const std::vector & at, - const std::vector & sph, - int nr_gen, int gen[3]) - : nAtoms_(nat), charges_(chg), masses_(m), geometry_(geo), atoms_(at), spheres_(sph) -{ - rotor_ = findRotorType(); - pointGroup_ = buildGroup(nr_gen, gen[0], gen[1], gen[2]); -} - -Molecule::Molecule(int nat, const Eigen::VectorXd & chg, const Eigen::VectorXd & m, - const Eigen::Matrix3Xd & geo, const std::vector & at, - const std::vector & sph, - const Symmetry & pg) - : nAtoms_(nat), charges_(chg), masses_(m), geometry_(geo), atoms_(at), spheres_(sph), - pointGroup_(pg) -{ - rotor_ = findRotorType(); -} - -Molecule::Molecule(const std::vector & sph) - : nAtoms_(sph.size()), spheres_(sph) -{ - charges_ = Eigen::VectorXd::Ones(nAtoms_); - masses_.resize(nAtoms_); - geometry_.resize(Eigen::NoChange, nAtoms_); - for (size_t i = 0; i < nAtoms_; ++i) { - masses_(i) = spheres_[i].radius; - geometry_.col(i) = spheres_[i].center; - double charge = charges_(i); - double mass = masses_(i); - atoms_.push_back( Atom("Dummy", "Du", charge, mass, mass, geometry_.col(i)) ); - } - rotor_ = findRotorType(); - pointGroup_ = buildGroup(0, 0, 0, 0); -} - -Molecule::Molecule(const Molecule &other) -{ - *this = other; -} - -Eigen::Vector3d Molecule::centerOfMass() -{ - Eigen::Vector3d com; - com << 0.0, 0.0, 0.0; - for (size_t i = 0; i < nAtoms_; ++i) { - com += masses_(i) * atoms_[i].position; - } - com *= 1.0/masses_.sum(); - return com; -} - -Eigen::Matrix3d Molecule::inertiaTensor() -{ - Eigen::Matrix3d inertia = Eigen::Matrix3d::Zero(); - - for (size_t i = 0; i < nAtoms_; ++i) { - // Diagonal - inertia(0,0) += masses_(i) * (geometry_(1,i) * geometry_(1,i) + geometry_(2, - i) * geometry_(2,i)); - inertia(1,1) += masses_(i) * (geometry_(0,i) * geometry_(0,i) + geometry_(2, - i) * geometry_(2,i)); - inertia(2,2) += masses_(i) * (geometry_(0,i) * geometry_(0,i) + geometry_(1, - i) * geometry_(1,i)); - - // Off-diagonal - inertia(0,1) -= masses_(i) * (geometry_(0,i) * geometry_(1,i)); - inertia(0,2) -= masses_(i) * (geometry_(0,i) * geometry_(2,i)); - inertia(1,2) -= masses_(i) * (geometry_(1,i) * geometry_(2,i)); - } - // Now symmetrize - hermitivitize(inertia); - - // Check elements for a numerical zero and make it a hard zero - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - if (fabs(inertia(i,j)) < 1.0e-14) { - inertia(i,j) = 0.0; - } - } - } - - return inertia; -} - -rotorType Molecule::findRotorType() -{ - rotorType type; - if (nAtoms_ == 1) { - type = rtAtom; - } else { - // Get inertia tensor - Eigen::Matrix3d inertia = inertiaTensor(); - // Diagonalize inertia tensor V^t * I * V - Eigen::SelfAdjointEigenSolver eigenSolver(inertia); - if (eigenSolver.info() != Eigen::Success) abort(); - // Determine the degeneracy of the eigenvalues. - int deg = 0; - double tmp, abs, rel; - for (int i = 0; i < 2; ++i) { - for (int j = i + 1; j < 3 && deg < 2; ++j) { // Check i and j != i - abs = fabs(eigenSolver.eigenvalues()[i] - eigenSolver.eigenvalues()[j]); - tmp = eigenSolver.eigenvalues()[j]; // Because the eigenvalues are already in ascending order. - if (abs > 1.0e-14) { - rel = abs/tmp; - } else { - rel = 0.0; - } - if (rel < 1.0e-8) { - ++deg; - } - } - } - // Get the rotor type based on the degeneracy. - if (eigenSolver.eigenvalues()[0] == 0.0) { - type = rtLinear; - } else if (deg == 2) { - type = rtSpherical; - } else if (deg == 1) { // We do not distinguish between prolate and oblate. - type = rtSymmetric; - } else { - type = rtAsymmetric; - } - } - - return type; -} - -void Molecule::translate(const Eigen::Vector3d &translationVector) -{ - // Translate the geometry_ matrix and update the geometric data in atoms_. - for (size_t i = 0; i < nAtoms_; ++i) { - geometry_.col(i) -= translationVector; - Eigen::Vector3d tmp = geometry_.col(i); - atoms_[i].position = tmp; - } -} - -void Molecule::moveToCOM() -{ - Eigen::Vector3d com = centerOfMass(); - this->translate(com); -} - -void Molecule::rotate(const Eigen::Matrix3d &rotationMatrix) -{ - // Rotate the geometry_ matrix and update the geometric data in atoms_. - geometry_ *= - rotationMatrix; // The power of Eigen: geometry_ = geometry_ * rotationMatrix; - for (size_t i = 0; i < nAtoms_; ++i) { - Eigen::Vector3d tmp = geometry_.col(i); - atoms_[i].position = tmp; - } -} - -void Molecule::moveToPAF() -{ - Eigen::Matrix3d inertia = inertiaTensor(); - // Diagonalize inertia tensor V^t * I * V - Eigen::SelfAdjointEigenSolver eigenSolver(inertia); - if (eigenSolver.info() != Eigen::Success) abort(); - // Rotate to Principal Axes Frame - this->rotate(eigenSolver.eigenvectors()); - std::cout << eigenSolver.eigenvalues() << std::endl; -} - -Molecule& Molecule::operator=(const Molecule& other) -{ - // Self assignment is bad - if (this == &other) - return *this; - - nAtoms_ = other.nAtoms_; - charges_ = other.charges_; - masses_ = other.masses_; - geometry_ = other.geometry_; - atoms_ = other.atoms_; - spheres_ = other.spheres_; - rotor_ = other.rotor_; - pointGroup_ = other.pointGroup_; - - return *this; -} - -std::ostream & operator<<(std::ostream &os, const Molecule &m) -{ - if (m.nAtoms_ != 0) { - os << " Geometry (in Angstrom)" << std::endl; - os << " Center X Y Z \n"; - os << "------------ ------------ ------------ ------------\n"; - for (size_t i = 0; i < m.nAtoms_; ++i) { - os << boost::format("%|=12s|") % m.atoms_[i].symbol; - os << boost::format(" %10.6f ") % (m.geometry_(0, i) * bohrToAngstrom()); - os << boost::format(" %10.6f ") % (m.geometry_(1, i) * bohrToAngstrom()); - os << boost::format(" %10.6f ") % (m.geometry_(2, i) * bohrToAngstrom()); - os << std::endl; - } - os << "Rotor type: " << rotorTypeList[m.rotor_]; - } else { - os << " No atoms in this molecule!" << std::endl; - } - - return os; - -} - -Eigen::VectorXd computeMEP(const Molecule & mol, const std::vector & el) -{ - Eigen::VectorXd mep = Eigen::VectorXd::Zero(el.size()); - for (size_t i = 0; i < mol.nAtoms(); ++i) { - for (size_t j = 0; j < el.size(); ++j) { - double dist = (mol.geometry().col(i) - el[j].center()).norm(); - mep(j) += mol.charges(i) / dist; - } - } - return mep; -} - -Eigen::VectorXd computeMEP(const std::vector & el, double charge, const Eigen::Vector3d & origin) -{ - Eigen::VectorXd mep = Eigen::VectorXd::Zero(el.size()); - for (size_t i = 0; i < el.size(); ++i) { - double dist = (origin - el[i].center()).norm(); - mep(i) += charge / dist; - } - return mep; -} diff --git a/tests/cpcm/cpcm_gepol-NH3.cpp.orig b/tests/cpcm/cpcm_gepol-NH3.cpp.orig deleted file mode 100644 index c2b9a08c4..000000000 --- a/tests/cpcm/cpcm_gepol-NH3.cpp.orig +++ /dev/null @@ -1,81 +0,0 @@ -/* pcmsolver_copyright_start */ -/* - * PCMSolver, an API for the Polarizable Continuum Model -<<<<<<< HEAD - * Copyright (C) 2013-2015 Roberto Di Remigio, Luca Frediani and contributors -======= - * Copyright (C) 2013-2016 Roberto Di Remigio, Luca Frediani and contributors ->>>>>>> master - * - * 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: - */ -/* pcmsolver_copyright_end */ - -#include "catch.hpp" - -#include - - -#include - -#include "bi_operators/CollocationIntegrator.hpp" -#include "green/DerivativeTypes.hpp" -#include "cavity/GePolCavity.hpp" -#include "utils/Molecule.hpp" -#include "green/Vacuum.hpp" -#include "green/UniformDielectric.hpp" -#include "solver/CPCMSolver.hpp" -#include "utils/Symmetry.hpp" -#include "TestingMolecules.hpp" - -/*! \class CPCMSolver - * \test \b NH3GePol tests CPCMSolver using ammonia and a GePol cavity - */ -TEST_CASE("Test solver for the C-PCM with NH3 molecule and a GePol cavity", "[solver][cpcm][cpcm_gepol-NH3]") -{ - Molecule molecule = NH3(); - - double area = 0.4; - double probeRadius = 0.0; - double minRadius = 100.0; - GePolCavity cavity = GePolCavity(molecule, area, probeRadius, minRadius); - cavity.saveCavity("nh3.npz"); - - double permittivity = 78.39; - Vacuum<> gf_i; - UniformDielectric<> gf_o(permittivity); - bool symm = true; - double correction = 0.8; - CPCMSolver solver(symm, correction); - solver.buildSystemMatrix(cavity, gf_i, gf_o); - - double Ncharge = 7.0; - double Hcharge = 1.0; - int size = cavity.size(); - Eigen::VectorXd fake_mep = computeMEP(molecule, 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/iefpcm/iefpcm_gepol-point.cpp.orig b/tests/iefpcm/iefpcm_gepol-point.cpp.orig deleted file mode 100644 index 2838b55b8..000000000 --- a/tests/iefpcm/iefpcm_gepol-point.cpp.orig +++ /dev/null @@ -1,139 +0,0 @@ -/* pcmsolver_copyright_start */ -/* - * PCMSolver, an API for the Polarizable Continuum Model -<<<<<<< HEAD - * Copyright (C) 2013-2015 Roberto Di Remigio, Luca Frediani and contributors -======= - * Copyright (C) 2013-2016 Roberto Di Remigio, Luca Frediani and contributors ->>>>>>> master - * - * 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 -<<<<<<< HEAD - * PCMSolver API, see: -======= - * PCMSolver API, see: ->>>>>>> master - */ -/* pcmsolver_copyright_end */ - -#include "catch.hpp" - -#include - -#include "Config.hpp" - -#include - -#include "bi_operators/CollocationIntegrator.hpp" -#include "green/DerivativeTypes.hpp" -#include "cavity/GePolCavity.hpp" -#include "green/Vacuum.hpp" -#include "green/UniformDielectric.hpp" -#include "solver/IEFSolver.hpp" -#include "TestingMolecules.hpp" - -SCENARIO("Test solver for the IEFPCM for a point charge and a GePol cavity", "[solver][iefpcm][iefpcm_gepol-point]") -{ - GIVEN("An isotropic environment and a point charge") - { - double permittivity = 78.39; - Vacuum<> gf_i; - UniformDielectric<> gf_o(permittivity); - bool symm = true; - - double charge = 8.0; - double totalASC = - charge * (permittivity - 1) / permittivity; - - /*! \class IEFSolver - * \test \b pointChargeGePol tests IEFSolver using a point charge with a GePol 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; - GePolCavity cavity = GePolCavity(point, area, probeRadius, minRadius); - cavity.saveCavity("point.npz"); - - IEFSolver solver(symm); - solver.buildSystemMatrix(cavity, gf_i, gf_o); - - int size = cavity.size(); - Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge); - Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); - fake_asc = solver.computeCharge(fake_mep); - - for (int i = 0; i < size; ++i) { - INFO("fake_mep(" << i << ") = " << fake_mep(i)); - } - for (int 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 pointChargeShiftedGePol tests IEFSolver using a point charge with a GePol 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; - GePolCavity cavity = GePolCavity(point, area, probeRadius, minRadius); - - IEFSolver solver(symm); - solver.buildSystemMatrix(cavity, gf_i, gf_o); - - double charge = 8.0; - int size = cavity.size(); - Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge, origin); - Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); - fake_asc = solver.computeCharge(fake_mep); - - for (int i = 0; i < size; ++i) { - INFO("fake_mep(" << i << ") = " << fake_mep(i)); - } - for (int 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)); - } - } - } -} From c444a363775bb548a96b8cfac9a0578af40f339c Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Fri, 31 Mar 2017 10:51:47 +0200 Subject: [PATCH 3/6] Update Autocmake --- cmake/autocmake.yml | 1 + cmake/autocmake/generate.py | 16 +++++++++++----- cmake/downloaded/autocmake_safeguards.cmake | 3 ++- cmake/update.py | 10 ++++++++-- setup.py | 2 +- 5 files changed, 23 insertions(+), 9 deletions(-) diff --git a/cmake/autocmake.yml b/cmake/autocmake.yml index c3eee42ea..a2037eeba 100644 --- a/cmake/autocmake.yml +++ b/cmake/autocmake.yml @@ -2,6 +2,7 @@ name: PCMSolver min_cmake_version: 2.8.10 setup_script: setup.py url_root: https://github.com/coderefinery/autocmake/raw/master/ +default_build_type: debug modules: - compilers: diff --git a/cmake/autocmake/generate.py b/cmake/autocmake/generate.py index ee0c4229e..3c24eca7a 100644 --- a/cmake/autocmake/generate.py +++ b/cmake/autocmake/generate.py @@ -45,7 +45,7 @@ def autogenerated_notice(): return '\n'.join(s) -def gen_setup(config, relative_path, setup_script_name): +def gen_setup(config, default_build_type, relative_path, setup_script_name): """ Generate setup script. """ @@ -76,7 +76,7 @@ def gen_setup(config, relative_path, setup_script_name): rest = ' '.join(opt.split()[1:]).strip() options.append([first, rest]) - options.append(['--type=', 'Set the CMake build type (debug, release, or relwithdeb) [default: release].']) + options.append(['--type=', 'Set the CMake build type (debug, release, relwithdebinfo, minsizerel) [default: {0}].'.format(default_build_type)]) options.append(['--generator=', 'Set the CMake build system generator [default: Unix Makefiles].']) options.append(['--show', 'Show CMake command and exit.']) options.append(['--cmake-executable=', 'Set the CMake executable [default: cmake].']) @@ -118,7 +118,7 @@ def gen_setup(config, relative_path, setup_script_name): return s -def gen_cmakelists(project_name, min_cmake_version, relative_path, modules): +def gen_cmakelists(project_name, min_cmake_version, default_build_type, relative_path, modules): """ Generate CMakeLists.txt. """ @@ -137,9 +137,15 @@ def gen_cmakelists(project_name, min_cmake_version, relative_path, modules): s.append('\n# do not rebuild if rules (compiler flags) change') s.append('set(CMAKE_SKIP_RULE_DEPENDENCY TRUE)') - s.append('\n# if CMAKE_BUILD_TYPE undefined, we set it to Debug') + build_type_capitalized = {'debug': 'Debug', + 'release': 'Release', + 'relwithdebinfo': 'RelWithDebInfo', + 'minsizerel': 'MinSizeRel'} + + _build_type = build_type_capitalized[default_build_type] + s.append('\n# if CMAKE_BUILD_TYPE undefined, we set it to {0}'.format(_build_type)) s.append('if(NOT CMAKE_BUILD_TYPE)') - s.append(' set(CMAKE_BUILD_TYPE "Debug")') + s.append(' set(CMAKE_BUILD_TYPE "{0}")'.format(_build_type)) s.append('endif()') if len(modules) > 0: diff --git a/cmake/downloaded/autocmake_safeguards.cmake b/cmake/downloaded/autocmake_safeguards.cmake index 39505a7db..a61389661 100644 --- a/cmake/downloaded/autocmake_safeguards.cmake +++ b/cmake/downloaded/autocmake_safeguards.cmake @@ -17,6 +17,7 @@ string(TOUPPER "${CMAKE_BUILD_TYPE}" cmake_build_type_toupper) if(NOT cmake_build_type_tolower STREQUAL "debug" AND NOT cmake_build_type_tolower STREQUAL "release" AND + NOT cmake_build_type_tolower STREQUAL "minsizerel" AND NOT cmake_build_type_tolower STREQUAL "relwithdebinfo") - message(FATAL_ERROR "Unknown build type \"${CMAKE_BUILD_TYPE}\". Allowed values are Debug, Release, RelWithDebInfo (case-insensitive).") + message(FATAL_ERROR "Unknown build type \"${CMAKE_BUILD_TYPE}\". Allowed values are Debug, Release, RelWithDebInfo, and MinSizeRel (case-insensitive).") endif() diff --git a/cmake/update.py b/cmake/update.py index 9508eddc5..e0d2f1a70 100644 --- a/cmake/update.py +++ b/cmake/update.py @@ -147,6 +147,12 @@ def process_yaml(argv): sys.stderr.write("ERROR: you have to specify min_cmake_version in autocmake.yml\n") sys.exit(-1) + if 'default_build_type' in config: + default_build_type = config['default_build_type'].lower() + else: + sys.stderr.write("ERROR: you have to specify default_build_type in autocmake.yml\n") + sys.exit(-1) + if 'setup_script' in config: setup_script_name = config['setup_script'] else: @@ -173,13 +179,13 @@ def process_yaml(argv): # create CMakeLists.txt print('- generating CMakeLists.txt') - s = gen_cmakelists(project_name, min_cmake_version, relative_path, modules) + s = gen_cmakelists(project_name, min_cmake_version, default_build_type, relative_path, modules) with open(os.path.join(project_root, 'CMakeLists.txt'), 'w') as f: f.write('{0}\n'.format('\n'.join(s))) # create setup script print('- generating setup script') - s = gen_setup(cleaned_config, relative_path, setup_script_name) + s = gen_setup(cleaned_config, default_build_type, relative_path, setup_script_name) file_path = os.path.join(project_root, setup_script_name) with open(file_path, 'w') as f: f.write('{0}\n'.format('\n'.join(s))) diff --git a/setup.py b/setup.py index 7ed691b10..e64f9f55c 100755 --- a/setup.py +++ b/setup.py @@ -40,7 +40,7 @@ --build-boost= Deactivate Boost detection and build on-the-fly [default: OFF]. --static Create only the static library [default: False]. --eigen= Root directory for Eigen3 [default: '']. - --type= Set the CMake build type (debug, release, or relwithdeb) [default: release]. + --type= Set the CMake build type (debug, release, relwithdebinfo, minsizerel) [default: debug]. --generator= Set the CMake build system generator [default: Unix Makefiles]. --show Show CMake command and exit. --cmake-executable= Set the CMake executable [default: cmake]. From c392bff7e71315484ff527ff3341e0d992d226fd Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Sun, 2 Apr 2017 18:17:25 +0200 Subject: [PATCH 4/6] Fix usage of minimal distance parameter in TsLess sources --- src/tsless/CMakeLists.txt | 2 +- .../{tsless_cavity.F90 => tsless_cavity.f90} | 12 ++++---- .../{tsless_lapack.F90 => tsless_lapack.f90} | 12 ++++---- ...ess_precision.F90 => tsless_precision.f90} | 0 ...sless_symmetry.F90 => tsless_symmetry.f90} | 0 ...unction.F90 => tsless_weight_function.f90} | 28 +++++++++---------- 6 files changed, 27 insertions(+), 27 deletions(-) rename src/tsless/{tsless_cavity.F90 => tsless_cavity.f90} (99%) rename src/tsless/{tsless_lapack.F90 => tsless_lapack.f90} (98%) rename src/tsless/{tsless_precision.F90 => tsless_precision.f90} (100%) rename src/tsless/{tsless_symmetry.F90 => tsless_symmetry.f90} (100%) rename src/tsless/{tsless_weight_function.F90 => tsless_weight_function.f90} (92%) diff --git a/src/tsless/CMakeLists.txt b/src/tsless/CMakeLists.txt index e02d983f7..737444425 100644 --- a/src/tsless/CMakeLists.txt +++ b/src/tsless/CMakeLists.txt @@ -1,5 +1,5 @@ # List of sources -list(APPEND sources_list tsless_cavity.F90 tsless_lapack.F90 tsless_precision.F90 tsless_symmetry.F90 tsless_weight_function.F90) +list(APPEND sources_list tsless_cavity.f90 tsless_lapack.f90 tsless_precision.f90 tsless_symmetry.f90 tsless_weight_function.f90) add_library(tsless OBJECT ${sources_list}) set_target_properties(tsless PROPERTIES POSITION_INDEPENDENT_CODE 1 ) diff --git a/src/tsless/tsless_cavity.F90 b/src/tsless/tsless_cavity.f90 similarity index 99% rename from src/tsless/tsless_cavity.F90 rename to src/tsless/tsless_cavity.f90 index b9e139e24..f209004f8 100644 --- a/src/tsless/tsless_cavity.F90 +++ b/src/tsless/tsless_cavity.f90 @@ -1,22 +1,22 @@ !pcmsolver_copyright_start ! PCMSolver, an API for the Polarizable Continuum Model ! Copyright (C) 2013-2016 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: !pcmsolver_copyright_end @@ -129,7 +129,7 @@ subroutine tsless_driver(maxts, maxsph, maxvert, nesfp, nts, ntsirr, addsph, & ! generation and weighting of sampling points ! Generate weight function - weights = create_weight_function(nderiv=nord, normalization=ifun, printer=print_unit) + 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) diff --git a/src/tsless/tsless_lapack.F90 b/src/tsless/tsless_lapack.f90 similarity index 98% rename from src/tsless/tsless_lapack.F90 rename to src/tsless/tsless_lapack.f90 index 846a26d99..9047d8963 100644 --- a/src/tsless/tsless_lapack.F90 +++ b/src/tsless/tsless_lapack.f90 @@ -1,22 +1,22 @@ !pcmsolver_copyright_start ! PCMSolver, an API for the Polarizable Continuum Model ! Copyright (C) 2013-2016 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: !pcmsolver_copyright_end @@ -51,7 +51,7 @@ pure subroutine solve_linear_system(ndim, A, b, x) real(kind=dp), intent(out) :: x(ndim) !> Local variables real(kind=dp), allocatable :: Ainverse(:, :) - integer(kind=regint_k) :: i, j, k + integer(kind=regint_k) :: i, j !> Compute inverse by LU decomposition allocate(Ainverse(ndim, ndim)); Ainverse = 0.0_dp diff --git a/src/tsless/tsless_precision.F90 b/src/tsless/tsless_precision.f90 similarity index 100% rename from src/tsless/tsless_precision.F90 rename to src/tsless/tsless_precision.f90 diff --git a/src/tsless/tsless_symmetry.F90 b/src/tsless/tsless_symmetry.f90 similarity index 100% rename from src/tsless/tsless_symmetry.F90 rename to src/tsless/tsless_symmetry.f90 diff --git a/src/tsless/tsless_weight_function.F90 b/src/tsless/tsless_weight_function.f90 similarity index 92% rename from src/tsless/tsless_weight_function.F90 rename to src/tsless/tsless_weight_function.f90 index 714ebee70..6768fd72a 100644 --- a/src/tsless/tsless_weight_function.F90 +++ b/src/tsless/tsless_weight_function.f90 @@ -1,22 +1,22 @@ !pcmsolver_copyright_start ! PCMSolver, an API for the Polarizable Continuum Model ! Copyright (C) 2013-2016 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: !pcmsolver_copyright_end @@ -29,7 +29,7 @@ module tsless_weight_function implicit none !> \brief Describes a weight function -!> Data fiels have dimensiont 2*nr_derivative+3 +!> Data fiels have dimensions 2*nr_derivative+3 type, public :: weight_function !> Number of continuous derivatives integer(kind=regint_k) :: nr_derivative @@ -49,13 +49,14 @@ module tsless_weight_function contains - function create_weight_function(nderiv, normalization, printer) result(wfun) + 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 @@ -69,7 +70,7 @@ function create_weight_function(nderiv, normalization, printer) result(wfun) end if call allocate_weight_function(wfun, nderiv) - call compute(wfun, normalization, print_out) + call compute(wfun, normalization, dmin, print_out) end function create_weight_function @@ -106,11 +107,11 @@ end subroutine destroy_weight_function !> \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) - !> d0: threshold value - subroutine compute(wfun, ifun, printer) + subroutine compute(wfun, ifun, dmin, printer) use, intrinsic :: iso_fortran_env, only: output_unit @@ -119,14 +120,13 @@ subroutine compute(wfun, ifun, printer) !> 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 - !> Parameters - real(kind=dp), parameter :: d0 = 0.0_dp if (present(printer)) then print_out = printer @@ -163,13 +163,13 @@ subroutine compute(wfun, ifun, printer) 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) * d0**(k-j) + 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 - d0**(k+1)) / (1.0_dp * (k+1)) + C(m, k+1) = (1.0_dp - dmin**(k+1)) / (1.0_dp * (k+1)) end do end if ! 1. Fill b vector (RHS) From a310fc699fdbabb899f6412e8268dee51f9d6c52 Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Sun, 2 Apr 2017 21:54:40 +0200 Subject: [PATCH 5/6] Document input/output paramteres in leopardi_grid subroutine --- src/tsless/tsless_cavity.f90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/tsless/tsless_cavity.f90 b/src/tsless/tsless_cavity.f90 index f209004f8..59fb1585c 100644 --- a/src/tsless/tsless_cavity.f90 +++ b/src/tsless/tsless_cavity.f90 @@ -210,15 +210,15 @@ end subroutine tsless_driver !> \brief Create and translate Leopardi points !> \author Christian S. Pomelli - !> \param[in] ntssp - !> \param[in] nts - !> \param[in] x - !> \param[in] y - !> \param[in] z - !> \param[in] r - !> \param[in] xtscor - !> \param[in] ytscor - !> \param[in] ztscor + !> \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) From 19c974f95b5e7978df701a0462717fc83636f6fc Mon Sep 17 00:00:00 2001 From: Roberto Di Remigio Date: Wed, 5 Apr 2017 14:42:40 +0200 Subject: [PATCH 6/6] Clarify in-code documentation --- src/tsless/tsless_weight_function.f90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/tsless/tsless_weight_function.f90 b/src/tsless/tsless_weight_function.f90 index 6768fd72a..dd792b9bf 100644 --- a/src/tsless/tsless_weight_function.f90 +++ b/src/tsless/tsless_weight_function.f90 @@ -147,12 +147,8 @@ subroutine compute(wfun, ifun, dmin, printer) else m = 2*n + 3 end if - !> Allocations and clean-up + !> Allocate, clean up and fill coefficient matrix C allocate(C(m, m)); C = 0.0_dp - allocate(b(m)); b = 0.0_dp - - ! Solve C * wfun%weight_0 = b - ! 0. Fill coefficient matrix ! l = 1 boundary conditions do j = 0_regint_k, n do k = j, m - 1_regint_k @@ -172,9 +168,10 @@ subroutine compute(wfun, ifun, dmin, printer) C(m, k+1) = (1.0_dp - dmin**(k+1)) / (1.0_dp * (k+1)) end do end if - ! 1. Fill b vector (RHS) + ! 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 + 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