MADNESS  0.10.1
molecularmask.h
Go to the documentation of this file.
1 #ifndef MOLECULAR_MASK_H
2 #define MOLECULAR_MASK_H
3 
4 #include <madness/mra/mra.h>
5 #include <madness/constants.h>
6 #include <cmath>
7 #include <vector>
8 
9 // Distance between two points in 3D
10 inline double distance(const madness::coord_3d& a, const madness::coord_3d& b) {
11  double x = a[0] - b[0];
12  double y = a[1] - b[1];
13  double z = a[2] - b[2];
14  return sqrt(x*x + y*y + z*z);
15 }
16 
17 // Basic functionality for the mask
19 protected:
20  const double sigma;
21  const std::vector<double> atomic_radii;
22  const std::vector<madness::coord_3d> atomic_coords;
23  const int natom;
24 
25  // signed distance function for point r relative to sphere radius R at center
26  double sdf(const madness::coord_3d& r, const madness::coord_3d& center, double R) const {
27  return distance(r,center) - R;
28  }
29 
30  // gradient of the signed distance function
32  return (r - center)*(1.0/distance(r,center));
33  }
34 
35  // Mask or characteristic function (argument s is the signed distance)
36  double mask(double s) const {
37  if (s > 6.0) return 0.0;
38  else if (s < -6.0) return 1.0;
39  else return 0.5*erfc(s);
40  }
41 
42  // Complement of the mask or characteristic function (argument s is the signed distance)
43  double cmask(double s) const {
44  return mask(-s);
45  }
46 
47  // Derivative of the mask w.r.t. s
48  double dmask(double s) const {
49  const double fac = 1.0/sqrt(madness::constants::pi);
50  if (fabs(s) > 6.0) return 0.0;
51  return -exp(-s*s)*fac;
52  }
53 
54  // Mask or characteristic function for atom i
55  double atomic_mask(const madness::coord_3d& r, unsigned int i) const {
56  double s = sdf(r, atomic_coords[i], atomic_radii[i]);
57  return mask(s/sigma);
58  }
59 
60  // Complement of the mask or characteristic function for atom i
61  // (we use this directly to avoid numerical cancellation)
62  double atomic_cmask(const madness::coord_3d& r, unsigned int i) const {
63  double s = sdf(r, atomic_coords[i], atomic_radii[i]);
64  return cmask(s/sigma);
65  }
66 
67  // Gradient of the atomic mask
68  madness::coord_3d grad_atomic_mask(const madness::coord_3d& r, unsigned int i) const {
69  double s = sdf(r, atomic_coords[i], atomic_radii[i]);
70  return grad_sdf(r,atomic_coords[i])*(dmask(s/sigma)/sigma);
71  }
72 
73  // Gradient of the molecular mask
75  // precompute the atomic masks
76  std::vector<double> m(natom);
77  double value = 1.0;
78  for (int i=0; i<natom; i++) {
79  m[i] = atomic_cmask(r,i);
80  value *= m[i];
81  }
82 
83  // return 0.0 if not in the surface
84  if (value<1e-12 || (1.0-value)<1e-12) return madness::coord_3d(0.0);
85 
87  for (int i=0; i<natom; i++) {
88  if (m[i] > 1e-12)
89  grad += grad_atomic_mask(r,i)*(value/m[i]);
90  }
91  return grad;
92  }
93 
94 public:
96  const std::vector<double> atomic_radii,
97  const std::vector<madness::coord_3d> atomic_coords)
98  : sigma(sigma)
101  , natom(atomic_coords.size())
102  {
103  MADNESS_ASSERT(atomic_radii.size() == atomic_coords.size());
104  }
105 
106  std::vector< madness::Vector<double,3> > special_points() const {
107  std::vector< madness::Vector<double,3> > v;
108  int npt = 4;
109  for (int t=0; t<=npt; t++) {
110  const double theta = madness::constants::pi*double(t)/double(npt);
111  for (int p=0; p<2*npt; p++) {
112  const double phi = 2.0*madness::constants::pi*double(p)/double(2*npt);
113  const double xx = std::sin(theta)*std::cos(phi);
114  const double yy = std::sin(theta)*std::sin(phi);
115  const double zz = std::cos(theta);
116  for (int i=0; i<natom; i++) {
117  const double x = atomic_coords[i][0] + atomic_radii[i]*xx;
118  const double y = atomic_coords[i][1] + atomic_radii[i]*yy;
119  const double z = atomic_coords[i][2] + atomic_radii[i]*zz;
120  v.push_back({x,y,z});
121  }
122  if (t==0 || t==npt) break;
123  }
124  }
125  return v;
126  }
127 };
128 
129 // This functor is one inside the molecule, 1/2 on the surface, and zero
130 // exterior to the molecule.
132  , public madness::FunctionFunctorInterface<double,3> {
133 public:
135  const std::vector<double> atomic_radii,
136  const std::vector<madness::coord_3d> atomic_coords)
138  {}
139 
140  virtual double operator()(const madness::coord_3d& r) const {
141  double value = 1.0;
142  for (int i=0; i<natom; i++) {
143  value *= atomic_cmask(r,i);
144  }
145  return 1.0 - value;
146  }
147 
148  std::vector< madness::Vector<double,3> > special_points() const {
150  }
151 };
152 
153 // This functor is zero inside the molecule, 1/2 on the surface, and one
154 // exterior to the molecule.
156  , public madness::FunctionFunctorInterface<double,3> {
157 public:
159  const std::vector<double> atomic_radii,
160  const std::vector<madness::coord_3d> atomic_coords)
162  {}
163 
164  virtual double operator()(const madness::coord_3d& r) const {
165  double value = 1.0;
166  for (int i=0; i<natom; i++) {
167  value *= atomic_cmask(r,i);
168  }
169  return value;
170  }
171 
172  std::vector< madness::Vector<double,3> > special_points() const {
174  }
175 
176 };
177 
178 /// Switches between \a positive values \c Vint and \c Vext with special log derivative
179 
180 
181 /// Switches between \a positive values \c Vint on the interior and \c Vext on the interior.
182 /// It has value
183 /// \f[
184 /// V_(r,\sigma) = \exp \left( \log V_{int} C(r,\sigma) + \log V_{ext} (1 - C(r,\sigma) \right)
185 /// = V_{ext} \exp \left( \log \frac{V_{int}}{V_{ext}} C(r,\sigma) \right)
186 /// = V_{int} \exp \left( \log \frac{V_{ext}}{V_{int}} \overline{C}(r,\sigma) \right)
187 /// \f]
188 /// where \f$ C(r,\sigma) \f$ is the regular volume mask provided by MolecularVolumeMask,
189 /// and \f$ \overline{C} \f$ is its complement.
190 /// Its log-derivative is precisely located in the surface with value
191 /// \f[
192 /// \nabla \log V = \frac{\nabla V}{V} = \log \frac{V_{int}}{V_{ext}} \nabla C
193 /// \f]
194 /// with \f$ \nabla C \f$ already being computed by MolecularVolumeMaskGrad.
195 /// The advantage of this is that if \f$ \| V_{ext} - V_{int} \| \f$
196 /// is big, the log derivative of the regular volume mask (\f$ C \f$) is
197 /// displaced from the surface (perhaps by multiple values of \f$
198 /// \sigma \f$). This leads to slow convergence w.r.t \f$ \sigma \f$ and
199 /// potential inaccuracies depending how the numerical representation is computed.
200 /// The surface charge
201 /// in dielectric problems is controlled by \f$ \nabla \log \epsilon
202 /// \f$ and hence the dielectric should employ this form of the switch.
205  const double Vint, Vext, fac;
206 public:
208  double Vint,
209  double Vext,
210  const std::vector<double> atomic_radii,
211  const std::vector<madness::coord_3d> atomic_coords)
212  : cmask(sigma, atomic_radii, atomic_coords)
213  , Vint(Vint)
214  , Vext(Vext)
215  , fac(log(Vext/Vint))
216  {
217  if (Vint <= 0 || Vext <= 0) throw "Only works for positive values";
218  }
219 
220  virtual double operator()(const madness::coord_3d& r) const {
221  double c = cmask(r);
222  if (c == 0.0) return Vint;
223  else if (c == 1.0) return Vext;
224  else return Vint * exp(fac * c);
225  }
226 
227  std::vector< madness::Vector<double,3> > special_points() const {
228  return cmask.special_points();
229  }
230 };
231 
232 /// Computes the reciprocal of MolecularVolumeExponentialSwitch
235 public:
237  double Vint,
238  double Vext,
239  const std::vector<double> atomic_radii,
240  const std::vector<madness::coord_3d> atomic_coords)
241  : s(sigma, Vint, Vext, atomic_radii, atomic_coords)
242  {}
243 
244  virtual double operator()(const madness::coord_3d& r) const {
245  return 1.0/s(r);
246  }
247 };
248 
249 
250 // This functor is a shell that limits to a delta function in the
251 // molecular surface and integrates to the molecular surface area.
253  , public madness::FunctionFunctorInterface<double,3> {
254 public:
256  const std::vector<double> atomic_radii,
257  const std::vector<madness::coord_3d> atomic_coords)
259  {}
260 
261  virtual double operator()(const madness::coord_3d& r) const {
263  return sqrt(grad[0]*grad[0] + grad[1]*grad[1] + grad[2]*grad[2]);
264  }
265 
266  std::vector< madness::Vector<double,3> > special_points() const {
268  }
269 };
270 
271 // Evaluates component i (0=x, 1=y, 2=z) of the gradient of the mask
273  , public madness::FunctionFunctorInterface<double,3> {
274  const int i;
275 public:
277  const std::vector<double> atomic_radii,
278  const std::vector<madness::coord_3d> atomic_coords,
279  int i)
281  , i(i)
282  {}
283 
284  virtual double operator()(const madness::coord_3d& r) const {
286  return grad[i];
287  }
288 
289  std::vector< madness::Vector<double,3> > special_points() const {
291  }
292 };
293 
294 /// Returns the requested component of the derivative of the log of MolecularVolumeExponentialSwitch
297  const double fac;
298 public:
300  double Vint,
301  double Vext,
302  const std::vector<double> atomic_radii,
303  const std::vector<madness::coord_3d> atomic_coords,
304  int i)
305  : g(sigma, atomic_radii, atomic_coords, i)
306  , fac(log(Vint/Vext))
307  {
308  if (Vint <= 0 || Vext <= 0) throw "Only works for positive values";
309  }
310 
311  virtual double operator()(const madness::coord_3d& r) const {
312  return fac * g(r);
313  }
314 
315  std::vector< madness::Vector<double,3> > special_points() const {
316  return g.special_points();
317  }
318 };
319 
320 #endif
Definition: molecularmask.h:18
const std::vector< double > atomic_radii
Definition: molecularmask.h:21
MolecularMaskBase(double sigma, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords)
Definition: molecularmask.h:95
madness::coord_3d grad_sdf(const madness::coord_3d &r, const madness::coord_3d &center) const
Definition: molecularmask.h:31
double sdf(const madness::coord_3d &r, const madness::coord_3d &center, double R) const
Definition: molecularmask.h:26
madness::coord_3d grad_atomic_mask(const madness::coord_3d &r, unsigned int i) const
Definition: molecularmask.h:68
double cmask(double s) const
Definition: molecularmask.h:43
std::vector< madness::Vector< double, 3 > > special_points() const
Definition: molecularmask.h:106
double atomic_cmask(const madness::coord_3d &r, unsigned int i) const
Definition: molecularmask.h:62
const double sigma
Definition: molecularmask.h:20
const std::vector< madness::coord_3d > atomic_coords
Definition: molecularmask.h:22
double atomic_mask(const madness::coord_3d &r, unsigned int i) const
Definition: molecularmask.h:55
const int natom
Definition: molecularmask.h:23
madness::coord_3d gradient(const madness::coord_3d &r) const
Definition: molecularmask.h:74
double mask(double s) const
Definition: molecularmask.h:36
double dmask(double s) const
Definition: molecularmask.h:48
Definition: molecularmask.h:253
virtual double operator()(const madness::coord_3d &r) const
Definition: molecularmask.h:261
MolecularSurface(double sigma, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords)
Definition: molecularmask.h:255
std::vector< madness::Vector< double, 3 > > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: molecularmask.h:266
Definition: molecularmask.h:156
MolecularVolumeComplementMask(double sigma, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords)
Definition: molecularmask.h:158
virtual double operator()(const madness::coord_3d &r) const
Definition: molecularmask.h:164
std::vector< madness::Vector< double, 3 > > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: molecularmask.h:172
Returns the requested component of the derivative of the log of MolecularVolumeExponentialSwitch.
Definition: molecularmask.h:295
virtual double operator()(const madness::coord_3d &r) const
Definition: molecularmask.h:311
const double fac
Definition: molecularmask.h:297
MolecularVolumeExponentialSwitchLogGrad(double sigma, double Vint, double Vext, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords, int i)
Definition: molecularmask.h:299
std::vector< madness::Vector< double, 3 > > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: molecularmask.h:315
const MolecularVolumeMaskGrad g
Definition: molecularmask.h:296
Computes the reciprocal of MolecularVolumeExponentialSwitch.
Definition: molecularmask.h:233
virtual double operator()(const madness::coord_3d &r) const
Definition: molecularmask.h:244
MolecularVolumeExponentialSwitchReciprocal(double sigma, double Vint, double Vext, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords)
Definition: molecularmask.h:236
const MolecularVolumeExponentialSwitch s
Definition: molecularmask.h:234
Switches between positive values Vint and Vext with special log derivative.
Definition: molecularmask.h:203
const double fac
Definition: molecularmask.h:205
virtual double operator()(const madness::coord_3d &r) const
Definition: molecularmask.h:220
MolecularVolumeExponentialSwitch(double sigma, double Vint, double Vext, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords)
Definition: molecularmask.h:207
std::vector< madness::Vector< double, 3 > > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: molecularmask.h:227
const MolecularVolumeComplementMask cmask
Definition: molecularmask.h:204
const double Vint
Definition: molecularmask.h:205
const double Vext
Definition: molecularmask.h:205
Definition: molecularmask.h:273
virtual double operator()(const madness::coord_3d &r) const
Definition: molecularmask.h:284
MolecularVolumeMaskGrad(double sigma, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords, int i)
Definition: molecularmask.h:276
std::vector< madness::Vector< double, 3 > > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: molecularmask.h:289
const int i
Definition: molecularmask.h:274
Definition: molecularmask.h:132
std::vector< madness::Vector< double, 3 > > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: molecularmask.h:148
virtual double operator()(const madness::coord_3d &r) const
Definition: molecularmask.h:140
MolecularVolumeMask(double sigma, const std::vector< double > atomic_radii, const std::vector< madness::coord_3d > atomic_coords)
Definition: molecularmask.h:134
Abstract base class interface required for functors used as input to Functions.
Definition: function_interface.h:68
Defines common mathematical and physical constants.
static const double R
Definition: csqrt.cc:46
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition: derivatives.cc:72
const double sigma
Definition: dielectric.cc:185
const double m
Definition: gfit.cc:199
static const double v
Definition: hatom_sf_dirac.cc:20
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
double distance(const madness::coord_3d &a, const madness::coord_3d &b)
Definition: molecularmask.h:10
Main include file for MADNESS and defines Function interface.
const double pi
Mathematical constant .
Definition: constants.h:48
std::vector< Function< T, NDIM > > grad(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
shorthand gradient operator
Definition: vmra.h:1818
Vector< double, 3 > coord_3d
Definition: funcplot.h:1042
static const double b
Definition: nonlinschro.cc:119
static const double a
Definition: nonlinschro.cc:118
static const double c
Definition: relops.cc:10
void e()
Definition: test_sig.cc:75