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