1 /* 2 * PCMSolver, an API for the Polarizable Continuum Model 3 * Copyright (C) 2020 Roberto Di Remigio, Luca Frediani and contributors. 4 * 5 * This file is part of PCMSolver. 6 * 7 * PCMSolver is free software: you can redistribute it and/or modify 8 * it under the terms of the GNU Lesser General Public License as published by 9 * the Free Software Foundation, either version 3 of the License, or 10 * (at your option) any later version. 11 * 12 * PCMSolver is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU Lesser General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General Public License 18 * along with PCMSolver. If not, see <http://www.gnu.org/licenses/>. 19 * 20 * For information on the complete list of contributors to the 21 * PCMSolver API, see: <http://pcmsolver.readthedocs.io/> 22 */ 23 24 #pragma once 25 26 #include <iosfwd> 27 #include <vector> 28 29 #include "Config.hpp" 30 31 #include <Eigen/Core> 32 33 #include "Element.hpp" 34 #include "utils/Molecule.hpp" 35 #include "utils/Sphere.hpp" 36 #include "utils/Symmetry.hpp" 37 38 /*! \file ICavity.hpp */ 39 40 namespace pcm { 41 using cavity::Element; 42 using utils::Sphere; 43 44 /*! \class ICavity 45 * \brief Abstract Base Class for cavities. 46 * \author Krzysztof Mozgawa 47 * \date 2011 48 * 49 * This class represents a cavity made of spheres, its surface being discretized in 50 * terms of finite elements. 51 */ 52 class ICavity { 53 protected: 54 /// List of spheres 55 std::vector<Sphere> spheres_; 56 /// The molecule to be wrapped by the cavity 57 Molecule molecule_; 58 /// Number of finite elements generated 59 PCMSolverIndex nElements_; 60 /// Number of irreducible finite elements 61 PCMSolverIndex nIrrElements_; 62 /// Whether the cavity has been built 63 bool built; 64 /// Coordinates of elements centers 65 Eigen::Matrix3Xd elementCenter_; 66 /// Outward-pointing normal vectors to the elements centers 67 Eigen::Matrix3Xd elementNormal_; 68 /// Elements areas 69 Eigen::VectorXd elementArea_; 70 /// Number of spheres 71 int nSpheres_; 72 /// Centers of the sphere the elements belong to 73 Eigen::Matrix3Xd elementSphereCenter_; 74 /// Radii of the sphere the elements belong to 75 Eigen::VectorXd elementRadius_; 76 /// Spheres centers 77 Eigen::Matrix3Xd sphereCenter_; 78 /// Spheres radii 79 Eigen::VectorXd sphereRadius_; 80 /// List of finite elements 81 std::vector<Element> elements_; 82 /// Molecular point group 83 Symmetry pointGroup_; 84 85 private: 86 /*! \brief Creates the cavity and discretizes its surface. 87 * 88 * Has to be implemented by classes lower down in the inheritance hierarchy 89 */ 90 virtual void makeCavity() = 0; 91 virtual std::ostream & printCavity(std::ostream & os) = 0; 92 93 public: 94 //! Default constructor 95 ICavity(); 96 /*! \brief Constructor from a single sphere 97 * \param[in] sph the sphere 98 * 99 * Only used when we have to deal with a single sphere, i.e. in the unit tests 100 */ 101 ICavity(const Sphere & sph); 102 /*! \brief Constructor from list of spheres 103 * \param[in] sph the list of spheres 104 * 105 * Only used when we have to deal with a single sphere, i.e. in the unit tests 106 */ 107 ICavity(const std::vector<Sphere> & sph); 108 /*! \brief Constructor from Molecule 109 * \param[in] molec the molecular aggregate 110 */ 111 ICavity(const Molecule & molec); ~ICavity()112 virtual ~ICavity() {} elementCenter()113 Eigen::Matrix3Xd & elementCenter() { return elementCenter_; } elementCenter() const114 const Eigen::Matrix3Xd & elementCenter() const { return elementCenter_; } elementCenter(int i)115 Eigen::Vector3d elementCenter(int i) { return elementCenter_.col(i); } elementCenter(int i) const116 Eigen::Vector3d elementCenter(int i) const { return elementCenter_.col(i); } elementNormal()117 Eigen::Matrix3Xd & elementNormal() { return elementNormal_; } elementNormal() const118 const Eigen::Matrix3Xd & elementNormal() const { return elementNormal_; } elementNormal(int i)119 Eigen::Vector3d elementNormal(int i) { return elementNormal_.col(i); } elementNormal(int i) const120 Eigen::Vector3d elementNormal(int i) const { return elementNormal_.col(i); } elementArea()121 Eigen::VectorXd & elementArea() { return elementArea_; } elementArea() const122 const Eigen::VectorXd & elementArea() const { return elementArea_; } elementArea(int i)123 double elementArea(int i) { return elementArea_(i); } elementArea(int i) const124 double elementArea(int i) const { return elementArea_(i); } size()125 PCMSolverIndex size() { return nElements_; } size() const126 PCMSolverIndex size() const { return nElements_; } irreducible_size()127 PCMSolverIndex irreducible_size() { return nIrrElements_; } irreducible_size() const128 PCMSolverIndex irreducible_size() const { return nIrrElements_; } pointGroup() const129 Symmetry pointGroup() const { return molecule_.pointGroup(); } spheres()130 std::vector<Sphere> & spheres() { return spheres_; } spheres() const131 const std::vector<Sphere> & spheres() const { return spheres_; } nSpheres()132 int nSpheres() { return nSpheres_; } nSpheres() const133 int nSpheres() const { return nSpheres_; } sphereRadius()134 Eigen::VectorXd & sphereRadius() { return sphereRadius_; } sphereRadius() const135 const Eigen::VectorXd & sphereRadius() const { return sphereRadius_; } sphereCenter()136 Eigen::Matrix3Xd & sphereCenter() { return sphereCenter_; } sphereCenter() const137 const Eigen::Matrix3Xd & sphereCenter() const { return sphereCenter_; } elementRadius()138 Eigen::VectorXd & elementRadius() { return elementRadius_; } elementRadius() const139 const Eigen::VectorXd & elementRadius() const { return elementRadius_; } elementRadius(int i)140 double elementRadius(int i) { return elementRadius_(i); } elementRadius(int i) const141 double elementRadius(int i) const { return elementRadius_(i); } elementSphereCenter()142 Eigen::Matrix3Xd & elementSphereCenter() { return elementSphereCenter_; } elementSphereCenter() const143 const Eigen::Matrix3Xd & elementSphereCenter() const { 144 return elementSphereCenter_; 145 } isBuilt()146 bool isBuilt() { return built; } elements() const147 const std::vector<Element> & elements() const { return elements_; } elements(int i) const148 const Element & elements(int i) const { return elements_[i]; } 149 /*! \brief Save cavity specification to file. 150 * 151 * The cavity specification contains: 152 * 0. the number of finite elements, nElements; 153 * 1. the weight of the finite elements, elementArea; 154 * 2. the radius of the finite elements, elementRadius; 155 * 3. the centers of the finite elements, elementCenter; 156 * 4. the normal vectors relative to the centers, elementNormal. 157 * Each of these objects is saved in a separate .npy binary file 158 * and compressed into one .npz file. 159 * Notice that this is just the minimal set of data needed to 160 * restart an energy calculation. 161 */ 162 virtual void saveCavity(const std::string & fname = "cavity.npz"); 163 /*! \brief Load cavity specification from file. 164 */ 165 virtual void loadCavity(const std::string & fname = "cavity.npz"); operator <<(std::ostream & os,ICavity & cavity)166 friend std::ostream & operator<<(std::ostream & os, ICavity & cavity) { 167 return cavity.printCavity(os); 168 } 169 }; 170 } // namespace pcm 171