APBS 3.0.0
Loading...
Searching...
No Matches
vopot.c
Go to the documentation of this file.
1
57#include "vopot.h"
58
59VEMBED(rcsid="$Id$")
60
61/* ///////////////////////////////////////////////////////////////////////////
62// Routine: Vopot_ctor
63// Author: Nathan Baker
65VPUBLIC 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
80VPUBLIC 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
94VPUBLIC 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
107VPUBLIC 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))
114VPUBLIC 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
214VPUBLIC 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
300VPUBLIC 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}
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition valist.c:115
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition valist.c:105
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
Definition vatom.c:105
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition vatom.c:63
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
Definition vatom.c:119
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition vhal.h:556
@ BCFL_FOCUS
Definition vhal.h:214
@ BCFL_MDH
Definition vhal.h:211
@ BCFL_UNUSED
Definition vhal.h:213
@ BCFL_SDH
Definition vhal.h:209
@ BCFL_ZERO
Definition vhal.h:208
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 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
VPUBLIC double Vpbe_getSoluteRadius(Vpbe *thee)
Get sphere radius which bounds biomolecule.
Definition vpbe.c:162
#define Vunit_ec
Charge of an electron in C.
Definition vunit.h:92
#define Vunit_kb
Boltzmann constant.
Definition vunit.h:96
#define Vunit_eps0
Vacuum permittivity.
Definition vunit.h:108
Potential oracle for Cartesian mesh data.