Skip to content

Commit 46c86f2

Browse files
author
Roberto Di Remigio
committed
Add stub for Purisima-style integration of diagonal of D
1 parent 888ca85 commit 46c86f2

File tree

4 files changed

+237
-1
lines changed

4 files changed

+237
-1
lines changed

‎doc/pcmsolver.bib‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -263,3 +263,14 @@ @Article{Cances1998
263263
Publisher = {Springer Netherlands},
264264
Timestamp = {17.02.2011}
265265
}
266+
267+
@article{Purisima1995,
268+
author = {{Purisima, E. O. and Nilar, S. H.}},
269+
journal = {{J. Comput. Chem.}},
270+
number = {6},
271+
pages = {681--689},
272+
title = {{A Simple Yet Accurate Boundary Element Method for Continuum Dielectric Calculations}},
273+
url = {http://onlinelibrary.wiley.com/doi/10.1002/jcc.540160604/abstract},
274+
volume = {16},
275+
year = {1995}
276+
}

‎src/bi_operators/IntegratorForward.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,5 +27,7 @@
2727
#define INTEGRATORFORWARD_HPP
2828

2929
struct CollocationIntegrator;
30+
struct NumericalIntegrator;
31+
struct PurisimaIntegrator;
3032

3133
#endif // INTEGRATORFORWARD_HPP

‎src/bi_operators/IntegratorTypes.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,9 @@
3131
#include <boost/mpl/vector.hpp>
3232

3333
#include "CollocationIntegrator.hpp"
34+
#include "PurisimaIntegator.hpp"
35+
#include "NumericalIntegrator.hpp"
3436

35-
typedef boost::mpl::vector<CollocationIntegrator> integrator_types;
37+
typedef boost::mpl::vector<CollocationIntegrator, PurisimaIntegrator, NumericalIntegrator> integrator_types;
3638

3739
#endif // INTEGRATORTYPES_HPP

‎src/bi_operators/PurisimaIntegrator.hpp‎

Lines changed: 221 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,221 @@
1+
/* pcmsolver_copyright_start */
2+
/*
3+
* PCMSolver, an API for the Polarizable Continuum Model
4+
* Copyright (C) 2013 Roberto Di Remigio, Luca Frediani and contributors
5+
*
6+
* This file is part of PCMSolver.
7+
*
8+
* PCMSolver is free software: you can redistribute it and/or modify
9+
* it under the terms of the GNU Lesser General Public License as published by
10+
* the Free Software Foundation, either version 3 of the License, or
11+
* (at your option) any later version.
12+
*
13+
* PCMSolver is distributed in the hope that it will be useful,
14+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
15+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16+
* GNU Lesser General Public License for more details.
17+
*
18+
* You should have received a copy of the GNU Lesser General Public License
19+
* along with PCMSolver. If not, see <http://www.gnu.org/licenses/>.
20+
*
21+
* For information on the complete list of contributors to the
22+
* PCMSolver API, see: <http://pcmsolver.github.io/pcmsolver-doc>
23+
*/
24+
/* pcmsolver_copyright_end */
25+
26+
#ifndef PURISIMAINTEGRATOR_HPP
27+
#define PURISIMAINTEGRATOR_HPP
28+
29+
#include <cmath>
30+
#include <functional>
31+
#include <iosfwd>
32+
#include <stdexcept>
33+
#include <tuple>
34+
#include <vector>
35+
36+
#include "Config.hpp"
37+
38+
#include <Eigen/Dense>
39+
40+
#include "IntegratorHelperFunctions.hpp"
41+
#include "Element.hpp"
42+
#include "AnisotropicLiquid.hpp"
43+
#include "IonicLiquid.hpp"
44+
#include "SphericalDiffuse.hpp"
45+
#include "UniformDielectric.hpp"
46+
#include "Vacuum.hpp"
47+
48+
/*! \file PurisimaIntegrator.hpp
49+
* \struct PurisimaIntegrator
50+
* \brief Implementation of the single and double layer operators matrix representation using
51+
* one-point collocation and Purisima's strategy for the diagonal of D
52+
* \author Roberto Di Remigio
53+
* \date 2015
54+
*
55+
* Calculates the diagonal elements of S as:
56+
* \f[
57+
* S_{ii} = factor_ * \sqrt{\frac{4\pi}{a_i}}
58+
* \f]
59+
* while the diagonal elements of D are:
60+
* \f[
61+
* D_{ii} = -\left(2\pi + \sum_{j\neq i}D_{ij}a_j \right)\frac{1}{a_i}
62+
* \f]
63+
* The original reference is \cite Purisima1995
64+
*/
65+
66+
using std::placeholders::_1;
67+
using std::placeholders::_2;
68+
using std::placeholders::_3;
69+
70+
struct PurisimaIntegrator
71+
{
72+
PurisimaIntegrator() : factor_(1.07) {}
73+
~PurisimaIntegrator() {}
74+
75+
/**@{ Single and double layer potentials for a Vacuum Green's function by collocation */
76+
/*! \tparam DerivativeTraits how the derivative of the Greens's function are calculated
77+
* \param[in] gf Green's function
78+
* \param[in] e list of finite elements
79+
*/
80+
template <typename DerivativeTraits>
81+
Eigen::MatrixXd singleLayer(const Vacuum<DerivativeTraits, PurisimaIntegrator> & gf, const std::vector<Element> & e) const {
82+
auto f = this->factor_;
83+
auto diagonal = [f] (const Element & el) -> double { return (f * std::sqrt(4 * M_PI / el.area())); };
84+
auto kernelS = std::bind(&Vacuum<DerivativeTraits, PurisimaIntegrator>::kernelS, gf, _1, _2);
85+
return integrator::singleLayer(e, diagonal, kernelS);
86+
}
87+
/*! \tparam DerivativeTraits how the derivative of the Greens's function are calculated
88+
* \param[in] gf Green's function
89+
* \param[in] e list of finite elements
90+
*/
91+
template <typename DerivativeTraits>
92+
Eigen::MatrixXd doubleLayer(const Vacuum<DerivativeTraits, PurisimaIntegrator> & gf, const std::vector<Element> & e) const {
93+
auto f = this->factor_;
94+
auto diagonal = [f] (const Element & el) -> double { return (-f * std::sqrt(M_PI/ el.area()) * (1.0 / el.sphere().radius())); };
95+
auto kernelD = std::bind(&Vacuum<DerivativeTraits, PurisimaIntegrator>::kernelD, gf, _1, _2, _3);
96+
return integrator::doubleLayer(e, diagonal, kernelD);
97+
}
98+
/**@}*/
99+
100+
/**@{ Single and double layer potentials for a UniformDielectric Green's function by collocation */
101+
/*! \tparam DerivativeTraits how the derivative of the Greens's function are calculated
102+
* \param[in] gf Green's function
103+
* \param[in] e list of finite elements
104+
*/
105+
template <typename DerivativeTraits>
106+
Eigen::MatrixXd singleLayer(const UniformDielectric<DerivativeTraits, PurisimaIntegrator> & gf, const std::vector<Element> & e) const {
107+
auto f = this->factor_;
108+
auto epsInv = 1.0 / gf.epsilon();
109+
auto diagonal = [f, epsInv] (const Element & el) -> double { return (f * std::sqrt(4 * M_PI / el.area()) * epsInv); };
110+
auto kernelS = std::bind(&UniformDielectric<DerivativeTraits, PurisimaIntegrator>::kernelS, gf, _1, _2);
111+
return integrator::singleLayer(e, diagonal, kernelS);
112+
}
113+
/*! \tparam DerivativeTraits how the derivative of the Greens's function are calculated
114+
* \param[in] gf Green's function
115+
* \param[in] e list of finite elements
116+
*/
117+
template <typename DerivativeTraits>
118+
Eigen::MatrixXd doubleLayer(const UniformDielectric<DerivativeTraits, PurisimaIntegrator> & gf, const std::vector<Element> & e) const {
119+
auto f = this->factor_;
120+
auto diagonal = [f] (const Element & el) -> double { return (-f * std::sqrt(M_PI/ el.area()) * (1.0 / el.sphere().radius())); };
121+
auto kernelD = std::bind(&UniformDielectric<DerivativeTraits, PurisimaIntegrator>::kernelD, gf, _1, _2, _3);
122+
return integrator::doubleLayer(e, diagonal, kernelD);
123+
}
124+
/**@}*/
125+
126+
/**@{ Single and double layer potentials for a IonicLiquid Green's function by collocation */
127+
template <typename DerivativeTraits>
128+
Eigen::MatrixXd singleLayer(const IonicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
129+
throw std::runtime_error("Not implemented yet!");
130+
}
131+
template <typename DerivativeTraits>
132+
Eigen::MatrixXd doubleLayer(const IonicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
133+
throw std::runtime_error("Not implemented yet!");
134+
}
135+
/**@}*/
136+
137+
/**@{ Single and double layer potentials for an AnisotropicLiquid Green's function by collocation */
138+
template <typename DerivativeTraits>
139+
Eigen::MatrixXd singleLayer(const AnisotropicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
140+
throw std::runtime_error("Not implemented yet!");
141+
}
142+
template <typename DerivativeTraits>
143+
Eigen::MatrixXd doubleLayer(const AnisotropicLiquid<DerivativeTraits, PurisimaIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
144+
throw std::runtime_error("Not implemented yet!");
145+
}
146+
/**@}*/
147+
148+
/**@{ Single and double layer potentials for a SphericalDiffuse Green's function by collocation */
149+
/*! \tparam ProfilePolicy the permittivity profile for the diffuse interface
150+
* \param[in] gf Green's function
151+
* \param[in] e list of finite elements
152+
*/
153+
template <typename ProfilePolicy>
154+
Eigen::MatrixXd singleLayer(const SphericalDiffuse<PurisimaIntegrator, ProfilePolicy> & gf, const std::vector<Element> & e) const {
155+
// The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
156+
size_t mat_size = e.size();
157+
Eigen::MatrixXd S = Eigen::MatrixXd::Zero(mat_size, mat_size);
158+
for (size_t i = 0; i < mat_size; ++i) {
159+
// Fill diagonal
160+
// Diagonal of S inside the cavity
161+
double Sii_I = factor_ * std::sqrt(4 * M_PI / e[i].area());
162+
// "Diagonal" of Coulomb singularity separation coefficient
163+
double coulomb_coeff = gf.coefficientCoulomb(e[i].center(), e[i].center());
164+
// "Diagonal" of the image Green's function
165+
double image = gf.imagePotential(e[i].center(), e[i].center());
166+
S(i, i) = Sii_I / coulomb_coeff + image;
167+
Eigen::Vector3d source = e[i].center();
168+
for (size_t j = 0; j < mat_size; ++j) {
169+
// Fill off-diagonal
170+
Eigen::Vector3d probe = e[j].center();
171+
if (i != j) S(i, j) = gf.kernelS(source, probe);
172+
}
173+
}
174+
return S;
175+
}
176+
/*! \tparam ProfilePolicy the permittivity profile for the diffuse interface
177+
* \param[in] gf Green's function
178+
* \param[in] e list of finite elements
179+
*/
180+
template <typename ProfilePolicy>
181+
Eigen::MatrixXd doubleLayer(const SphericalDiffuse<PurisimaIntegrator, ProfilePolicy> & gf, const std::vector<Element> & e) const {
182+
// The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
183+
size_t mat_size = e.size();
184+
Eigen::MatrixXd D = Eigen::MatrixXd::Zero(mat_size, mat_size);
185+
for (size_t i = 0; i < mat_size; ++i) {
186+
// Fill diagonal
187+
double area = e[i].area();
188+
double radius = e[i].sphere().radius();
189+
// Diagonal of S inside the cavity
190+
double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
191+
// Diagonal of D inside the cavity
192+
double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
193+
// "Diagonal" of Coulomb singularity separation coefficient
194+
double coulomb_coeff = gf.coefficientCoulomb(e[i].center(), e[i].center());
195+
// "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
196+
double coeff_grad = gf.coefficientCoulombDerivative(e[i].normal(), e[i].center(), e[i].center()) / std::pow(coulomb_coeff, 2);
197+
// "Diagonal" of the directional derivative of the image Green's function
198+
double image_grad = gf.imagePotentialDerivative(e[i].normal(), e[i].center(), e[i].center());
199+
200+
double eps_r2 = 0.0;
201+
std::tie(eps_r2, std::ignore) = gf.epsilon(e[i].center());
202+
203+
D(i, i) = eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad);
204+
Eigen::Vector3d source = e[i].center();
205+
for (size_t j = 0; j < mat_size; ++j) {
206+
// Fill off-diagonal
207+
Eigen::Vector3d probe = e[j].center();
208+
Eigen::Vector3d probeNormal = e[j].normal();
209+
probeNormal.normalize();
210+
if (i != j) D(i, j) = gf.kernelD(probeNormal, source, probe);
211+
}
212+
}
213+
return D;
214+
}
215+
/**@}*/
216+
217+
/// Scaling factor for the collocation formulas
218+
double factor_;
219+
};
220+
221+
#endif // PURISIMAINTEGRATOR_HPP

0 commit comments

Comments
 (0)