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 #include <catch.hpp> 25 26 #include <cmath> 27 #include <vector> 28 29 #include <Eigen/Core> 30 31 #include "TestingMolecules.hpp" 32 #include "cavity/GePolCavity.hpp" 33 #include "utils/Molecule.hpp" 34 #include "utils/Symmetry.hpp" 35 36 using namespace pcm; 37 using cavity::GePolCavity; 38 39 SCENARIO("GePol cavity for the C2H4 molecule in D2h symmetry", 40 "[gepol][gepol_C2H4_D2h]") { 41 GIVEN("The C2H4 molecule in D2h symmetry") { 42 Molecule molec = C2H4(); 43 44 WHEN("the addition of spheres is enabled") { 45 double area = 0.2 / bohr2ToAngstrom2(); 46 double probeRadius = 1.385 / bohrToAngstrom(); 47 double minRadius = 0.2 / bohrToAngstrom(); 48 GePolCavity cavity(molec, area, probeRadius, minRadius, "d2h"); 49 cavity.saveCavity("c2h4_d2h.npz"); 50 51 /*! \class GePolCavity 52 * \test \b GePolCavityD2hAddTest_size tests GePol cavity size for C2H4 in D2h 53 * symmetry with added spheres 54 */ 55 THEN("the size of the cavity is") { 56 int size = 576; 57 int actualSize = cavity.size(); 58 REQUIRE(size == actualSize); 59 } 60 /*! \class GePolCavity 61 * \test \b GePolCavityD2hAddTest_irreducible_size tests GePol cavity 62 * irreducible size for C2H4 in D2h symmetry with added spheres 63 */ 64 AND_THEN("the irreducible size of the cavity is") { 65 int size = 72; 66 int actualSize = cavity.irreducible_size(); 67 REQUIRE(size == actualSize); 68 } 69 /*! \class GePolCavity 70 * \test \b GePolCavityD2hAddTest_area tests GePol cavity surface area for 71 * C2H4 in D2h symmetry with added spheres 72 */ 73 AND_THEN("the surface area of the cavity is") { 74 double area = 281.81993683500656; 75 double actualArea = cavity.elementArea().sum(); 76 REQUIRE(area == Approx(actualArea)); 77 } 78 /*! \class GePolCavity 79 * \test \b GePolCavityD2hAddTest_volume tests GePol cavity volume for C2H4 in 80 * D2h symmetry with added spheres 81 */ 82 AND_THEN("the volume of the cavity is") { 83 double volume = 406.54737252764619; 84 Eigen::Matrix3Xd elementCenter = cavity.elementCenter(); 85 Eigen::Matrix3Xd elementNormal = cavity.elementNormal(); 86 double actualVolume = 0; 87 for (int i = 0; i < cavity.size(); ++i) { 88 actualVolume += 89 cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i)); 90 } 91 actualVolume /= 3; 92 REQUIRE(volume == Approx(actualVolume)); 93 } 94 } 95 96 WHEN("the addition of spheres is disabled") { 97 double area = 0.2 / bohr2ToAngstrom2(); 98 double probeRadius = 1.385 / bohrToAngstrom(); 99 double minRadius = 100.0 / bohrToAngstrom(); 100 GePolCavity cavity(molec, area, probeRadius, minRadius, "d2h_noadd"); 101 cavity.saveCavity("c2h4_d2h_noadd.npz"); 102 103 /*! \class GePolCavity 104 * \test \b GePolCavityD2hTest_size tests GePol cavity size for C2H4 in D2h 105 * symmetry without added spheres 106 */ 107 THEN("the size of the cavity is") { 108 int size = 576; 109 int actualSize = cavity.size(); 110 REQUIRE(size == actualSize); 111 } 112 /*! \class GePolCavity 113 * \test \b GePolCavityD2hTest_irreducible_size tests GePol cavity irreducible 114 * size for C2H4 in D2h symmetry without added spheres 115 */ 116 AND_THEN("the irreducible size of the cavity is") { 117 int size = 72; 118 int actualSize = cavity.irreducible_size(); 119 REQUIRE(size == actualSize); 120 } 121 /*! \class GePolCavity 122 * \test \b GePolCavityD2hTest_area tests GePol cavity surface area for C2H4 123 * in D2h symmetry without added spheres 124 */ 125 AND_THEN("the surface area of the cavity is") { 126 double area = 281.81993683500656; 127 double actualArea = cavity.elementArea().sum(); 128 REQUIRE(area == Approx(actualArea)); 129 } 130 /*! \class GePolCavity 131 * \test \b GePolCavityD2hTest_volume tests GePol cavity volume for C2H4 in 132 * D2h symmetry without added spheres 133 */ 134 AND_THEN("the volume of the cavity is") { 135 double volume = 406.54737252764619; 136 Eigen::Matrix3Xd elementCenter = cavity.elementCenter(); 137 Eigen::Matrix3Xd elementNormal = cavity.elementNormal(); 138 double actualVolume = 0; 139 for (int i = 0; i < cavity.size(); ++i) { 140 actualVolume += 141 cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i)); 142 } 143 actualVolume /= 3; 144 REQUIRE(volume == Approx(actualVolume)); 145 } 146 } 147 } 148 } 149