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 <iostream>
27 
28 #include <Eigen/Core>
29 
30 #include "bi_operators/Collocation.hpp"
31 #include "cavity/GePolCavity.hpp"
32 #include "green/DerivativeTypes.hpp"
33 #include "green/UniformDielectric.hpp"
34 #include "green/Vacuum.hpp"
35 #include "solver/IEFSolver.hpp"
36 #include "utils/Molecule.hpp"
37 
38 using namespace pcm;
39 using bi_operators::Collocation;
40 using cavity::GePolCavity;
41 using green::UniformDielectric;
42 using green::Vacuum;
43 using solver::IEFSolver;
44 
45 /*! \class IEFSolver
46  *  \test \b pointChargeGePolRestart tests IEFSolver using a point charge with a
47  * GePol cavity read from .npz file
48  */
49 TEST_CASE(
50     "Test solver for the IEFPCM for a point charge and a restarted GePol cavity",
51     "[solver][iefpcm][iefpcm_gepol-point_from-file]") {
52   // Load cavity
53   GePolCavity cavity;
54   cavity.loadCavity("point.npz");
55 
56   double permittivity = 78.39;
57   Vacuum<> gf_i;
58   UniformDielectric<> gf_o(permittivity);
59 
60   Collocation op;
61 
62   bool symm = true;
63   IEFSolver solver(symm);
64   solver.buildSystemMatrix(cavity, gf_i, gf_o, op);
65 
66   double charge = 8.0;
67   int size = cavity.size();
68   Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge);
69   for (int i = 0; i < size; ++i) {
70     Eigen::Vector3d center = cavity.elementCenter(i);
71     double distance = center.norm();
72     fake_mep(i) = charge / distance;
73   }
74   Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size);
75   fake_asc = solver.computeCharge(fake_mep);
76 
77   double totalASC = -charge * (permittivity - 1) / permittivity;
78   double totalFakeASC = fake_asc.sum();
79   CAPTURE(totalASC - totalFakeASC);
80   REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03));
81 }
82