APBS  1.4.1
vopot.c
Go to the documentation of this file.
1 
57 #include "vopot.h"
58 
59 VEMBED(rcsid="$Id$")
60 
61 /* ///////////////////////////////////////////////////////////////////////////
62 // Routine: Vopot_ctor
63 // Author: Nathan Baker
65 VPUBLIC Vopot* Vopot_ctor(Vmgrid *mgrid, Vpbe *pbe, Vbcfl bcfl) {
66 
67  Vopot *thee = VNULL;
68 
69  thee = Vmem_malloc(VNULL, 1, sizeof(Vopot));
70  VASSERT(thee != VNULL);
71  VASSERT(Vopot_ctor2(thee, mgrid, pbe, bcfl));
72 
73  return thee;
74 }
75 
76 /* ///////////////////////////////////////////////////////////////////////////
77 // Routine: Vopot_ctor2
78 // Author: Nathan Baker
80 VPUBLIC int Vopot_ctor2(Vopot *thee, Vmgrid *mgrid, Vpbe *pbe, Vbcfl bcfl) {
81 
82  if (thee == VNULL) return 0;
83  thee->bcfl = bcfl;
84  thee->mgrid = mgrid;
85  thee->pbe = pbe;
86 
87  return 1;
88 }
89 
90 /* ///////////////////////////////////////////////////////////////////////////
91 // Routine: Vopot_dtor
92 // Author: Nathan Baker
94 VPUBLIC void Vopot_dtor(Vopot **thee) {
95 
96  if ((*thee) != VNULL) {
97  Vopot_dtor2(*thee);
98  Vmem_free(VNULL, 1, sizeof(Vopot), (void **)thee);
99  (*thee) = VNULL;
100  }
101 }
102 
103 /* ///////////////////////////////////////////////////////////////////////////
104 // Routine: Vopot_dtor2
105 // Author: Nathan Baker
107 VPUBLIC void Vopot_dtor2(Vopot *thee) { return; }
108 
109 /* ///////////////////////////////////////////////////////////////////////////
110 // Routine: Vopot_pot
111 // Author: Nathan Baker
113 #define IJK(i,j,k) (((k)*(nx)*(ny))+((j)*(nx))+(i))
114 VPUBLIC int Vopot_pot(Vopot *thee, double pt[3], double *value) {
115 
116  Vatom *atom;
117  int i, iatom;
118  double u, T, charge, eps_w, xkappa, dist, size, val, *position;
119  Valist *alist;
120 
121  VASSERT(thee != VNULL);
122 
123  eps_w = Vpbe_getSolventDiel(thee->pbe);
124  xkappa = (1.0e10)*Vpbe_getXkappa(thee->pbe);
125  T = Vpbe_getTemperature(thee->pbe);
126  alist = Vpbe_getValist(thee->pbe);
127 
128  u = 0;
129 
130  /* See if we're on the mesh */
131  if (Vmgrid_value(thee->mgrid, pt, &u)) {
132 
133  *value = u;
134 
135  } else {
136 
137  switch (thee->bcfl) {
138 
139  case BCFL_ZERO:
140  u = 0;
141  break;
142 
143  case BCFL_SDH:
144  size = (1.0e-10)*Vpbe_getSoluteRadius(thee->pbe);
145  position = Vpbe_getSoluteCenter(thee->pbe);
146  charge = Vunit_ec*Vpbe_getSoluteCharge(thee->pbe);
147  dist = 0;
148  for (i=0; i<3; i++)
149  dist += VSQR(position[i] - pt[i]);
150  dist = (1.0e-10)*VSQRT(dist);
151  val = (charge)/(4*VPI*Vunit_eps0*eps_w*dist);
152  if (xkappa != 0.0)
153  val = val*(exp(-xkappa*(dist-size))/(1+xkappa*size));
154  val = val*Vunit_ec/(Vunit_kb*T);
155  u = val;
156  break;
157 
158  case BCFL_MDH:
159  u = 0;
160  for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++) {
161  atom = Valist_getAtom(alist, iatom);
162  position = Vatom_getPosition(atom);
163  charge = Vunit_ec*Vatom_getCharge(atom);
164  size = (1e-10)*Vatom_getRadius(atom);
165  dist = 0;
166  for (i=0; i<3; i++)
167  dist += VSQR(position[i] - pt[i]);
168  dist = (1.0e-10)*VSQRT(dist);
169  val = (charge)/(4*VPI*Vunit_eps0*eps_w*dist);
170  if (xkappa != 0.0)
171  val = val*(exp(-xkappa*(dist-size))/(1+xkappa*size));
172  val = val*Vunit_ec/(Vunit_kb*T);
173  u = u + val;
174  }
175  break;
176 
177  case BCFL_UNUSED:
178  Vnm_print(2, "Vopot_pot: Invalid bcfl flag (%d)!\n",
179  thee->bcfl);
180  return 0;
181 
182  case BCFL_FOCUS:
183  Vnm_print(2, "Vopot_pot: Invalid bcfl flag (%d)!\n",
184  thee->bcfl);
185  return 0;
186 
187  default:
188  Vnm_print(2, "Vopot_pot: Bogus thee->bcfl flag (%d)!\n",
189  thee->bcfl);
190  return 0;
191  break;
192  }
193 
194  *value = u;
195 
196  }
197 
198  return 1;
199 
200 }
201 
202 /* ///////////////////////////////////////////////////////////////////////////
203 // Routine: Vopot_curvature
204 //
205 // Notes: cflag=0 ==> Reduced Maximal Curvature
206 // cflag=1 ==> Mean Curvature (Laplace)
207 // cflag=2 ==> Gauss Curvature
208 // cflag=3 ==> True Maximal Curvature
209 // If we are off the grid, we can still evaluate the Laplacian; assuming, we
210 // are away from the molecular surface, it is simply equal to the DH factor.
211 //
212 // Authors: Nathan Baker
214 VPUBLIC int Vopot_curvature(Vopot *thee, double pt[3], int cflag,
215  double *value) {
216 
217  Vatom *atom;
218  int i, iatom;
219  double u, T, charge, eps_w, xkappa, dist, size, val, *position, zkappa2;
220  Valist *alist;
221 
222  VASSERT(thee != VNULL);
223 
224  eps_w = Vpbe_getSolventDiel(thee->pbe);
225  xkappa = (1.0e10)*Vpbe_getXkappa(thee->pbe);
226  zkappa2 = Vpbe_getZkappa2(thee->pbe);
227  T = Vpbe_getTemperature(thee->pbe);
228  alist = Vpbe_getValist(thee->pbe);
229 
230  u = 0;
231 
232  if (Vmgrid_curvature(thee->mgrid, pt, cflag, value)) return 1;
233  else if (cflag != 1) {
234  Vnm_print(2, "Vopot_curvature: Off mesh!\n");
235  return 1;
236  } else {
237 
238  switch (thee->bcfl) {
239 
240  case BCFL_ZERO:
241  u = 0;
242  break;
243 
244  case BCFL_SDH:
245  size = (1.0e-10)*Vpbe_getSoluteRadius(thee->pbe);
246  position = Vpbe_getSoluteCenter(thee->pbe);
247  charge = Vunit_ec*Vpbe_getSoluteCharge(thee->pbe);
248  dist = 0;
249  for (i=0; i<3; i++)
250  dist += VSQR(position[i] - pt[i]);
251  dist = (1.0e-10)*VSQRT(dist);
252  if (xkappa != 0.0)
253  u = zkappa2*(exp(-xkappa*(dist-size))/(1+xkappa*size));
254  break;
255 
256  case BCFL_MDH:
257  u = 0;
258  for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++) {
259  atom = Valist_getAtom(alist, iatom);
260  position = Vatom_getPosition(atom);
261  charge = Vunit_ec*Vatom_getCharge(atom);
262  size = (1e-10)*Vatom_getRadius(atom);
263  dist = 0;
264  for (i=0; i<3; i++)
265  dist += VSQR(position[i] - pt[i]);
266  dist = (1.0e-10)*VSQRT(dist);
267  if (xkappa != 0.0)
268  val = zkappa2*(exp(-xkappa*(dist-size))/(1+xkappa*size));
269  u = u + val;
270  }
271  break;
272 
273  case BCFL_UNUSED:
274  Vnm_print(2, "Vopot_pot: Invlid bcfl (%d)!\n", thee->bcfl);
275  return 0;
276 
277  case BCFL_FOCUS:
278  Vnm_print(2, "Vopot_pot: Invlid bcfl (%d)!\n", thee->bcfl);
279  return 0;
280 
281  default:
282  Vnm_print(2, "Vopot_pot: Bogus thee->bcfl flag (%d)!\n",
283  thee->bcfl);
284  return 0;
285  break;
286  }
287 
288  *value = u;
289  }
290 
291  return 1;
292 
293 }
294 
295 /* ///////////////////////////////////////////////////////////////////////////
296 // Routine: Vopot_gradient
297 //
298 // Authors: Nathan Baker
300 VPUBLIC int Vopot_gradient(Vopot *thee, double pt[3], double grad[3]) {
301 
302  Vatom *atom;
303  int iatom;
304  double T, charge, eps_w, xkappa, size, val, *position;
305  double dx, dy, dz, dist;
306  Valist *alist;
307 
308  VASSERT(thee != VNULL);
309 
310  eps_w = Vpbe_getSolventDiel(thee->pbe);
311  xkappa = (1.0e10)*Vpbe_getXkappa(thee->pbe);
312  T = Vpbe_getTemperature(thee->pbe);
313  alist = Vpbe_getValist(thee->pbe);
314 
315 
316  if (!Vmgrid_gradient(thee->mgrid, pt, grad)) {
317 
318  switch (thee->bcfl) {
319 
320  case BCFL_ZERO:
321  grad[0] = 0.0;
322  grad[1] = 0.0;
323  grad[2] = 0.0;
324  break;
325 
326  case BCFL_SDH:
327  grad[0] = 0.0;
328  grad[1] = 0.0;
329  grad[2] = 0.0;
330  size = (1.0e-10)*Vpbe_getSoluteRadius(thee->pbe);
331  position = Vpbe_getSoluteCenter(thee->pbe);
332  charge = Vunit_ec*Vpbe_getSoluteCharge(thee->pbe);
333  dx = position[0] - pt[0];
334  dy = position[1] - pt[1];
335  dz = position[2] - pt[2];
336  dist = VSQR(dx) + VSQR(dy) + VSQR(dz);
337  dist = (1.0e-10)*VSQRT(dist);
338  val = (charge)/(4*VPI*Vunit_eps0*eps_w);
339  if (xkappa != 0.0)
340  val = val*(exp(-xkappa*(dist-size))/(1+xkappa*size));
341  val = val*Vunit_ec/(Vunit_kb*T);
342  grad[0] = val*dx/dist*(-1.0/dist/dist + xkappa/dist);
343  grad[1] = val*dy/dist*(-1.0/dist/dist + xkappa/dist);
344  grad[2] = val*dz/dist*(-1.0/dist/dist + xkappa/dist);
345  break;
346 
347  case BCFL_MDH:
348  grad[0] = 0.0;
349  grad[1] = 0.0;
350  grad[2] = 0.0;
351  for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++) {
352  atom = Valist_getAtom(alist, iatom);
353  position = Vatom_getPosition(atom);
354  charge = Vunit_ec*Vatom_getCharge(atom);
355  size = (1e-10)*Vatom_getRadius(atom);
356  dx = position[0] - pt[0];
357  dy = position[1] - pt[1];
358  dz = position[2] - pt[2];
359  dist = VSQR(dx) + VSQR(dy) + VSQR(dz);
360  dist = (1.0e-10)*VSQRT(dist);
361  val = (charge)/(4*VPI*Vunit_eps0*eps_w);
362  if (xkappa != 0.0)
363  val = val*(exp(-xkappa*(dist-size))/(1+xkappa*size));
364  val = val*Vunit_ec/(Vunit_kb*T);
365  grad[0] += (val*dx/dist*(-1.0/dist/dist + xkappa/dist));
366  grad[1] += (val*dy/dist*(-1.0/dist/dist + xkappa/dist));
367  grad[2] += (val*dz/dist*(-1.0/dist/dist + xkappa/dist));
368  }
369  break;
370 
371  case BCFL_UNUSED:
372  Vnm_print(2, "Vopot: Invalid bcfl (%d)!\n", thee->bcfl);
373  return 0;
374 
375  case BCFL_FOCUS:
376  Vnm_print(2, "Vopot: Invalid bcfl (%d)!\n", thee->bcfl);
377  return 0;
378 
379  default:
380  Vnm_print(2, "Vopot_pot: Bogus thee->bcfl flag (%d)!\n",
381  thee->bcfl);
382  return 0;
383  break;
384  }
385 
386  return 1;
387  }
388 
389  return 1;
390 
391 }
392 
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition: valist.c:105
VPUBLIC double Vpbe_getSoluteRadius(Vpbe *thee)
Get sphere radius which bounds biomolecule.
Definition: vpbe.c:162
Definition: vhal.h:211
VPUBLIC double Vpbe_getSoluteCharge(Vpbe *thee)
Get total solute charge.
Definition: vpbe.c:186
VPUBLIC double * Vpbe_getSoluteCenter(Vpbe *thee)
Get coordinates of solute center.
Definition: vpbe.c:107
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition: vhal.h:563
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition: valist.c:115
#define Vunit_kb
Boltzmann constant.
Definition: vunit.h:96
Definition: vhal.h:213
VPUBLIC int Vmgrid_value(Vmgrid *thee, double pt[3], double *value)
Get potential value (from mesh or approximation) at a point.
Definition: vmgrid.c:107
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition: vatom.c:63
#define Vunit_eps0
Vacuum permittivity.
Definition: vunit.h:108
Potential oracle for Cartesian mesh data.
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
Definition: vatom.c:105
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
Definition: vatom.c:119
#define Vunit_ec
Charge of an electron in C.
Definition: vunit.h:92