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 "CPCMSolver.hpp"
25
26 #include <fstream>
27 #include <iostream>
28 #include <string>
29
30 #include "Config.hpp"
31
32 #include <Eigen/Core>
33 #include <Eigen/LU>
34
35 #include "SolverData.hpp"
36 #include "bi_operators/IBoundaryIntegralOperator.hpp"
37 #include "cavity/Element.hpp"
38 #include "cavity/ICavity.hpp"
39 #include "green/IGreensFunction.hpp"
40 #include "utils/MathUtils.hpp"
41
42 namespace pcm {
43 namespace solver {
buildSystemMatrix_impl(const ICavity & cavity,const IGreensFunction & gf_i,const IGreensFunction & gf_o,const IBoundaryIntegralOperator & op)44 void CPCMSolver::buildSystemMatrix_impl(const ICavity & cavity,
45 const IGreensFunction & gf_i,
46 const IGreensFunction & gf_o,
47 const IBoundaryIntegralOperator & op) {
48 if (!isotropic_)
49 PCMSOLVER_ERROR("C-PCM is defined only for isotropic environments!");
50 TIMER_ON("Computing S");
51 double epsilon = gf_o.permittivity();
52 S_ = op.computeS(cavity, gf_i);
53 S_ /= (epsilon - 1.0) / (epsilon + correction_);
54 // Get in Hermitian form
55 if (hermitivitize_)
56 utils::hermitivitize(S_);
57 TIMER_OFF("Computing S");
58
59 // Symmetry-pack
60 // The number of irreps in the group
61 int nrBlocks = cavity.pointGroup().nrIrrep();
62 // The size of the irreducible portion of the cavity
63 int dimBlock = cavity.irreducible_size();
64 utils::symmetryPacking(blockS_, S_, dimBlock, nrBlocks);
65
66 built_ = true;
67 }
68
computeCharge_impl(const Eigen::VectorXd & potential,int irrep) const69 Eigen::VectorXd CPCMSolver::computeCharge_impl(const Eigen::VectorXd & potential,
70 int irrep) const {
71 // The potential and charge vector are of dimension equal to the
72 // full dimension of the cavity. We have to select just the part
73 // relative to the irrep needed.
74 int fullDim = S_.rows();
75 Eigen::VectorXd charge = Eigen::VectorXd::Zero(fullDim);
76 int nrBlocks = blockS_.size();
77 int irrDim = fullDim / nrBlocks;
78 charge.segment(irrep * irrDim, irrDim) =
79 -blockS_[irrep].lu().solve(potential.segment(irrep * irrDim, irrDim));
80
81 return charge;
82 }
83
printSolver(std::ostream & os)84 std::ostream & CPCMSolver::printSolver(std::ostream & os) {
85 os << "Solver Type: C-PCM" << std::endl;
86 if (hermitivitize_) {
87 os << "PCM matrix hermitivitized" << std::endl;
88 } else {
89 os << "PCM matrix NOT hermitivitized (matches old DALTON)" << std::endl;
90 }
91 os << "Correction = " << correction_;
92
93 return os;
94 }
95
createCPCMSolver(const SolverData & data)96 ISolver * createCPCMSolver(const SolverData & data) {
97 return new CPCMSolver(data.hermitivitize, data.correction);
98 }
99 } // namespace solver
100 } // namespace pcm
101