114VPUBLIC int Vopot_pot(Vopot *thee, double pt[3], double *value) {
118 double u, T, charge, eps_w, xkappa, dist, size, val, *position;
121 VASSERT(thee != VNULL);
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);
130 /* See if we're on the mesh */
137 switch (thee->bcfl) {
149 dist += VSQR(position[i] - pt[i]);
150 dist = (1.0e-10)*VSQRT(dist);
153 val = val*(exp(-xkappa*(dist-size))/(1+xkappa*size));
167 dist += VSQR(position[i] - pt[i]);
168 dist = (1.0e-10)*VSQRT(dist);
171 val = val*(exp(-xkappa*(dist-size))/(1+xkappa*size));
178 Vnm_print(2,
"Vopot_pot: Invalid bcfl flag (%d)!\n",
183 Vnm_print(2,
"Vopot_pot: Invalid bcfl flag (%d)!\n",
188 Vnm_print(2,
"Vopot_pot: Bogus thee->bcfl flag (%d)!\n",
214VPUBLIC int Vopot_curvature(Vopot *thee, double pt[3], int cflag,
219 double u, T, charge, eps_w, xkappa, dist, size, val, *position, zkappa2;
222 VASSERT(thee != VNULL);
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);
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");
238 switch (thee->bcfl) {
245 size = (1.0e-10)*Vpbe_getSoluteRadius(thee->pbe);
246 position = Vpbe_getSoluteCenter(thee->pbe);
247 charge = Vunit_ec*Vpbe_getSoluteCharge(thee->pbe);
250 dist += VSQR(position[i] - pt[i]);
251 dist = (1.0e-10)*VSQRT(dist);
253 u = zkappa2*(exp(-xkappa*(dist-size))/(1+xkappa*size));
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);
265 dist += VSQR(position[i] - pt[i]);
266 dist = (1.0e-10)*VSQRT(dist);
268 val = zkappa2*(exp(-xkappa*(dist-size))/(1+xkappa*size));
274 Vnm_print(2, "Vopot_pot: Invlid bcfl (%d)!\n", thee->bcfl);
278 Vnm_print(2, "Vopot_pot: Invlid bcfl (%d)!\n", thee->bcfl);
282 Vnm_print(2, "Vopot_pot: Bogus thee->bcfl flag (%d)!\n",
300VPUBLIC int Vopot_gradient(Vopot *thee, double pt[3], double grad[3]) {
304 double T, charge, eps_w, xkappa, size, val, *position;
305 double dx, dy, dz, dist;
308 VASSERT(thee != VNULL);
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);
316 if (!Vmgrid_gradient(thee->mgrid, pt, grad)) {
318 switch (thee->bcfl) {
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);
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);
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);
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));
372 Vnm_print(2, "Vopot: Invalid bcfl (%d)!\n", thee->bcfl);
376 Vnm_print(2, "Vopot: Invalid bcfl (%d)!\n", thee->bcfl);
380 Vnm_print(2, "Vopot_pot: Bogus thee->bcfl flag (%d)!\n",