1 //
2 // bem.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #include <stdio.h>
29 #include <util/misc/math.h>
30 #include <util/misc/formio.h>
31 #include <util/misc/timer.h>
32 #include <math/scmat/matrix.h>
33 #include <math/scmat/vector3.h>
34 #include <math/scmat/local.h>
35 #include <chemistry/solvent/bem.h>
36
37 using namespace std;
38 using namespace sc;
39
40 static ClassDesc BEMSolvent_cd(
41 typeid(BEMSolvent),"BEMSolvent",1,"public DescribedClass",
42 0, create<BEMSolvent>, 0);
43
BEMSolvent(const Ref<KeyVal> & keyval)44 BEMSolvent::BEMSolvent(const Ref<KeyVal>& keyval)
45 {
46 vertex_area_ = 0;
47
48 matrixkit_ = new LocalSCMatrixKit;
49
50 debug_ = keyval->intvalue("debug");
51
52 solute_ << keyval->describedclassvalue("solute");
53
54 solvent_ << keyval->describedclassvalue("solvent");
55 // Use the aug-cc-pVQZ MP2 optimum geometry for H2O as default
56 if (solvent_.null()) {
57 solvent_ = new Molecule;
58 solvent_->add_atom(8, 0.0000000000, 0.0000000000, -0.1265941233);
59 solvent_->add_atom(1, 0.0000000000, 1.4304840085, 0.9856159541);
60 solvent_->add_atom(1, 0.0000000000, -1.4304840085, 0.9856159541);
61 }
62
63 solvent_density_ = keyval->doublevalue("solvent_density");
64 // use as default the number density of water in au^-3, T=25 C, P=101325 Pa
65 if (keyval->error() != KeyVal::OK) solvent_density_ = 0.004938887;
66
67 surf_ << keyval->describedclassvalue("surface");
68
69 dielectric_constant_ = keyval->doublevalue("dielectric_constant");
70 if (keyval->error() != KeyVal::OK) dielectric_constant_ = 78.0;
71
72 grp_ = MessageGrp::get_default_messagegrp();
73 }
74
~BEMSolvent()75 BEMSolvent::~BEMSolvent()
76 {
77 }
78
79 double**
alloc_array(int n,int m)80 BEMSolvent::alloc_array(int n, int m)
81 {
82 double ** result = new double*[n];
83 result[0] = new double[n*m];
84 for (int i=1; i<n; i++) {
85 result[i] = &result[i-1][m];
86 }
87 return result;
88 }
89
90 void
free_array(double ** array)91 BEMSolvent::free_array(double** array)
92 {
93 if (!array) return;
94 delete[] array[0];
95 delete[] array;
96 }
97
98 void
charge_positions(double ** pos)99 BEMSolvent::charge_positions(double**pos)
100 {
101 int i,j;
102 int n = ncharge();
103 for (i=0; i<n; i++) {
104 const SCVector3& p = surf_->vertex(i)->point();
105 for (j=0; j<3; j++) {
106 pos[i][j] = p[j];
107 }
108 }
109 }
110
111 void
normals(double ** norms)112 BEMSolvent::normals(double**norms)
113 {
114 int i,j;
115 int n = ncharge();
116 for (i=0; i<n; i++) {
117 const SCVector3& p = surf_->vertex(i)->normal();
118 for (j=0; j<3; j++) {
119 norms[i][j] = p[j];
120 }
121 }
122 }
123
124 void
init()125 BEMSolvent::init()
126 {
127 surf_->clear();
128 surf_->init();
129 system_matrix_i_ = 0;
130
131 f_ = (1.0-dielectric_constant_)/(2.0*M_PI*(1.0+dielectric_constant_));
132
133 if (vertex_area_) delete[] vertex_area_;
134 vertex_area_ = new double[ncharge()];
135 for (int i=0; i<ncharge(); i++) vertex_area_[i] = 0.0;
136 TriangulatedSurfaceIntegrator triint(surf_.pointer());
137 for (triint = 0; triint.update(); triint++) {
138 int j0 = triint.vertex_number(0);
139 int j1 = triint.vertex_number(1);
140 int j2 = triint.vertex_number(2);
141 double r = triint.r();
142 double s = triint.s();
143 double dA = triint.w();
144 vertex_area_[j0] += dA * (1 - r - s);
145 vertex_area_[j1] += dA * r;
146 vertex_area_[j2] += dA * s;
147 }
148 }
149
150 void
done(int clear_surface)151 BEMSolvent::done(int clear_surface)
152 {
153 if (clear_surface) surf_->clear();
154 system_matrix_i_ = 0;
155
156 if (vertex_area_) delete[] vertex_area_;
157 vertex_area_ = 0;
158 }
159
160 void
charges_to_surface_charge_density(double * charges)161 BEMSolvent::charges_to_surface_charge_density(double *charges)
162 {
163 for (int i=0; i<ncharge(); i++) charges[i] /= vertex_area_[i];
164 }
165
166 void
surface_charge_density_to_charges(double * charges)167 BEMSolvent::surface_charge_density_to_charges(double *charges)
168 {
169 for (int i=0; i<ncharge(); i++) charges[i] *= vertex_area_[i];
170 }
171
172 double
polarization_charge(double * charges)173 BEMSolvent::polarization_charge(double *charges)
174 {
175 double charge = 0.0;
176 int n = ncharge();
177 for (int i=0; i<n; i++) charge += charges[i];
178 return charge;
179 }
180
181 // the passed enclosed_charge is determined by the called and
182 // might different from the enclosed charge computed by Gauss's
183 // law, which is stored as computed_enclosed_charge_
184 void
normalize_charge(double enclosed_charge,double * charges)185 BEMSolvent::normalize_charge(double enclosed_charge, double* charges)
186 {
187 int i;
188 double expected_charge = enclosed_charge
189 * (1.0/dielectric_constant_ - 1.0);
190 double charge = 0.0;
191 double charge_pos = 0.0;
192 double charge_neg = 0.0;
193 int n = ncharge();
194 for (i=0; i<n; i++) {
195 charge += charges[i];
196 if (charges[i] > 0.0) charge_pos += charges[i];
197 else charge_neg += charges[i];
198 }
199
200 double scale_pos = 1.0;
201 double scale_neg = 1.0;
202 if (charge_pos > 1.0e-4 && charge_neg < -1.0e-4) {
203 scale_pos += (expected_charge-charge)/(2.0*charge_pos);
204 scale_neg += (expected_charge-charge)/(2.0*charge_neg);
205 }
206 else if (charge_pos > 1.0e-4) {
207 scale_pos += (expected_charge-charge)/charge_pos;
208 }
209 else if (charge_neg < -1.0e-4) {
210 scale_neg += (expected_charge-charge)/charge_neg;
211 }
212
213 double new_charge = 0.0;
214 for (i=0; i<n; i++) {
215 if (charges[i] > 0.0) charges[i] *= scale_pos;
216 else charges[i] *= scale_neg;
217 new_charge += charges[i];
218 }
219
220 if (fabs(new_charge - expected_charge) > 1.0e-3) {
221 ExEnv::outn() << "BEMSolvent:normalize_charge: failed:" << endl
222 << "new_charge = " << new_charge << endl
223 << "expected_charge = " << expected_charge << endl;
224 abort();
225 }
226
227 if (debug_) {
228 ExEnv::out0() << indent
229 << "BEMSolvent:normalize_charge:"
230 << endl << indent
231 << scprintf(" integrated surface charge = %20.15f", charge)
232 << endl << indent
233 << scprintf(" expected surface charge = %20.15f", expected_charge)
234 << endl;
235 }
236 }
237
238 void
init_system_matrix()239 BEMSolvent::init_system_matrix()
240 {
241 int i, j;
242 int n = ncharge();
243
244 RefSCDimension d = new SCDimension(n);
245 RefSCMatrix system_matrix(d,d,matrixkit());
246 system_matrix.assign(0.0);
247
248 tim_enter("precomp");
249 // precompute some arrays
250 TriangulatedSurfaceIntegrator triint(surf_.pointer());
251 int n_integration_points = triint.n();
252 SCVector3 *surfpv = new SCVector3[n_integration_points];
253 double *rfdA = new double[n_integration_points];
254 double *sfdA = new double[n_integration_points];
255 double *rsfdA = new double[n_integration_points];
256 int *j0 = new int[n_integration_points];
257 int *j1 = new int[n_integration_points];
258 int *j2 = new int[n_integration_points];
259 for (triint=0, i=0; i<n_integration_points&&triint.update(); i++,triint++) {
260 surfpv[i] = triint.current()->point();
261 j0[i] = triint.vertex_number(0);
262 j1[i] = triint.vertex_number(1);
263 j2[i] = triint.vertex_number(2);
264 double r = triint.r();
265 double s = triint.s();
266 double rs = 1 - r - s;
267 double dA = triint.w();
268 double fdA = - f_ * dA;
269 rfdA[i] = r * fdA;
270 sfdA[i] = s * fdA;
271 rsfdA[i] = rs * fdA;
272 }
273 tim_exit("precomp");
274
275 tim_enter("sysmat");
276 double *sysmati = new double[n];
277 RefSCVector vsysmati(system_matrix->rowdim(),system_matrix->kit());
278 // loop thru all the vertices
279 for (i = 0; i<n; i++) {
280 memset(sysmati,0,sizeof(double)*n);
281 Ref<Vertex> v = surf_->vertex(i);
282 const SCVector3& pv = v->point();
283 const SCVector3& nv = v->normal();
284 // integrate over the surface
285 for (j = 0; j < n_integration_points; j++) {
286 SCVector3 diff(pv - surfpv[j]);
287 double normal_component = diff.dot(nv);
288 double diff2 = diff.dot(diff);
289 if (diff2 <= 1.0e-8) {
290 // The self term must not be included here. This
291 // case shouldn't occur for the usual integrators
292 // so abort.
293 ExEnv::errn() << "BEMSolvent: integrator gave the self term" << endl;
294 abort();
295 }
296 double denom = diff2*sqrt(diff2);
297 double common_factor = normal_component/denom;
298 sysmati[j0[j]] += common_factor * rsfdA[j];
299 sysmati[j1[j]] += common_factor * rfdA[j];
300 sysmati[j2[j]] += common_factor * sfdA[j];
301 }
302 vsysmati->assign(sysmati);
303 system_matrix->assign_row(vsysmati,i);
304 }
305 tim_exit("sysmat");
306
307 delete[] surfpv;
308 delete[] rfdA;
309 delete[] sfdA;
310 delete[] rsfdA;
311 delete[] j0;
312 delete[] j1;
313 delete[] j2;
314 delete[] sysmati;
315
316 tim_enter("AV");
317 double A = 0.0;
318 double V = 0.0;
319 for (triint = 0; triint.update(); triint++) {
320 V += triint.weight()*triint.dA()[2]*triint.current()->point()[2];
321 A += triint.w();
322 }
323 area_ = A;
324 volume_ = V;
325 tim_exit("AV");
326
327 ExEnv::out0() << indent
328 << scprintf("Solvent Accessible Surface:") << endl
329 << indent
330 << scprintf(" Area = %15.10f ", A)
331 << scprintf("Volume = %15.10f ", V)
332 << scprintf("Nvertex = %3d", n) << endl;
333
334 // Add I to the system matrix.
335 system_matrix->shift_diagonal(1.0);
336
337 //system_matrix->print("System Matrix");
338
339 tim_enter("inv");
340 system_matrix->invert_this();
341 system_matrix_i_ = system_matrix;
342 tim_exit("inv");
343
344 //system_matrix_i_->print("System Matrix Inverse");
345 }
346
347 void
compute_charges(double * efield_dot_normals,double * charges)348 BEMSolvent::compute_charges(double* efield_dot_normals, double* charges)
349 {
350 if (system_matrix_i_.null()) {
351 tim_enter("sysmat");
352 init_system_matrix();
353 tim_exit("sysmat");
354 }
355
356 tim_enter("qenq");
357 double efield_dot_normal = 0.0;
358 int n = ncharge();
359 for (int i=0; i<n; i++)
360 efield_dot_normal += efield_dot_normals[i] * vertex_area_[i];
361 tim_exit("qenq");
362
363 computed_enclosed_charge_ = efield_dot_normal/(4.0*M_PI);
364
365 if (debug_) {
366 double computed_expected_charge = computed_enclosed_charge_
367 * (1.0/dielectric_constant_ - 1.0);
368
369 ExEnv::out0() << indent
370 << scprintf("BEMSolvent:compute_charges: encl q = %20.15f",
371 computed_enclosed_charge_)
372 << endl << indent
373 << scprintf("BEMSolvent:compute_charges: exp surface q = %20.15f",
374 computed_expected_charge) << endl;
375 }
376
377 tim_enter("scomp");
378 RefSCVector edotn(system_matrix_i_.coldim(),matrixkit());
379 edotn.assign(efield_dot_normals);
380 //edotn.print("E dot normals");
381 edotn.scale(f_);
382 RefSCVector chrg = system_matrix_i_ * edotn;
383 //chrg.print("Charges");
384 chrg.convert(charges);
385 tim_exit("scomp");
386
387 tim_enter("stoq");
388 surface_charge_density_to_charges(charges);
389 tim_exit("stoq");
390 }
391
392 double
nuclear_charge_interaction_energy(double * nuclear_charge,double ** charge_positions,double * charge)393 BEMSolvent::nuclear_charge_interaction_energy(double *nuclear_charge,
394 double** charge_positions,
395 double* charge)
396 {
397 double energy = 0.0;
398 int natom = solute_->natom();
399 for (int i=0; i<natom; i++) {
400 for (int j=0; j<ncharge(); j++) {
401 double r2 = 0.0;
402 for (int k=0; k<3; k++) {
403 double r = charge_positions[j][k] - solute_->r(i,k);
404 r2 += r*r;
405 }
406 energy += nuclear_charge[i] * charge[j] / sqrt(r2);
407 }
408 }
409 return energy;
410 }
411
412 double
nuclear_interaction_energy(double ** charge_positions,double * charge)413 BEMSolvent::nuclear_interaction_energy(double** charge_positions,
414 double* charge)
415 {
416 double energy = 0.0;
417 int natom = solute_->natom();
418 for (int i=0; i<natom; i++) {
419 for (int j=0; j<ncharge(); j++) {
420 double r2 = 0.0;
421 for (int k=0; k<3; k++) {
422 double r = charge_positions[j][k] - solute_->r(i,k);
423 r2 += r*r;
424 }
425 energy += double(solute_->Z(i)) * charge[j] / sqrt(r2);
426 }
427 }
428 return energy;
429 }
430
431 double
self_interaction_energy(double ** charge_positions,double * charge)432 BEMSolvent::self_interaction_energy(double** charge_positions,
433 double* charge)
434 {
435 int i,j;
436
437 charges_to_surface_charge_density(charge);
438
439 TriangulatedSurfaceIntegrator triint(surf_.pointer());
440 int n_integration_points = triint.n();
441 SCVector3 *points = new SCVector3[n_integration_points];
442 double *charges = new double[n_integration_points];
443
444 double energy = 0.0;
445 for (triint=0, i=0; i<n_integration_points&&triint.update(); i++,triint++) {
446 points[i] = triint.current()->point();
447 int v0 = triint.vertex_number(0);
448 int v1 = triint.vertex_number(1);
449 int v2 = triint.vertex_number(2);
450 double r = triint.r();
451 double s = triint.s();
452 double rs = 1.0 - r - s;
453 double dA = triint.w();
454 charges[i] = (charge[v0]*rs + charge[v1]*r + charge[v2]*s)*dA;
455 energy += 0.0; // is this good enough for the self term?
456 }
457 for (i=0; i<n_integration_points; i++) {
458 double chargesi = charges[i];
459 SCVector3 pointsi(points[i]);
460 for (j = 0; j<i; j++) {
461 energy += chargesi*charges[j]/pointsi.dist(points[j]);
462 }
463 }
464
465 delete[] points;
466 delete[] charges;
467
468 surface_charge_density_to_charges(charge);
469
470 return energy;
471 }
472
473 /////////////////////////////////////////////////////////////////////////////
474
475 // Local Variables:
476 // mode: c++
477 // c-file-style: "CLJ"
478 // End:
479