Add an additional molecule (H2O) for testing

Author: Roberto Di Remigio

tests/TestingMolecules.hpp

@@ -2,22 +2,22 @@
 /*
  *     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 <http://www.gnu.org/licenses/>.
- *     
+ *
  *     For information on the complete list of contributors to the
  *     PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
  */
@@ -73,6 +73,10 @@ inline Molecule C2H4();
  */
 inline Molecule C6H6();
 
+/*! Returns the water molecule with C1 symmetry
+ */
+inline Molecule H2O();
+
 template <int group>
 Molecule dummy(double radius, const Eigen::Vector3d & center)
 {
@@ -508,10 +512,47 @@ Molecule C6H6()
     spheres.push_back(sph11);
     spheres.push_back(sph12);
 
-    // D2h as generated by Oxy, Oxz, Oyz
     Symmetry pGroup = buildGroup(0, 0, 0, 0);
 
     return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
 };
 
+Molecule H2O()
+{
+    int nAtoms = 3;
+
+    // Geometry in atomic units
+    Eigen::Vector3d H1( 0.0000000000, 0.0000000000, -0.2249058930);
+    Eigen::Vector3d H2(-1.4523499293, 0.0000000000,  0.8996235720);
+    Eigen::Vector3d O(  1.4523499293, 0.0000000000,  0.8996235720);
+
+    Eigen::MatrixXd geom(3, nAtoms);
+    geom.col(0) = H1.transpose();
+    geom.col(1) = H2.transpose();
+    geom.col(2) = O.transpose();
+    Eigen::Vector3d charges, masses;
+    charges << 1.0, 1.0, 8.0;
+    masses  << 1.0078250, 1.0078250, 15.9949150;
+
+    double radiusO = (1.52 * 1.20) / bohrToAngstrom();
+    double radiusH = (1.20 * 1.20) / bohrToAngstrom();
+    std::vector<Atom> atoms;
+    atoms.push_back( Atom("Hydrogen", "H", charges(0), masses(0), radiusH, H1, 1.0) );
+    atoms.push_back( Atom("Hydrogen", "H", charges(1), masses(1), radiusH, H2, 1.0) );
+    atoms.push_back( Atom("Oxygen",   "O", charges(2), masses(2), radiusO, O, 1.0) );
+
+    std::vector<Sphere> spheres;
+    Sphere sph1(H1, radiusH);
+    Sphere sph2(H2, radiusH);
+    Sphere sph3(O, radiusO);
+    spheres.push_back(sph1);
+    spheres.push_back(sph2);
+    spheres.push_back(sph3);
+
+    Symmetry pGroup = buildGroup(0, 0, 0, 0);
+
+    return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
+};
+
+
 #endif // TESTINGMOLECULES_HPP