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/CPCMSolver.hpp" 36 37 using namespace pcm; 38 using bi_operators::Collocation; 39 using cavity::GePolCavity; 40 using green::UniformDielectric; 41 using green::Vacuum; 42 using solver::CPCMSolver; 43 44 /*! \class CPCMSolver 45 * \test \b NH3GePolRestart tests CPCMSolver using ammonia with a GePol cavity read 46 * from .npz file 47 */ 48 TEST_CASE("Test solver for the C-PCM for NH3 and a restarted GePol cavity", 49 "[solver][cpcm][cpcm_gepol-NH3_from-file]") { 50 // Set up cavity 51 Eigen::Vector3d N(-0.000000000, -0.104038047, 0.000000000); 52 Eigen::Vector3d H1(-0.901584415, 0.481847022, -1.561590016); 53 Eigen::Vector3d H2(-0.901584415, 0.481847022, 1.561590016); 54 Eigen::Vector3d H3(1.803168833, 0.481847022, 0.000000000); 55 GePolCavity cavity; 56 cavity.loadCavity("nh3.npz"); 57 58 double permittivity = 78.39; 59 Vacuum<> gf_i; 60 UniformDielectric<> gf_o(permittivity); 61 bool symm = true; 62 double correction = 0.0; 63 64 Collocation S; 65 66 CPCMSolver solver(symm, correction); 67 solver.buildSystemMatrix(cavity, gf_i, gf_o, S); 68 69 double Ncharge = 7.0; 70 double Hcharge = 1.0; 71 int size = cavity.size(); 72 Eigen::VectorXd fake_mep = Eigen::VectorXd::Zero(size); 73 for (int i = 0; i < size; ++i) { 74 Eigen::Vector3d center = cavity.elementCenter(i); 75 double Ndistance = (center - N).norm(); 76 double H1distance = (center - H1).norm(); 77 double H2distance = (center - H2).norm(); 78 double H3distance = (center - H3).norm(); 79 fake_mep(i) = Ncharge / Ndistance + Hcharge / H1distance + Hcharge / H2distance + 80 Hcharge / H3distance; 81 } 82 // The total ASC for a conductor is -Q 83 // for CPCM it will be -Q*(epsilon-1)/(epsilon + correction) 84 Eigen::VectorXd fake_asc = Eigen::VectorXd::Zero(size); 85 fake_asc = solver.computeCharge(fake_mep); 86 double totalASC = -(Ncharge + 3.0 * Hcharge) * (permittivity - 1) / permittivity; 87 double totalFakeASC = fake_asc.sum(); 88 CAPTURE(totalASC - totalFakeASC); 89 REQUIRE(totalASC == Approx(totalFakeASC).epsilon(1.0e-03)); 90 } 91