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