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 <string> 28 #include <vector> 29 30 #include "Config.hpp" 31 #include "PCMSolverExport.h" 32 33 #include <Eigen/Core> 34 35 /*! \file Molecule.hpp */ 36 37 namespace pcm { 38 namespace cavity { 39 class Element; 40 } // namespace cavity 41 } // namespace pcm 42 43 #include "Atom.hpp" 44 #include "Sphere.hpp" 45 #include "Symmetry.hpp" 46 47 namespace pcm { 48 using utils::Atom; 49 using utils::Sphere; 50 51 enum rotorType { rtAsymmetric, rtSymmetric, rtSpherical, rtLinear, rtAtom }; 52 const std::string rotorTypeList[] = {"Asymmetric", 53 "Symmetric", 54 "Spherical", 55 "Linear", 56 "Atom"}; 57 58 /*! \class Molecule 59 * \brief Class representing a molecule or general aggregate of atoms. 60 * \author Roberto Di Remigio 61 * \date 2014 62 * 63 * This class is based on the similar class available in the Mints library 64 * of Psi4 65 */ 66 67 class Molecule { 68 private: 69 /// The number of atoms in the molecule 70 size_t nAtoms_; 71 /// A vector of dimension (# atoms) containing the charges 72 Eigen::VectorXd charges_; 73 /// A vector of dimension (# atoms) containing the masses 74 Eigen::VectorXd masses_; 75 /// Molecular geometry, in cartesian coordinates. The dimensions are (# atoms * 3) 76 /// Units are Bohr. 77 Eigen::Matrix3Xd geometry_; 78 /// A container for all the atoms composing the molecule 79 std::vector<Atom> atoms_; 80 /// A container for the spheres composing the molecule 81 std::vector<Sphere> spheres_; 82 /// The molecular rotor type 83 rotorType rotor_; 84 /// The molecular point group 85 Symmetry pointGroup_; 86 87 public: 88 /*! \brief Default constructor 89 * Initialize a dummy molecule, e.g. as placeholder, see ICavity.cpp loadCavity 90 * method 91 */ Molecule()92 Molecule() : rotor_(rtAsymmetric) {} 93 /*! \brief Constructor from full molecular data 94 * \param[in] nat number of atoms 95 * \param[in] chg vector of atomic charges 96 * \param[in] masses vector of atomic masses 97 * \param[in] geo molecular geometry (format nat*3) 98 * \param[in] at vector of Atom objects 99 * \param[in] sph vector of Sphere objects 100 * 101 * This initializes the molecule in C1 symmetry 102 */ 103 Molecule(int nat, 104 const Eigen::VectorXd & chg, 105 const Eigen::VectorXd & masses, 106 const Eigen::Matrix3Xd & geo, 107 const std::vector<Atom> & at, 108 const std::vector<Sphere> & sph); 109 /*! \brief Constructor from full molecular data, plus number of generators and 110 *generators 111 * \param[in] nat number of atoms 112 * \param[in] chg vector of atomic charges 113 * \param[in] masses vector of atomic masses 114 * \param[in] geo molecular geometry (format nat*3) 115 * \param[in] at vector of Atom objects 116 * \param[in] sph vector of Sphere objects 117 * \param[in] nr_gen number of molecular point group generators 118 * \param[in] gen molecular point group generators 119 * 120 * This initializes the molecule in the symmetry prescribed by nr_gen and gen. 121 * See documentation of the Symmetry object for the conventions. 122 */ 123 Molecule(int nat, 124 const Eigen::VectorXd & chg, 125 const Eigen::VectorXd & masses, 126 const Eigen::Matrix3Xd & geo, 127 const std::vector<Atom> & at, 128 const std::vector<Sphere> & sph, 129 int nr_gen, 130 std::array<int, 3> gens); 131 /*! \brief Constructor from full molecular data and point group 132 * \param[in] nat number of atoms 133 * \param[in] chg vector of atomic charges 134 * \param[in] masses vector of atomic masses 135 * \param[in] geo molecular geometry (format nat*3) 136 * \param[in] at vector of Atom objects 137 * \param[in] sph vector of Sphere objects 138 * \param[in] pg the molecular point group (a Symmetry object) 139 * 140 * This initializes the molecule in the symmetry prescribed by pg. 141 */ 142 Molecule(int nat, 143 const Eigen::VectorXd & chg, 144 const Eigen::VectorXd & masses, 145 const Eigen::Matrix3Xd & geo, 146 const std::vector<Atom> & at, 147 const std::vector<Sphere> & sph, 148 const Symmetry & pg); 149 /*! \brief Constructor from list of spheres 150 * \param[in] sph list of spheres 151 * 152 * \warning This constructor is to be used **exclusively** when initializing the 153 *Molecule 154 * in EXPLICIT mode, i.e. when the user specifies explicitly spheres centers and 155 *radii. 156 * 157 * Molecule is treated as an aggregate of spheres. We do not have information on 158 *the atomic 159 * species involved in the aggregate. 160 * Charges are set to 1.0; masses are set based on the radii; geometry is set from 161 *the list of spheres. 162 * All the atoms are dummy atoms. The point group is C1. 163 */ 164 Molecule(const std::vector<Sphere> & sph); 165 /// Copy constructor. 166 Molecule(const Molecule & other); ~Molecule()167 ~Molecule() {} 168 nAtoms() const169 size_t nAtoms() const { return nAtoms_; } charges() const170 Eigen::VectorXd charges() const { return charges_; } charges(int i) const171 double charges(int i) const { return charges_(i); } masses() const172 Eigen::VectorXd masses() const { return masses_; } masses(int i) const173 double masses(int i) const { return masses_(i); } geometry() const174 Eigen::Matrix3Xd geometry() const { return geometry_; } geometry(int i,int j) const175 double geometry(int i, int j) const { return geometry_(i, j); } atoms() const176 std::vector<Atom> atoms() const { return atoms_; } atoms(int i) const177 Atom atoms(int i) const { return atoms_[i]; } spheres() const178 std::vector<Sphere> spheres() const { return spheres_; } spheres(int i) const179 Sphere spheres(int i) const { return spheres_[i]; } 180 181 rotorType rotor(); 182 rotorType findRotorType(); 183 pointGroup() const184 Symmetry pointGroup() const { return pointGroup_; } pointGroup(const Symmetry & pg)185 void pointGroup(const Symmetry & pg) { pointGroup_ = pg; } 186 187 Eigen::Vector3d centerOfMass(); 188 Eigen::Matrix3d inertiaTensor(); 189 190 /*! 191 * Given a vector, carries out translation of the molecule. 192 * \param translationVector The translation vector. 193 */ 194 void translate(const Eigen::Vector3d & translationVector); 195 /// Performs translation to the Center of Mass Frame. 196 void moveToCOM(); 197 198 /*! 199 * Given a matrix, carries out rotation of the molecule. 200 * \param rotationMatrix The matrix representing the rotation. 201 */ 202 void rotate(const Eigen::Matrix3d & rotationMatrix); 203 /// Performs rotation to the Principal Axes Frame. 204 void moveToPAF(); 205 206 /// @{ 207 /// Operators 208 /// Assignment operator. 209 Molecule & operator=(const Molecule & other); 210 friend std::ostream & operator<<(std::ostream & os, const Molecule & molecule); 211 /// @} 212 EIGEN_MAKE_ALIGNED_OPERATOR_NEW 213 }; 214 215 /*! \brief Compute MEP from molecule object and list of finite elements 216 * \param[in] mol molecule object 217 * \param[in] el list of finite elements 218 * \return MEP at finite elements center 219 */ 220 Eigen::VectorXd computeMEP(const Molecule & mol, 221 const std::vector<cavity::Element> & el); 222 223 /*! \brief Compute MEP from molecule object and a grid of points 224 * \param[in] mol molecule object 225 * \param[in] grid grid points coordinates 226 * \return MEP at grid points 227 */ 228 PCMSolver_EXPORT Eigen::VectorXd computeMEP(const Molecule & mol, 229 const Eigen::Matrix3Xd & grid); 230 231 /*! \brief Compute MEP for a single point charge 232 * \param[in] el list of finite elements 233 * \param[in] charge value of the charge 234 * \param[in] origin location of the point charge 235 * \return MEP at finite elements center 236 */ 237 Eigen::VectorXd computeMEP(const std::vector<cavity::Element> & el, 238 double charge = 1.0, 239 const Eigen::Vector3d & origin = Eigen::Vector3d::Zero()); 240 241 /*! Gauss' theorem estimate of the total ASC for a set of point charges 242 * \param[in] charges Vector of point charges 243 * \param[in] permittivity Solvent permittivity 244 * \param[in] correction The CPCM correction factor 245 */ 246 PCMSolver_EXPORT double GaussEstimate(const Eigen::VectorXd & charges, 247 double permittivity, 248 double correction = 0.0); 249 } // namespace pcm 250