APBS  1.5
vfetk.c
Go to the documentation of this file.
1 
57 #include "vfetk.h"
58 
59 /* Define the macro DONEUMANN to run with all-Neumann boundary conditions.
60  * Set this macro at your own risk! */
61 /* #define DONEUMANN 1 */
62 
63 /*
64  * @brief Calculuate the contribution to the charge-potential energy from one
65  * atom
66  * @ingroup Vfetk
67  * @author Nathan Baker
68  * @param thee current Vfetk object
69  * @param iatom current atom index
70  * @param color simplex subset (partition) under consideration
71  * @param sol current solution
72  * @returns Per-atom energy
73  */
74 VPRIVATE double Vfetk_qfEnergyAtom(
75  Vfetk *thee,
76  int iatom,
77  int color,
78  double *sol
79  );
80 
81 /*
82  * @brief Container for local variables
83  * @ingroup Vfetk
84  * @bug Not thread-safe
85  */
86 VPRIVATE Vfetk_LocalVar var;
87 
88 /*
89  * @brief MCSF-format cube mesh (all Dirichlet)
90  * @ingroup Vfetk
91  * @author Based on mesh by Mike Holst
92  */
93 VPRIVATE char *diriCubeString =
94 "mcsf_begin=1;\n\
95 \n\
96 dim=3;\n\
97 dimii=3;\n\
98 vertices=8;\n\
99 simplices=6;\n\
100 \n\
101 vert=[\n\
102 0 0 -0.5 -0.5 -0.5\n\
103 1 0 0.5 -0.5 -0.5\n\
104 2 0 -0.5 0.5 -0.5\n\
105 3 0 0.5 0.5 -0.5\n\
106 4 0 -0.5 -0.5 0.5\n\
107 5 0 0.5 -0.5 0.5\n\
108 6 0 -0.5 0.5 0.5\n\
109 7 0 0.5 0.5 0.5\n\
110 ];\n\
111 \n\
112 simp=[\n\
113 0 0 0 0 1 0 1 0 5 1 2\n\
114 1 0 0 0 1 1 0 0 5 2 4\n\
115 2 0 0 0 1 0 1 1 5 3 2\n\
116 3 0 0 0 1 0 1 3 5 7 2\n\
117 4 0 0 1 1 0 0 2 5 7 6\n\
118 5 0 0 1 1 0 0 2 5 6 4\n\
119 ];\n\
120 \n\
121 mcsf_end=1;\n\
122 \n\
123 ";
124 
125 /*
126  * @brief MCSF-format cube mesh (all Neumann)
127  * @ingroup Vfetk
128  * @author Based on mesh by Mike Holst
129  */
130 VPRIVATE char *neumCubeString =
131 "mcsf_begin=1;\n\
132 \n\
133 dim=3;\n\
134 dimii=3;\n\
135 vertices=8;\n\
136 simplices=6;\n\
137 \n\
138 vert=[\n\
139 0 0 -0.5 -0.5 -0.5\n\
140 1 0 0.5 -0.5 -0.5\n\
141 2 0 -0.5 0.5 -0.5\n\
142 3 0 0.5 0.5 -0.5\n\
143 4 0 -0.5 -0.5 0.5\n\
144 5 0 0.5 -0.5 0.5\n\
145 6 0 -0.5 0.5 0.5\n\
146 7 0 0.5 0.5 0.5\n\
147 ];\n\
148 \n\
149 simp=[\n\
150 0 0 0 0 2 0 2 0 5 1 2\n\
151 1 0 0 0 2 2 0 0 5 2 4\n\
152 2 0 0 0 2 0 2 1 5 3 2\n\
153 3 0 0 0 2 0 2 3 5 7 2\n\
154 4 0 0 2 2 0 0 2 5 7 6\n\
155 5 0 0 2 2 0 0 2 5 6 4\n\
156 ];\n\
157 \n\
158 mcsf_end=1;\n\
159 \n\
160 ";
161 
162 /*
163  * @brief Return the smoothed value of the dielectric coefficient at the
164  * current point using a fast, chart-based method
165  * @ingroup Vfetk
166  * @author Nathan Baker
167  * @returns Value of dielectric coefficient
168  * @bug Not thread-safe
169  */
170 VPRIVATE double diel();
171 
172 /*
173  * @brief Return the smoothed value of the ion accessibility at the
174  * current point using a fast, chart-based method
175  * @ingroup Vfetk
176  * @author Nathan Baker
177  * @returns Value of mobile ion coefficient
178  * @bug Not thread-safe
179  */
180 VPRIVATE double ionacc();
181 
182 /*
183  * @brief Smooths a mesh-based coefficient with a simple harmonic function
184  * @ingroup Vfetk
185  * @author Nathan Baker
186  * @param meth Method for smoothing
187  * \li 0 ==> arithmetic mean (gives bad results)
188  * \li 1 ==> geometric mean
189  * @param nverts Number of vertices
190  * @param dist distance from point to each vertex
191  * @param coeff coefficient value at each vertex
192  * @note Thread-safe
193  * @return smoothed value of coefficieent at point of interest */
194 VPRIVATE double smooth(
195  int nverts,
196  double dist[VAPBS_NVS],
197  double coeff[VAPBS_NVS],
198  int meth
199  );
200 
201 
202 /*
203  * @brief Return the analytical multi-sphere Debye-Huckel approximation (in
204  * kT/e) at the specified point
205  * @ingroup Vfetk
206  * @author Nathan Baker
207  * @param pbe Vpbe object
208  * @param d Dimension of x
209  * @param x Coordinates of point of interest (in Å)
210  * @note Thread-safe
211  * @returns Multi-sphere Debye-Huckel potential in kT/e
212  */
213 VPRIVATE double debye_U(
214  Vpbe *pbe,
215  int d,
216  double x[]
217  );
218 
219 /*
220  * @brief Return the difference between the analytical multi-sphere
221  * Debye-Huckel approximation and Coulomb's law (in kT/e) at the specified
222  * point
223  * @ingroup Vfetk
224  * @author Nathan Baker
225  * @param pbe Vpbe object
226  * @param d Dimension of x
227  * @param x Coordinates of point of interest (in Å)
228  * @note Thread-safe
229  * @returns Multi-sphere Debye-Huckel potential in kT/e */
230 VPRIVATE double debye_Udiff(
231  Vpbe *pbe,
232  int d,
233  double x[]
234  );
235 
236 /*
237  * @brief Calculate the Coulomb's
238  * Debye-Huckel approximation and Coulomb's law (in kT/e) at the specified
239  * point
240  * @ingroup Vfetk
241  * @author Nathan Baker
242  * @param pbe Vpbe object
243  * @param d Dimension of x
244  * @param x Coordinates of point of interest (in Å)
245  * @param eps Dielectric constant
246  * @param U Set to potential (in kT/e)
247  * @param dU Set to potential gradient (in kT/e/Å)
248  * @param d2U Set to Laplacian of potential (in \f$kT e^{-1} \AA^{-2}\f$)
249  * @returns Multi-sphere Debye-Huckel potential in kT/e */
250 VPRIVATE void coulomb(
251  Vpbe *pbe,
252  int d,
253  double x[],
254  double eps,
255  double *U,
256  double dU[],
257  double *d2U
258  );
259 
260 /*
261  * @brief 2D linear master simplex information generator
262  * @ingroup Vfetk
263  * @author Mike Holst
264  * @param dimIS dunno
265  * @param ndof dunno
266  * @param dof dunno
267  * @param c dunno
268  * @param cx dunno
269  * @note Trust in Mike */
270 VPRIVATE void init_2DP1(
271  int dimIS[],
272  int *ndof,
273  int dof[],
274  double c[][VMAXP],
275  double cx[][VMAXP],
276  double cy[][VMAXP],
277  double cz[][VMAXP]
278  );
279 
280 /*
281  * @brief 3D linear master simplex information generator
282  * @ingroup Vfetk
283  * @author Mike Holst
284  * @param dimIS dunno
285  * @param ndof dunno
286  * @param dof dunno
287  * @param c dunno
288  * @param cx dunno
289  * @param cy dunno
290  * @param cz dunno
291  * @note Trust in Mike */
292 VPRIVATE void init_3DP1(
293  int dimIS[],
294  int *ndof,
295  int dof[],
296  double c[][VMAXP],
297  double cx[][VMAXP],
298  double cy[][VMAXP],
299  double cz[][VMAXP]
300  );
301 
302 /*
303  * @brief Setup coefficients of polynomials from integer table data
304  * @ingroup Vfetk
305  * @author Mike Holst
306  * @param numP dunno
307  * @param c dunno
308  * @param cx dunno
309  * @param cy dunno
310  * @param cz dunno
311  * @param ic dunno
312  * @param icx dunno
313  * @param icy dunno
314  * @param icz dunno
315  * @note Trust in Mike */
316 VPRIVATE void setCoef(
317  int numP,
318  double c[][VMAXP],
319  double cx[][VMAXP],
320  double cy[][VMAXP],
321  double cz[][VMAXP],
322  int ic[][VMAXP],
323  int icx[][VMAXP],
324  int icy[][VMAXP],
325  int icz[][VMAXP]
326  );
327 
328 /*
329  * @brief Evaluate a collection of at most cubic polynomials at a
330  * specified point in at most R^3.
331  * @ingroup Vfetk
332  * @author Mike Holst
333  * @param numP the number of polynomials to evaluate
334  * @param p the results of the evaluation
335  * @param c the coefficients of each polynomial
336  * @param xv the point (x,y,z) to evaluate the polynomials.
337  * @note Mike says:
338  * <pre>
339  * Note that "VMAXP" must be >= 19 for cubic polynomials.
340  * The polynomials are build from the coefficients c[][] as
341  * follows. To build polynomial "k", fix k and set:
342  *
343  * c0=c[k][0], c1=c[k][1], .... , cp=c[k][p]
344  *
345  * Then evaluate as:
346  *
347  * p3(x,y,z) = c0 + c1*x + c2*y + c3*z
348  * + c4*x*x + c5*y*y + c6*z*z + c7*x*y + c8*x*z + c9*y*z
349  * + c10*x*x*x + c11*y*y*y + c12*z*z*z
350  * + c13*x*x*y + c14*x*x*z + c15*x*y*y
351  * + c16*y*y*z + c17*x*z*z + c18*y*z*z
352  * </pre>
353  */
354 VPRIVATE void polyEval(
355  int numP,
356  double p[],
357  double c[][VMAXP],
358  double xv[]
359  );
360 
361 /*
362  * @brief I have no clue what this variable does, but we need it to initialize
363  * the simplices
364  * @ingroup Vfetk
365  * @author Mike Holst */
366 VPRIVATE int dim_2DP1 = 3;
367 
368 /*
369  * @brief I have no clue what these variable do, but we need it to initialize
370  * the simplices
371  * @ingroup Vfetk
372  * @author Mike Holst
373  * @note Mike says:
374  * <pre>
375  * 2D-P1 Basis:
376  *
377  * p1(x,y) = c0 + c1*x + c2*y
378  *
379  * Lagrange Point Lagrange Basis Function Definition
380  * -------------- ----------------------------------
381  * (0, 0) p[0](x,y) = 1 - x - y
382  * (1, 0) p[1](x,y) = x
383  * (0, 1) p[2](x,y) = y
384  * </pre>
385  */
386 VPRIVATE int lgr_2DP1[3][VMAXP] = {
387 /*c0 c1 c2 c3
388 * ---------------------------------------------------------- */
389 /* 1 x y z
390 * ---------------------------------------------------------- */
391 { 2, -2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
392 { 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
393 { 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
394 };
395 VPRIVATE int lgr_2DP1x[3][VMAXP] = {
396 /*c0 ---------------------------------------------------------------------- */
397 /* 1 ---------------------------------------------------------------------- */
398 { -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
399 { 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
400 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
401 };
402 VPRIVATE int lgr_2DP1y[3][VMAXP] = {
403 /*c0 ---------------------------------------------------------------------- */
404 /* 1 ---------------------------------------------------------------------- */
405 { -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
406 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
407 { 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
408 };
409 VPRIVATE int lgr_2DP1z[3][VMAXP] = {
410 /*c0 ---------------------------------------------------------------------- */
411 /* 1 ---------------------------------------------------------------------- */
412 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
413 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
414 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
415 };
416 
417 
418 /*
419  * @brief I have no clue what these variable do, but we need it to initialize
420  * the simplices
421  * @ingroup Vfetk
422  * @author Mike Holst
423  * @note Mike says:
424  * <pre>
425  * 3D-P1 Basis:
426  *
427  * p1(x,y,z) = c0 + c1*x + c2*y + c3*z
428  *
429  * Lagrange Point Lagrange Basis Function Definition
430  * -------------- ----------------------------------
431  * (0, 0, 0) p[0](x,y,z) = 1 - x - y - z
432  * (1, 0, 0) p[1](x,y,z) = x
433  * (0, 1, 0) p[2](x,y,z) = y
434  * (0, 0, 1) p[3](x,y,z) = z
435  * </pre>
436  */
437 VPRIVATE int dim_3DP1 = VAPBS_NVS;
438 VPRIVATE int lgr_3DP1[VAPBS_NVS][VMAXP] = {
439 /*c0 c1 c2 c3 ---------------------------------------------------------- */
440 /* 1 x y z ---------------------------------------------------------- */
441 { 2, -2, -2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
442 { 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
443 { 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
444 { 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
445 };
446 VPRIVATE int lgr_3DP1x[VAPBS_NVS][VMAXP] = {
447 /*c0 ---------------------------------------------------------------------- */
448 /* 1 ---------------------------------------------------------------------- */
449 { -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
450 { 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
451 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
452 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
453 };
454 VPRIVATE int lgr_3DP1y[VAPBS_NVS][VMAXP] = {
455 /*c0 ---------------------------------------------------------------------- */
456 /* 1 ---------------------------------------------------------------------- */
457 { -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
458 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
459 { 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
460 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
461 };
462 VPRIVATE int lgr_3DP1z[VAPBS_NVS][VMAXP] = {
463 /*c0 ---------------------------------------------------------------------- */
464 /* 1 ---------------------------------------------------------------------- */
465 { -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
466 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
467 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
468 { 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
469 };
470 
471 /*
472  * @brief Another Holst variable
473  * @ingroup Vfetk
474  * @author Mike Holst
475  * @note Mike says: 1 = linear, 2 = quadratic */
476 VPRIVATE const int P_DEG=1;
477 
478 /*
479  * @brief Another Holst variable
480  * @ingroup Vfetk
481  * @author Mike Holst */
482 VPRIVATE int numP;
483 VPRIVATE double c[VMAXP][VMAXP];
484 VPRIVATE double cx[VMAXP][VMAXP];
485 VPRIVATE double cy[VMAXP][VMAXP];
486 VPRIVATE double cz[VMAXP][VMAXP];
487 
488 #if !defined(VINLINE_VFETK)
489 
490 VPUBLIC Gem* Vfetk_getGem(Vfetk *thee) {
491 
492  VASSERT(thee != VNULL);
493  return thee->gm;
494 
495 }
496 
497 VPUBLIC AM* Vfetk_getAM(Vfetk *thee) {
498 
499  VASSERT(thee != VNULL);
500  return thee->am;
501 }
502 
503 VPUBLIC Vpbe* Vfetk_getVpbe(Vfetk *thee) {
504 
505  VASSERT(thee != VNULL);
506  return thee->pbe;
507 
508 }
509 
510 VPUBLIC Vcsm* Vfetk_getVcsm(Vfetk *thee) {
511 
512  VASSERT(thee != VNULL);
513  return thee->csm;
514 
515 }
516 
517 VPUBLIC int Vfetk_getAtomColor(Vfetk *thee,
518  int iatom
519  ) {
520 
521  int natoms;
522 
523  VASSERT(thee != VNULL);
524 
525  natoms = Valist_getNumberAtoms(Vpbe_getValist(thee->pbe));
526  VASSERT(iatom < natoms);
527 
528  return Vatom_getPartID(Valist_getAtom(Vpbe_getValist(thee->pbe), iatom));
529 }
530 #endif /* if !defined(VINLINE_VFETK) */
531 
532 VPUBLIC Vfetk* Vfetk_ctor(Vpbe *pbe,
533  Vhal_PBEType type
534  ) {
535 
536  /* Set up the structure */
537  Vfetk *thee = VNULL;
538  thee = (Vfetk*)Vmem_malloc(VNULL, 1, sizeof(Vfetk) );
539  VASSERT(thee != VNULL);
540  VASSERT(Vfetk_ctor2(thee, pbe, type));
541 
542  return thee;
543 }
544 
545 VPUBLIC int Vfetk_ctor2(Vfetk *thee,
546  Vpbe *pbe,
547  Vhal_PBEType type
548  ) {
549 
550  int i;
551  double center[VAPBS_DIM];
552 
553  /* Make sure things have been properly initialized & store them */
554  VASSERT(pbe != VNULL);
555  thee->pbe = pbe;
556  VASSERT(pbe->alist != VNULL);
557  VASSERT(pbe->acc != VNULL);
558 
559  /* Store PBE type */
560  thee->type = type;
561 
562  /* Set up memory management object */
563  thee->vmem = Vmem_ctor("APBS::VFETK");
564 
565  /* Set up FEtk objects */
566  Vnm_print(0, "Vfetk_ctor2: Constructing PDE...\n");
567  thee->pde = Vfetk_PDE_ctor(thee);
568  Vnm_print(0, "Vfetk_ctor2: Constructing Gem...\n");
569  thee->gm = Gem_ctor(thee->vmem, thee->pde);
570  Vnm_print(0, "Vfetk_ctor2: Constructing Aprx...\n");
571  thee->aprx = Aprx_ctor(thee->vmem, thee->gm, thee->pde);
572  Vnm_print(0, "Vfetk_ctor2: Constructing Aprx...\n");
573  thee->am = AM_ctor(thee->vmem, thee->aprx);
574 
575  /* Reset refinement level */
576  thee->level = 0;
577 
578  /* Set default solver variables */
579  thee->lkey = VLT_MG;
580  thee->lmax = 1000000;
581  thee->ltol = 1e-5;
582  thee->lprec = VPT_MG;
583  thee->nkey = VNT_NEW;
584  thee->nmax = 1000000;
585  thee->ntol = 1e-5;
586  thee->gues = VGT_ZERO;
587  thee->pjac = -1;
588 
589  /* Store local copy of myself */
590  var.fetk = thee;
591  var.initGreen = 0;
592 
593  /* Set up the external Gem subdivision hook */
595 
596  /* Set up ion-related variables */
597  var.zkappa2 = Vpbe_getZkappa2(var.fetk->pbe);
599  if (var.ionstr > 0.0) var.zks2 = 0.5*var.zkappa2/var.ionstr;
600  else var.zks2 = 0.0;
601  Vpbe_getIons(var.fetk->pbe, &(var.nion), var.ionConc, var.ionRadii,
602  var.ionQ);
603  for (i=0; i<var.nion; i++) {
604  var.ionConc[i] = var.zks2 * var.ionConc[i] * var.ionQ[i];
605  }
606 
607  /* Set uninitialized objects to NULL */
608  thee->pbeparm = VNULL;
609  thee->feparm = VNULL;
610  thee->csm = VNULL;
611 
612  return 1;
613 }
614 
615 VPUBLIC void Vfetk_setParameters(Vfetk *thee,
616  PBEparm *pbeparm,
617  FEMparm *feparm
618  ) {
619 
620  VASSERT(thee != VNULL);
621  thee->feparm = feparm;
622  thee->pbeparm = pbeparm;
623 }
624 
625 VPUBLIC void Vfetk_dtor(Vfetk **thee) {
626  if ((*thee) != VNULL) {
627  Vfetk_dtor2(*thee);
628  //Vmem_free(VNULL, 1, sizeof(Vfetk), (void **)thee);
629  (*thee) = VNULL;
630  }
631 }
632 
633 VPUBLIC void Vfetk_dtor2(Vfetk *thee) {
634  Vcsm_dtor(&(thee->csm));
635  AM_dtor(&(thee->am));
636  Aprx_dtor(&(thee->aprx));
637  Vfetk_PDE_dtor(&(thee->pde));
638  Vmem_dtor(&(thee->vmem));
639 }
640 
641 VPUBLIC double* Vfetk_getSolution(Vfetk *thee,
642  int *length
643  ) {
644 
645  int i;
646  double *solution,
647  *theAnswer;
648  AM *am;
649 
650  VASSERT(thee != VNULL);
651 
652  /* Get the AM object */
653  am = thee->am;
654  /* Copy the solution into the w0 vector */
655  Bvec_copy(am->w0, am->u);
656  /* Add the Dirichlet conditions */
657  Bvec_axpy(am->w0, am->ud, 1.);
658  /* Get the data from the Bvec */
659  solution = Bvec_addr(am->w0);
660  /* Get the length of the data from the Bvec */
661  *length = Bvec_numRT(am->w0);
662  /* Make sure that we got scalar data (only one block) for the solution
663  * to the FETK */
664  VASSERT(1 == Bvec_numB(am->w0));
665  /* Allocate space for the returned vector and copy the solution into it */
666  theAnswer = VNULL;
667  theAnswer = (double*)Vmem_malloc(VNULL, *length, sizeof(double));
668  VASSERT(theAnswer != VNULL);
669  for (i=0; i<(*length); i++) theAnswer[i] = solution[i];
670 
671  return theAnswer;
672 }
673 
674 
693 VPUBLIC double Vfetk_energy(Vfetk *thee,
694  int color,
698  int nonlin
700  ) {
701 
702  double totEnergy = 0.0,
703  qfEnergy = 0.0,
704  dqmEnergy = 0.0;
706  VASSERT(thee != VNULL);
707 
708  if (nonlin && (Vpbe_getBulkIonicStrength(thee->pbe) > 0.)) {
709  Vnm_print(0, "Vfetk_energy: calculating full PBE energy\n");
710  Vnm_print(0, "Vfetk_energy: bulk ionic strength = %g M\n",
712  dqmEnergy = Vfetk_dqmEnergy(thee, color);
713  Vnm_print(0, "Vfetk_energy: dqmEnergy = %g kT\n", dqmEnergy);
714  qfEnergy = Vfetk_qfEnergy(thee, color);
715  Vnm_print(0, "Vfetk_energy: qfEnergy = %g kT\n", qfEnergy);
716 
717  totEnergy = qfEnergy - dqmEnergy;
718  } else {
719  Vnm_print(0, "Vfetk_energy: calculating only q-phi energy\n");
720  dqmEnergy = Vfetk_dqmEnergy(thee, color);
721  Vnm_print(0, "Vfetk_energy: dqmEnergy = %g kT (NOT USED)\n", dqmEnergy);
722  qfEnergy = Vfetk_qfEnergy(thee, color);
723  Vnm_print(0, "Vfetk_energy: qfEnergy = %g kT\n", qfEnergy);
724  totEnergy = 0.5*qfEnergy;
725  }
726 
727  return totEnergy;
728 
729 }
730 
731 
732 VPUBLIC double Vfetk_qfEnergy(Vfetk *thee,
733  int color
734  ) {
735 
736  double *sol,
737  energy = 0.0;
738  int nsol,
739  iatom,
740  natoms;
741  AM *am;
742 
743  VASSERT(thee != VNULL);
744  am = thee->am;
745 
746  /* Get the finest level solution */
747  sol= VNULL;
748  sol = Vfetk_getSolution(thee, &nsol);
749  VASSERT(sol != VNULL);
750 
751  /* Make sure the number of entries in the solution array matches the
752  * number of vertices currently in the mesh */
753  if (nsol != Gem_numVV(thee->gm)) {
754  Vnm_print(2, "Vfetk_qfEnergy: Number of unknowns in solution does not match\n");
755  Vnm_print(2, "Vfetk_qfEnergy: number of vertices in mesh!!! Bailing out!\n");
756  VASSERT(0);
757  }
758 
759  /* Now we do the sum over atoms... */
760  natoms = Valist_getNumberAtoms(thee->pbe->alist);
761  for (iatom=0; iatom<natoms; iatom++) {
762 
763  energy = energy + Vfetk_qfEnergyAtom(thee, iatom, color, sol);
764 
765  } /* end for iatom */
766 
767  /* Destroy the finest level solution */
768  Vmem_free(VNULL, nsol, sizeof(double), (void **)&sol);
769 
770  /* Return the energy */
771  return energy;
772 }
773 
774 VPRIVATE double Vfetk_qfEnergyAtom(
775  Vfetk *thee,
776  int iatom,
777  int color,
778  double *sol) {
779 
780  Vatom *atom;
781  double charge,
782  phi[VAPBS_NVS],
783  phix[VAPBS_NVS][3],
784  *position,
785  uval,
786  energy = 0.0;
787  int isimp,
788  nsimps,
789  icolor,
790  ivert,
791  usingColor;
792  SS *simp;
793 
794 
795  /* Get atom information */
796  atom = Valist_getAtom(thee->pbe->alist, iatom);
797  icolor = Vfetk_getAtomColor(thee, iatom);
798  charge = Vatom_getCharge(atom);
799  position = Vatom_getPosition(atom);
800 
801  /* Find out if we're using colors */
802  usingColor = (color >= 0);
803 
804  if (usingColor && (icolor<0)) {
805  Vnm_print(2, "Vfetk_qfEnergy: Atom colors not set!\n");
806  VASSERT(0);
807  }
808 
809  /* Check if this atom belongs to the specified partition */
810  if ((icolor==color) || (!usingColor)) {
811  /* Loop over the simps associated with this atom */
812  nsimps = Vcsm_getNumberSimplices(thee->csm, iatom);
813 
814  /* Get the first simp of the correct color; we can use just one
815  * simplex for energy evaluations, but not for force
816  * evaluations */
817  for (isimp=0; isimp<nsimps; isimp++) {
818 
819  simp = Vcsm_getSimplex(thee->csm, isimp, iatom);
820 
821  /* If we've asked for a particular partition AND if the atom
822  * is our partition, then compute the energy */
823  if ((SS_chart(simp)==color)||(color<0)) {
824  /* Get the value of each basis function evaluated at this
825  * point */
826  Gem_pointInSimplexVal(thee->gm, simp, position, phi, phix);
827  for (ivert=0; ivert<SS_dimVV(simp); ivert++) {
828  uval = sol[VV_id(SS_vertex(simp,ivert))];
829  energy += (charge*phi[ivert]*uval);
830  } /* end for ivert */
831  /* We only use one simplex of the appropriate color for
832  * energy calculations, so break here */
833  break;
834  } /* endif (color) */
835  } /* end for isimp */
836  }
837 
838  return energy;
839 }
840 
841 
842 VPUBLIC double Vfetk_dqmEnergy(Vfetk *thee,
843  int color) {
844 
845  return AM_evalJ(thee->am);
846 
847 }
848 
849 VPUBLIC void Vfetk_setAtomColors(Vfetk *thee) {
850 
851  SS *simp;
852  Vatom *atom;
853  int i,
854  natoms;
855 
856  VASSERT(thee != VNULL);
857 
858  natoms = Valist_getNumberAtoms(thee->pbe->alist);
859  for (i=0; i<natoms; i++) {
860  atom = Valist_getAtom(thee->pbe->alist, i);
861  simp = Vcsm_getSimplex(thee->csm, 0, i);
862  Vatom_setPartID(atom, SS_chart(simp));
863  }
864 
865 }
866 
867 VPUBLIC unsigned long int Vfetk_memChk(Vfetk *thee) {
868 
869  int memUse = 0;
870 
871  if (thee == VNULL) return 0;
872 
873  memUse = memUse + sizeof(Vfetk);
874  memUse = memUse + Vcsm_memChk(thee->csm);
875 
876  return memUse;
877 }
878 
885 VPUBLIC Vrc_Codes Vfetk_genCube(Vfetk *thee,
886  double center[3],
887  double length[3],
888  Vfetk_MeshLoad meshType
889  ) {
890 
891  VASSERT(thee != VNULL);
892 
893  AM *am = VNULL; /* @todo - no idea what this is */
894  Gem *gm = VNULL; /* Geometry manager */
895 
896  int skey = 0, /* Simplex format */
897  bufsize = 0, /* Buffer size */
898  i, /* Loop counter */
899  j; /* Loop counter */
900  char *key = "r", /* Read */
901  *iodev = "BUFF", /* Buffer */
902  *iofmt = "ASC", /* ASCII */
903  *iohost = "localhost", /* localhost (dummy) */
904  *iofile = "0", /*< socket 0 (dummy) */
905  buf[VMAX_BUFSIZE]; /* Socket buffer */
906  Vio *sock = VNULL; /* Socket object */
907  VV *vx = VNULL; /* @todo - no idea what this is */
908  double x;
909 
910  am = thee->am;
911  VASSERT(am != VNULL);
912  gm = thee->gm;
913  VASSERT(gm != VNULL);
914 
915  /* @note This code is based on Gem_makeCube by Mike Holst */
916  /* Write mesh string to buffer and read back */
917  switch (meshType) {
918  case VML_DIRICUBE:
919  /* Create a new copy of the DIRICUBE mesh (see globals higher in this file) */
920  bufsize = strlen(diriCubeString);
921  VASSERT( bufsize <= VMAX_BUFSIZE );
922  strncpy(buf, diriCubeString, VMAX_BUFSIZE);
923  break;
924  case VML_NEUMCUBE:
925  /* Create a new copy of the NEUMCUBE mesh (see globals higher in this file) */
926  bufsize = strlen(neumCubeString);
927  Vnm_print(2, "Vfetk_genCube: WARNING! USING EXPERIMENTAL NEUMANN BOUNDARY CONDITIONS!\n");
928  VASSERT( bufsize <= VMAX_BUFSIZE );
929  strncpy(buf, neumCubeString, VMAX_BUFSIZE);
930  break;
931  case VML_EXTERNAL:
932  Vnm_print(2, "Vfetk_genCube: Got request for external mesh!\n");
933  Vnm_print(2, "Vfetk_genCube: How did we get here?\n");
934  return VRC_FAILURE;
935  default:
936  Vnm_print(2, "Vfetk_genCube: Unknown mesh type (%d)\n", meshType);
937  return VRC_FAILURE;
938  }
939 
940  VASSERT( VNULL != (sock=Vio_socketOpen(key,iodev,iofmt,iohost,iofile)) ); /* Open socket */
941  Vio_bufTake(sock, buf, bufsize); /* Initialize internal buffer for socket */
942  AM_read(am, skey, sock); /* Take the initial mesh from the socket and load
943  into internal AM data structure with simplex
944  format */
945  Vio_connectFree(sock); /* Purge output buffers */
946  Vio_bufGive(sock); /* Get pointer to output buffer? No assignment of return value... */
947  Vio_dtor(&sock); /* Destroy output buffer */
948 
949  /* @todo - could the following be done in a single pass? - PCE */
950  /* Scale (unit) cube - for each vertex, set the new coordinates of that
951  vertex based on the vertex length */
952  for (i=0; i<Gem_numVV(gm); i++) {
953  vx = Gem_VV(gm, i);
954  for (j=0; j<3; j++) {
955  x = VV_coord(vx, j);
956  x *= length[j];
957  VV_setCoord(vx, j, x);
958  }
959  }
960 
961  /* Add new center - for each vertex, set a new center for the vertex */
962  for (i=0; i<Gem_numVV(gm); i++) {
963  vx = Gem_VV(gm, i);
964  for (j=0; j<3; j++) {
965  x = VV_coord(vx, j);
966  x += center[j];
967  VV_setCoord(vx, j, x);
968  }
969  }
970 
971  return VRC_SUCCESS;
972 }
973 
980 VPUBLIC Vrc_Codes Vfetk_loadMesh(Vfetk *thee, /* Vfetk object to load into */
981  double center[3], /* Center for mesh (if constructed) */
982  double length[3], /* Mesh lengths (if constructed) */
983  Vfetk_MeshLoad meshType, /* Type of mesh to load */
984  Vio *sock /* Socket for external mesh data (NULL otherwise) */
985  ) {
986 
987  Vrc_Codes vrc; /* Function return codes - see vhal.h for enum */
988  int skey = 0; /* Simplex format */
989 
990  /* Load mesh from socket if external mesh, otherwise generate mesh */
991  switch (meshType) {
992  case VML_EXTERNAL:
993  if (sock == VNULL) {
994  Vnm_print(2, "Vfetk_loadMesh: Got NULL socket!\n");
995  return VRC_FAILURE;
996  }
997  AM_read(thee->am, skey, sock);
998  Vio_connectFree(sock);
999  Vio_bufGive(sock);
1000  Vio_dtor(&sock);
1001  break;
1002  case VML_DIRICUBE:
1003  case VML_NEUMCUBE:
1004  /* Create new mesh and store in thee */
1005  vrc = Vfetk_genCube(thee, center, length, meshType);
1006  if (vrc == VRC_FAILURE) return VRC_FAILURE;
1007  break;
1008  default:
1009  Vnm_print(2, "Vfetk_loadMesh: unrecognized mesh type (%d)!\n",
1010  meshType);
1011  return VRC_FAILURE;
1012  };
1013 
1014  /* Setup charge-simplex map */
1015  Vnm_print(0, "Vfetk_ctor2: Constructing Vcsm...\n");
1016  thee->csm = VNULL;
1017  /* Construct a new Vcsm with the atom list and gem data */
1018  thee->csm = Vcsm_ctor(Vpbe_getValist(thee->pbe), thee->gm);
1019  VASSERT(thee->csm != VNULL);
1020  Vcsm_init(thee->csm);
1021 
1022  return VRC_SUCCESS;
1023 }
1024 
1025 
1026 VPUBLIC void Bmat_printHB(Bmat *thee,
1027  char *fname
1028  ) {
1029 
1030  Mat *Ablock;
1031  MATsym pqsym;
1032  int i, j, jj;
1033  int *IA, *JA;
1034  double *D, *L, *U;
1035  FILE *fp;
1036 
1037  char mmtitle[72];
1038  char mmkey[] = {"8charkey"};
1039  int totc = 0, ptrc = 0, indc = 0, valc = 0;
1040  char mxtyp[] = {"RUA"}; /* Real Unsymmetric Assembled */
1041  int nrow = 0, ncol = 0, numZ = 0;
1042  int numZdigits = 0, nrowdigits = 0;
1043  int nptrline = 8, nindline = 8, nvalline = 5;
1044  char ptrfmt[] = {"(8I10) "}, ptrfmtstr[] = {"%10d"};
1045  char indfmt[] = {"(8I10) "}, indfmtstr[] = {"%10d"};
1046  char valfmt[] = {"(5E16.8) "}, valfmtstr[] = {"%16.8E"};
1047 
1048  VASSERT( thee->numB == 1 ); /* HARDWIRE FOR NOW */
1049  Ablock = thee->AD[0][0];
1050 
1051  VASSERT( Mat_format( Ablock ) == DRC_FORMAT ); /* HARDWIRE FOR NOW */
1052 
1053  pqsym = Mat_sym( Ablock );
1054 
1055  if ( pqsym == IS_SYM ) {
1056  mxtyp[1] = 'S';
1057  } else if ( pqsym == ISNOT_SYM ) {
1058  mxtyp[1] = 'U';
1059  } else {
1060  VASSERT( 0 ); /* NOT VALID */
1061  }
1062 
1063  nrow = Bmat_numRT( thee ); /* Number of rows */
1064  ncol = Bmat_numCT( thee ); /* Number of cols */
1065  numZ = Bmat_numZT( thee ); /* Number of entries */
1066 
1067  nrowdigits = (int) (log( nrow )/log( 10 )) + 1;
1068  numZdigits = (int) (log( numZ )/log( 10 )) + 1;
1069 
1070  nptrline = (int) ( 80 / (numZdigits + 1) );
1071  nindline = (int) ( 80 / (nrowdigits + 1) );
1072 
1073  sprintf(ptrfmt,"(%dI%d)",nptrline,numZdigits+1);
1074  sprintf(ptrfmtstr,"%%%dd",numZdigits+1);
1075  sprintf(indfmt,"(%dI%d)",nindline,nrowdigits+1);
1076  sprintf(indfmtstr,"%%%dd",nrowdigits+1);
1077 
1078  ptrc = (int) ( ( (ncol + 1) - 1 ) / nptrline ) + 1;
1079  indc = (int) ( (numZ - 1) / nindline ) + 1;
1080  valc = (int) ( (numZ - 1) / nvalline ) + 1;
1081 
1082  totc = ptrc + indc + valc;
1083 
1084  sprintf( mmtitle, "Sparse '%s' Matrix - Harwell-Boeing Format - '%s'",
1085  thee->name, fname );
1086 
1087  /* Step 0: Open the file for writing */
1088 
1089  fp = fopen( fname, "w" );
1090  if (fp == VNULL) {
1091  Vnm_print(2,"Bmat_printHB: Ouch couldn't open file <%s>\n",fname);
1092  return;
1093  }
1094 
1095  /* Step 1: Print the header information */
1096 
1097  fprintf( fp, "%-72s%-8s\n", mmtitle, mmkey );
1098  fprintf( fp, "%14d%14d%14d%14d%14d\n", totc, ptrc, indc, valc, 0 );
1099  fprintf( fp, "%3s%11s%14d%14d%14d\n", mxtyp, " ", nrow, ncol, numZ );
1100  fprintf( fp, "%-16s%-16s%-20s%-20s\n", ptrfmt, indfmt, valfmt, "6E13.5" );
1101 
1102  IA = Ablock->IA;
1103  JA = Ablock->JA;
1104  D = Ablock->diag;
1105  L = Ablock->offL;
1106  U = Ablock->offU;
1107 
1108  if ( pqsym == IS_SYM ) {
1109 
1110  /* Step 2: Print the pointer information */
1111 
1112  for (i=0; i<(ncol+1); i++) {
1113  fprintf( fp, ptrfmtstr, Ablock->IA[i] + (i+1) );
1114  if ( ( (i+1) % nptrline ) == 0 ) {
1115  fprintf( fp, "\n" );
1116  }
1117  }
1118 
1119  if ( ( (ncol+1) % nptrline ) != 0 ) {
1120  fprintf( fp, "\n" );
1121  }
1122 
1123  /* Step 3: Print the index information */
1124 
1125  j = 0;
1126  for (i=0; i<ncol; i++) {
1127  fprintf( fp, indfmtstr, i+1); /* diagonal */
1128  if ( ( (j+1) % nindline ) == 0 ) {
1129  fprintf( fp, "\n" );
1130  }
1131  j++;
1132  for (jj=IA[i]; jj<IA[i+1]; jj++) {
1133  fprintf( fp, indfmtstr, JA[jj] + 1 ); /* lower triangle */
1134  if ( ( (j+1) % nindline ) == 0 ) {
1135  fprintf( fp, "\n" );
1136  }
1137  j++;
1138  }
1139  }
1140 
1141  if ( ( j % nindline ) != 0 ) {
1142  fprintf( fp, "\n" );
1143  }
1144 
1145  /* Step 4: Print the value information */
1146 
1147  j = 0;
1148  for (i=0; i<ncol; i++) {
1149  fprintf( fp, valfmtstr, D[i] );
1150  if ( ( (j+1) % nvalline ) == 0 ) {
1151  fprintf( fp, "\n" );
1152  }
1153  j++;
1154  for (jj=IA[i]; jj<IA[i+1]; jj++) {
1155  fprintf( fp, valfmtstr, L[jj] );
1156  if ( ( (j+1) % nvalline ) == 0 ) {
1157  fprintf( fp, "\n" );
1158  }
1159  j++;
1160  }
1161  }
1162 
1163  if ( ( j % nvalline ) != 0 ) {
1164  fprintf( fp, "\n" );
1165  }
1166 
1167  } else { /* ISNOT_SYM */
1168 
1169  VASSERT( 0 ); /* NOT CODED YET */
1170  }
1171 
1172  /* Step 5: Close the file */
1173  fclose( fp );
1174 }
1175 
1176 VPUBLIC PDE* Vfetk_PDE_ctor(Vfetk *fetk) {
1177 
1178  PDE *thee = VNULL;
1179 
1180  thee = (PDE*)Vmem_malloc(fetk->vmem, 1, sizeof(PDE));
1181  VASSERT(thee != VNULL);
1182  VASSERT(Vfetk_PDE_ctor2(thee, fetk));
1183 
1184  return thee;
1185 }
1186 
1187 VPUBLIC int Vfetk_PDE_ctor2(PDE *thee,
1188  Vfetk *fetk
1189  ) {
1190 
1191  int i;
1192 
1193  if (thee == VNULL) {
1194  Vnm_print(2, "Vfetk_PDE_ctor2: Got NULL thee!\n");
1195  return 0;
1196  }
1197 
1198  /* Store a local copy of the Vfetk class */
1199  var.fetk = fetk;
1200 
1201  /* PDE-specific parameters and function pointers */
1202  thee->initAssemble = Vfetk_PDE_initAssemble;
1203  thee->initElement = Vfetk_PDE_initElement;
1204  thee->initFace = Vfetk_PDE_initFace;
1205  thee->initPoint = Vfetk_PDE_initPoint;
1206  thee->Fu = Vfetk_PDE_Fu;
1207  thee->Fu_v = Vfetk_PDE_Fu_v;
1208  thee->DFu_wv = Vfetk_PDE_DFu_wv;
1209  thee->delta = Vfetk_PDE_delta;
1210  thee->u_D = Vfetk_PDE_u_D;
1211  thee->u_T = Vfetk_PDE_u_T;
1212  thee->Ju = Vfetk_PDE_Ju;
1213  thee->vec = 1; /* FIX! */
1214  thee->sym[0][0] = 1;
1215  thee->est[0] = 1.0;
1216  for (i=0; i<VMAX_BDTYPE; i++) thee->bmap[0][i] = i;
1217 
1218  /* Manifold-specific function pointers */
1219  thee->bisectEdge = Vfetk_PDE_bisectEdge;
1220  thee->mapBoundary = Vfetk_PDE_mapBoundary;
1221  thee->markSimplex = Vfetk_PDE_markSimplex;
1222  thee->oneChart = Vfetk_PDE_oneChart;
1223 
1224  /* Element-specific function pointers */
1225  thee->simplexBasisInit = Vfetk_PDE_simplexBasisInit;
1226  thee->simplexBasisForm = Vfetk_PDE_simplexBasisForm;
1227 
1228  return 1;
1229 }
1230 
1231 VPUBLIC void Vfetk_PDE_dtor(PDE **thee) {
1232 
1233  if ((*thee) != VNULL) {
1234  Vfetk_PDE_dtor2(*thee);
1235  /* TODO: The following line is commented out because at the moment,
1236  there is a seg fault when deallocating at the end of a run. Since
1237  this routine is called only once at the very end, we'll leave it
1238  commented out. However, this could be a memory leak.
1239  */
1240  /* Vmem_free(var.fetk->vmem, 1, sizeof(PDE), (void **)thee); */
1241  (*thee) = VNULL;
1242  }
1243 
1244 }
1245 
1246 VPUBLIC void Vfetk_PDE_dtor2(PDE *thee) {
1247  var.fetk = VNULL;
1248 }
1249 
1250 VPRIVATE double smooth(int nverts, double dist[VAPBS_NVS], double coeff[VAPBS_NVS], int meth) {
1251 
1252  int i;
1253  double weight;
1254  double num = 0.0;
1255  double den = 0.0;
1256 
1257  for (i=0; i<nverts; i++) {
1258  if (dist[i] < VSMALL) return coeff[i];
1259  weight = 1.0/dist[i];
1260  if (meth == 0) {
1261  num += (weight * coeff[i]);
1262  den += weight;
1263  } else if (meth == 1) {
1264  /* Small coefficients reset the average to 0; we need to break out
1265  * of the loop */
1266  if (coeff[i] < VSMALL) {
1267  num = 0.0;
1268  break;
1269  } else {
1270  num += weight; den += (weight/coeff[i]);
1271  }
1272  } else VASSERT(0);
1273  }
1274 
1275  return (num/den);
1276 
1277 }
1278 
1279 VPRIVATE double diel() {
1280 
1281  int i, j;
1282  double eps, epsp, epsw, dist[5], coeff[5], srad, swin, *vx;
1283  Vsurf_Meth srfm;
1284  Vacc *acc;
1285  PBEparm *pbeparm;
1286 
1287  epsp = Vpbe_getSoluteDiel(var.fetk->pbe);
1288  epsw = Vpbe_getSolventDiel(var.fetk->pbe);
1289  VASSERT(var.fetk->pbeparm != VNULL);
1290  pbeparm = var.fetk->pbeparm;
1291  srfm = pbeparm->srfm;
1292  srad = pbeparm->srad;
1293  swin = pbeparm->swin;
1294  acc = var.fetk->pbe->acc;
1295 
1296  eps = 0;
1297 
1298  if (VABS(epsp - epsw) < VSMALL) return epsp;
1299  switch (srfm) {
1300  case VSM_MOL:
1301  eps = ((epsw-epsp)*Vacc_molAcc(acc, var.xq, srad) + epsp);
1302  break;
1303  case VSM_MOLSMOOTH:
1304  for (i=0; i<var.nverts; i++) {
1305  dist[i] = 0;
1306  vx = var.vx[i];
1307  for (j=0; j<3; j++) {
1308  dist[i] += VSQR(var.xq[j] - vx[j]);
1309  }
1310  dist[i] = VSQRT(dist[i]);
1311  coeff[i] = (epsw-epsp)*Vacc_molAcc(acc, var.xq, srad) + epsp;
1312  }
1313  eps = smooth(var.nverts, dist, coeff, 1);
1314  break;
1315  case VSM_SPLINE:
1316  eps = ((epsw-epsp)*Vacc_splineAcc(acc, var.xq, swin, 0.0) + epsp);
1317  break;
1318  default:
1319  Vnm_print(2, "Undefined surface method (%d)!\n", srfm);
1320  VASSERT(0);
1321  }
1322 
1323  return eps;
1324 }
1325 
1326 VPRIVATE double ionacc() {
1327 
1328  int i, j;
1329  double dist[5], coeff[5], irad, swin, *vx, accval;
1330  Vsurf_Meth srfm;
1331  Vacc *acc = VNULL;
1332  PBEparm *pbeparm = VNULL;
1333 
1334  VASSERT(var.fetk->pbeparm != VNULL);
1335  pbeparm = var.fetk->pbeparm;
1336  srfm = pbeparm->srfm;
1337  irad = Vpbe_getMaxIonRadius(var.fetk->pbe);
1338  swin = pbeparm->swin;
1339  acc = var.fetk->pbe->acc;
1340 
1341  if (var.zks2 < VSMALL) return 0.0;
1342  switch (srfm) {
1343  case VSM_MOL:
1344  accval = Vacc_ivdwAcc(acc, var.xq, irad);
1345  break;
1346  case VSM_MOLSMOOTH:
1347  for (i=0; i<var.nverts; i++) {
1348  dist[i] = 0;
1349  vx = var.vx[i];
1350  for (j=0; j<3; j++) {
1351  dist[i] += VSQR(var.xq[j] - vx[j]);
1352  }
1353  dist[i] = VSQRT(dist[i]);
1354  coeff[i] = Vacc_ivdwAcc(acc, var.xq, irad);
1355  }
1356  accval = smooth(var.nverts, dist, coeff, 1);
1357  break;
1358  case VSM_SPLINE:
1359  accval = Vacc_splineAcc(acc, var.xq, swin, irad);
1360  break;
1361  default:
1362  Vnm_print(2, "Undefined surface method (%d)!\n", srfm);
1363  VASSERT(0);
1364  }
1365 
1366  return accval;
1367 }
1368 
1369 VPRIVATE double debye_U(Vpbe *pbe, int d, double x[]) {
1370 
1371  double size, *position, charge, xkappa, eps_w, dist, T, pot, val;
1372  int iatom, i;
1373  Valist *alist;
1374  Vatom *atom;
1375 
1376  eps_w = Vpbe_getSolventDiel(pbe);
1377  xkappa = (1.0e10)*Vpbe_getXkappa(pbe);
1378  T = Vpbe_getTemperature(pbe);
1379  alist = Vpbe_getValist(pbe);
1380  val = 0;
1381  pot = 0;
1382 
1383  for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++) {
1384  atom = Valist_getAtom(alist, iatom);
1385  position = Vatom_getPosition(atom);
1386  charge = Vunit_ec*Vatom_getCharge(atom);
1387  size = (1e-10)*Vatom_getRadius(atom);
1388  dist = 0;
1389  for (i=0; i<d; i++) {
1390  dist += VSQR(position[i] - x[i]);
1391  }
1392  dist = (1.0e-10)*VSQRT(dist);
1393  val = (charge)/(4*VPI*Vunit_eps0*eps_w*dist);
1394  if (xkappa != 0.0) {
1395  val = val*(exp(-xkappa*(dist-size))/(1+xkappa*size));
1396  }
1397  pot = pot + val;
1398  }
1399  pot = pot*Vunit_ec/(Vunit_kb*T);
1400 
1401  return pot;
1402 }
1403 
1404 VPRIVATE double debye_Udiff(Vpbe *pbe, int d, double x[]) {
1405 
1406  double size, *position, charge, eps_p, dist, T, pot, val;
1407  double Ufull;
1408  int iatom, i;
1409  Valist *alist;
1410  Vatom *atom;
1411 
1412  Ufull = debye_U(pbe, d, x);
1413 
1414  eps_p = Vpbe_getSoluteDiel(pbe);
1415  T = Vpbe_getTemperature(pbe);
1416  alist = Vpbe_getValist(pbe);
1417  val = 0;
1418  pot = 0;
1419 
1420  for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++) {
1421  atom = Valist_getAtom(alist, iatom);
1422  position = Vatom_getPosition(atom);
1423  charge = Vunit_ec*Vatom_getCharge(atom);
1424  size = (1e-10)*Vatom_getRadius(atom);
1425  dist = 0;
1426  for (i=0; i<d; i++) {
1427  dist += VSQR(position[i] - x[i]);
1428  }
1429  dist = (1.0e-10)*VSQRT(dist);
1430  val = (charge)/(4*VPI*Vunit_eps0*eps_p*dist);
1431  pot = pot + val;
1432  }
1433  pot = pot*Vunit_ec/(Vunit_kb*T);
1434 
1435  pot = Ufull - pot;
1436 
1437  return pot;
1438 }
1439 
1440 VPRIVATE void coulomb(Vpbe *pbe, int d, double pt[], double eps, double *U,
1441  double dU[], double *d2U) {
1442 
1443  int iatom, i;
1444  double T, pot, fx, fy, fz, x, y, z, scale;
1445  double *position, charge, dist, dist2, val, vec[3], dUold[3], Uold;
1446  Valist *alist;
1447  Vatom *atom;
1448 
1449  /* Initialize variables */
1450  T = Vpbe_getTemperature(pbe);
1451  alist = Vpbe_getValist(pbe);
1452  pot = 0; fx = 0; fy = 0; fz = 0;
1453  x = pt[0]; y = pt[1]; z = pt[2];
1454 
1455  /* Calculate */
1456  if (!Vgreen_coulombD(var.green, 1, &x, &y, &z, &pot, &fx, &fy, &fz)) {
1457  Vnm_print(2, "Error calculating Green's function!\n");
1458  VASSERT(0);
1459  }
1460 
1461 
1462  /* Scale the results */
1463  scale = Vunit_ec/(eps*Vunit_kb*T);
1464  *U = pot*scale;
1465  *d2U = 0.0;
1466  dU[0] = -fx*scale;
1467  dU[1] = -fy*scale;
1468  dU[2] = -fz*scale;
1469 
1470 #if 0
1471  /* Compare with old results */
1472  val = 0.0;
1473  Uold = 0.0; dUold[0] = 0.0; dUold[1] = 0.0; dUold[2] = 0.0;
1474  for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++) {
1475  atom = Valist_getAtom(alist, iatom);
1476  position = Vatom_getPosition(atom);
1477  charge = Vatom_getCharge(atom);
1478  dist2 = 0;
1479  for (i=0; i<d; i++) {
1480  vec[i] = (position[i] - pt[i]);
1481  dist2 += VSQR(vec[i]);
1482  }
1483  dist = VSQRT(dist2);
1484 
1485  /* POTENTIAL */
1486  Uold = Uold + charge/dist;
1487 
1488  /* GRADIENT */
1489  for (i=0; i<d; i++) dUold[i] = dUold[i] + vec[i]*charge/(dist2*dist);
1490 
1491  }
1492  Uold = Uold*VSQR(Vunit_ec)*(1.0e10)/(4*VPI*Vunit_eps0*eps*Vunit_kb*T);
1493  for (i=0; i<d; i++) {
1494  dUold[i] = dUold[i]*VSQR(Vunit_ec)*(1.0e10)/(4*VPI*Vunit_eps0*eps*Vunit_kb*T);
1495  }
1496 
1497  printf("Unew - Uold = %g - %g = %g\n", *U, Uold, (*U - Uold));
1498  printf("||dUnew - dUold||^2 = %g\n", (VSQR(dU[0] - dUold[0])
1499  + VSQR(dU[1] - dUold[1]) + VSQR(dU[2] - dUold[2])));
1500  printf("dUnew[0] = %g, dUold[0] = %g\n", dU[0], dUold[0]);
1501  printf("dUnew[1] = %g, dUold[1] = %g\n", dU[1], dUold[1]);
1502  printf("dUnew[2] = %g, dUold[2] = %g\n", dU[2], dUold[2]);
1503 
1504 #endif
1505 
1506 }
1507 
1508 VPUBLIC void Vfetk_PDE_initAssemble(PDE *thee, int ip[], double rp[]) {
1509 
1510 #if 1
1511  /* Re-initialize the Green's function oracle in case the atom list has
1512  * changed */
1513  if (var.initGreen) {
1514  Vgreen_dtor(&(var.green));
1515  var.initGreen = 0;
1516  }
1517  var.green = Vgreen_ctor(var.fetk->pbe->alist);
1518  var.initGreen = 1;
1519 #else
1520  if (!var.initGreen) {
1521  var.green = Vgreen_ctor(var.fetk->pbe->alist);
1522  var.initGreen = 1;
1523  }
1524 #endif
1525 
1526 }
1527 
1528 VPUBLIC void Vfetk_PDE_initElement(PDE *thee, int elementType, int chart,
1529  double tvx[][3], void *data) {
1530 
1531  int i, j;
1532  double epsp, epsw;
1533 
1534  /* We assume that the simplex has been passed in as the void *data * *
1535  * argument. Store it */
1536  VASSERT(data != NULL);
1537  var.simp = (SS *)data;
1538 
1539  /* save the element type */
1540  var.sType = elementType;
1541 
1542  /* Grab the vertices from this simplex */
1543  var.nverts = thee->dim+1;
1544  for (i=0; i<thee->dim+1; i++) var.verts[i] = SS_vertex(var.simp, i);
1545 
1546  /* Vertex locations of this simplex */
1547  for (i=0; i<thee->dim+1; i++) {
1548  for (j=0; j<thee->dim; j++) {
1549  var.vx[i][j] = tvx[i][j];
1550  }
1551  }
1552 
1553  /* Set the dielectric constant for this element for use in the jump term *
1554  * of the residual-based error estimator. The value is set to the average
1555  * * value of the vertices */
1556  var.jumpDiel = 0; /* NOT IMPLEMENTED YET! */
1557 }
1558 
1559 VPUBLIC void Vfetk_PDE_initFace(PDE *thee, int faceType, int chart,
1560  double tnvec[]) {
1561 
1562  int i;
1563 
1564  /* unit normal vector of this face */
1565  for (i=0; i<thee->dim; i++) var.nvec[i] = tnvec[i];
1566 
1567  /* save the face type */
1568  var.fType = faceType;
1569 }
1570 
1571 VPUBLIC void Vfetk_PDE_initPoint(PDE *thee, int pointType, int chart,
1572  double txq[], double tU[], double tdU[][3]) {
1573 
1574  int i, j, ichop;
1575  double u2, coef2, eps_p;
1576  Vhal_PBEType pdetype;
1577  Vpbe *pbe = VNULL;
1578 
1579  eps_p = Vpbe_getSoluteDiel(var.fetk->pbe);
1580  pdetype = var.fetk->type;
1581  pbe = var.fetk->pbe;
1582 
1583  /* the point, the solution value and gradient, and the Coulomb value and *
1584  * gradient at the point */
1585  if ((pdetype == PBE_LRPBE) || (pdetype == PBE_NRPBE)) {
1586  coulomb(pbe, thee->dim, txq, eps_p, &(var.W), var.dW, &(var.d2W));
1587  }
1588  for (i=0; i<thee->vec; i++) {
1589  var.U[i] = tU[i];
1590  for (j=0; j<thee->dim; j++) {
1591  var.xq[j] = txq[j];
1592  var.dU[i][j] = tdU[i][j];
1593  }
1594  }
1595 
1596  /* interior form case */
1597  if (pointType == 0) {
1598 
1599  /* Get the dielectric values */
1600  var.diel = diel();
1601  var.ionacc = ionacc();
1602  var.A = var.diel;
1603  var.F = (var.diel - eps_p);
1604 
1605  switch (pdetype) {
1606 
1607  case PBE_LPBE:
1608  var.DB = var.ionacc*var.zkappa2*var.ionstr;
1609  var.B = var.DB*var.U[0];
1610  break;
1611 
1612  case PBE_NPBE:
1613 
1614  var.B = 0;
1615  var.DB = 0;
1616  if ((var.ionacc > VSMALL) && (var.zks2 > VSMALL)) {
1617  for (i=0; i<var.nion; i++) {
1618  u2 = -1.0 * var.U[0] * var.ionQ[i];
1619 
1620  /* NONLINEAR TERM */
1621  coef2 = -1.0 * var.ionacc * var.zks2 * var.ionConc[i];
1622  var.B += (coef2 * Vcap_exp(u2, &ichop));
1623  /* LINEARIZED TERM */
1624  coef2 = -1.0 * var.ionQ[i] * coef2;
1625  var.DB += (coef2 * Vcap_exp(u2, &ichop));
1626  }
1627  }
1628  break;
1629 
1630  case PBE_LRPBE:
1631  var.DB = var.ionacc*var.zkappa2*var.ionstr;
1632  var.B = var.DB*(var.U[0]+var.W);
1633  break;
1634 
1635  case PBE_NRPBE:
1636 
1637  var.B = 0;
1638  var.DB = 0;
1639  if ((var.ionacc > VSMALL) && (var.zks2 > VSMALL)) {
1640  for (i=0; i<var.nion; i++) {
1641  u2 = -1.0 * (var.U[0] + var.W) * var.ionQ[i];
1642 
1643  /* NONLINEAR TERM */
1644  coef2 = -1.0 * var.ionacc * var.zks2 * var.ionConc[i];
1645  var.B += (coef2 * Vcap_exp(u2, &ichop));
1646 
1647  /* LINEARIZED TERM */
1648  coef2 = -1.0 * var.ionQ[i] * coef2;
1649  var.DB += (coef2 * Vcap_exp(u2, &ichop));
1650  }
1651  }
1652  break;
1653 
1654  case PBE_SMPBE: /* SMPBE Temp */
1655 
1656  var.B = 0;
1657  var.DB = 0;
1658  if ((var.ionacc > VSMALL) && (var.zks2 > VSMALL)) {
1659  for (i=0; i<var.nion; i++) {
1660  u2 = -1.0 * var.U[0] * var.ionQ[i];
1661 
1662  /* NONLINEAR TERM */
1663  coef2 = -1.0 * var.ionacc * var.zks2 * var.ionConc[i];
1664  var.B += (coef2 * Vcap_exp(u2, &ichop));
1665  /* LINEARIZED TERM */
1666  coef2 = -1.0 * var.ionQ[i] * coef2;
1667  var.DB += (coef2 * Vcap_exp(u2, &ichop));
1668  }
1669  }
1670  break;
1671  default:
1672  Vnm_print(2, "Vfetk_PDE_initPoint: Unknown PBE type (%d)!\n",
1673  pdetype);
1674  VASSERT(0);
1675  break;
1676  }
1677 
1678 
1679  /* boundary form case */
1680  } else {
1681 #ifdef DONEUMANN
1682  ;
1683 #else
1684  Vnm_print(2, "Vfetk: Whoa! I just got a boundary point to evaluate (%d)!\n", pointType);
1685  Vnm_print(2, "Vfetk: Did you do that on purpose?\n");
1686 #endif
1687  }
1688 
1689 #if 0 /* THIS IS VERY NOISY! */
1691 #endif
1692 
1693 }
1694 
1695 VPUBLIC void Vfetk_PDE_Fu(PDE *thee, int key, double F[]) {
1696 
1697  //Vnm_print(2, "Vfetk_PDE_Fu: Setting error to zero!\n");
1698 
1699  F[0] = 0.;
1700 
1701 }
1702 
1703 VPUBLIC double Vfetk_PDE_Fu_v(
1704  PDE *thee,
1705  int key,
1706  double V[],
1707  double dV[][VAPBS_DIM]
1708  ) {
1709 
1710  Vhal_PBEType type;
1711  int i;
1712  double value = 0.;
1713 
1714  type = var.fetk->type;
1715 
1716  /* interior form case */
1717  if (key == 0) {
1718 
1719  for (i=0; i<thee->dim; i++) value += ( var.A * var.dU[0][i] * dV[0][i] );
1720  value += var.B * V[0];
1721 
1722  if ((type == PBE_LRPBE) || (type == PBE_NRPBE)) {
1723  for (i=0; i<thee->dim; i++) {
1724  if (var.F > VSMALL) value += (var.F * var.dW[i] * dV[0][i]);
1725  }
1726  }
1727 
1728  /* boundary form case */
1729  } else {
1730 #ifdef DONEUMANN
1731  value = 0.0;
1732 #else
1733  Vnm_print(2, "Vfetk: Whoa! I was just asked to evaluate a boundary weak form for point type %d!\n", key);
1734 #endif
1735  }
1736 
1737  var.Fu_v = value;
1738  return value;
1739 }
1740 
1741 VPUBLIC double Vfetk_PDE_DFu_wv(
1742  PDE *thee,
1743  int key,
1744  double W[],
1745  double dW[][VAPBS_DIM],
1746  double V[],
1747  double dV[][3]
1748  ) {
1749 
1750  Vhal_PBEType type;
1751  int i;
1752  double value = 0.;
1753 
1754  type = var.fetk->type;
1755 
1756  /* Interior form */
1757  if (key == 0) {
1758  value += var.DB * W[0] * V[0];
1759  for (i=0; i<thee->dim; i++) value += ( var.A * dW[0][i] * dV[0][i] );
1760 
1761  /* boundary form case */
1762  } else {
1763 #ifdef DONEUMANN
1764  value = 0.0;
1765 #else
1766  Vnm_print(2, "Vfetk: Whoa! I was just asked to evaluate a boundary weak form for point type %d!\n", key);
1767 #endif
1768  }
1769 
1770  var.DFu_wv = value;
1771  return value;
1772 }
1773 
1776 #define VRINGMAX 1000
1777 
1779 #define VATOMMAX 1000000
1780 VPUBLIC void Vfetk_PDE_delta(PDE *thee, int type, int chart, double txq[],
1781  void *user, double F[]) {
1782 
1783  int iatom, jatom, natoms, atomIndex, atomList[VATOMMAX], nAtomList;
1784  int gotAtom, numSring, isimp, ivert, sid;
1785  double *position, charge, phi[VAPBS_NVS], phix[VAPBS_NVS][3], value;
1786  Vatom *atom;
1787  Vhal_PBEType pdetype;
1788  SS *sring[VRINGMAX];
1789  VV *vertex = (VV *)user;
1790 
1791  pdetype = var.fetk->type;
1792 
1793  F[0] = 0.0;
1794 
1795  if ((pdetype == PBE_LPBE) || (pdetype == PBE_NPBE) || (pdetype == PBE_SMPBE) /* SMPBE Added */) {
1796  VASSERT( vertex != VNULL);
1797  numSring = 0;
1798  sring[numSring] = VV_firstSS(vertex);
1799  while (sring[numSring] != VNULL) {
1800  numSring++;
1801  sring[numSring] = SS_link(sring[numSring-1], vertex);
1802  }
1803  VASSERT( numSring > 0 );
1804  VASSERT( numSring <= VRINGMAX );
1805 
1806  /* Move around the simplex ring and determine the charge locations */
1807  F[0] = 0.;
1808  charge = 0.;
1809  nAtomList = 0;
1810  for (isimp=0; isimp<numSring; isimp++) {
1811  sid = SS_id(sring[isimp]);
1812  natoms = Vcsm_getNumberAtoms(Vfetk_getVcsm(var.fetk), sid);
1813  for (iatom=0; iatom<natoms; iatom++) {
1814  /* Get the delta function information * */
1815  atomIndex = Vcsm_getAtomIndex(Vfetk_getVcsm(var.fetk),
1816  iatom, sid);
1817  gotAtom = 0;
1818  for (jatom=0; jatom<nAtomList; jatom++) {
1819  if (atomList[jatom] == atomIndex) {
1820  gotAtom = 1;
1821  break;
1822  }
1823  }
1824  if (!gotAtom) {
1825  VASSERT(nAtomList < VATOMMAX);
1826  atomList[nAtomList] = atomIndex;
1827  nAtomList++;
1828 
1829  atom = Vcsm_getAtom(Vfetk_getVcsm(var.fetk), iatom, sid);
1830  charge = Vatom_getCharge(atom);
1831  position = Vatom_getPosition(atom);
1832 
1833  /* Get the test function value at the delta function I
1834  * used to do a VASSERT to make sure the point was in the
1835  * simplex (i.e., make sure round-off error isn't an
1836  * issue), but round off errors became an issue */
1837  if (!Gem_pointInSimplexVal(Vfetk_getGem(var.fetk),
1838  sring[isimp], position, phi, phix)) {
1839  if (!Gem_pointInSimplex(Vfetk_getGem(var.fetk),
1840  sring[isimp], position)) {
1841  Vnm_print(2, "delta: Both Gem_pointInSimplexVal \
1842 and Gem_pointInSimplex detected misplaced point charge!\n");
1843  Vnm_print(2, "delta: I think you have problems: \
1844 phi = {");
1845  for (ivert=0; ivert<Gem_dimVV(Vfetk_getGem(var.fetk)); ivert++) Vnm_print(2, "%e ", phi[ivert]);
1846  Vnm_print(2, "}\n");
1847  }
1848  }
1849  value = 0;
1850  for (ivert=0; ivert<Gem_dimVV(Vfetk_getGem(var.fetk)); ivert++) {
1851  if (VV_id(SS_vertex(sring[isimp], ivert)) == VV_id(vertex)) value += phi[ivert];
1852  }
1853 
1854  F[0] += (value * Vpbe_getZmagic(var.fetk->pbe) * charge);
1855  } /* if !gotAtom */
1856  } /* for iatom */
1857  } /* for isimp */
1858 
1859  } else if ((pdetype == PBE_LRPBE) || (pdetype == PBE_NRPBE)) {
1860  F[0] = 0.0;
1861  } else { VASSERT(0); }
1862 
1863  var.delta = F[0];
1864 
1865 }
1866 
1867 VPUBLIC void Vfetk_PDE_u_D(PDE *thee, int type, int chart, double txq[],
1868  double F[]) {
1869 
1870  if ((var.fetk->type == PBE_LPBE) || (var.fetk->type == PBE_NPBE) || (var.fetk->type == PBE_SMPBE) /* SMPBE Added */) {
1871  F[0] = debye_U(var.fetk->pbe, thee->dim, txq);
1872  } else if ((var.fetk->type == PBE_LRPBE) || (var.fetk->type == PBE_NRPBE)) {
1873  F[0] = debye_Udiff(var.fetk->pbe, thee->dim, txq);
1874  } else VASSERT(0);
1875 
1876  var.u_D = F[0];
1877 
1878 }
1879 
1886 VPUBLIC void Vfetk_PDE_u_T(PDE *thee, int type, int chart, double txq[],
1887  double F[]) {
1888 /*VPUBLIC void Vfetk_PDE_u_T(sPDE *thee,
1889  int type,
1890  int chart,
1891  double txq[],
1892  double F[],
1893  double dF[][3]
1894  ) { */
1895 
1896  F[0] = 0.0;
1897  var.u_T = F[0];
1898 
1899 }
1900 
1901 
1902 VPUBLIC void Vfetk_PDE_bisectEdge(int dim, int dimII, int edgeType,
1903  int chart[], double vx[][3]) {
1904 
1905  int i;
1906 
1907  for (i=0; i<dimII; i++) vx[2][i] = .5 * (vx[0][i] + vx[1][i]);
1908  chart[2] = chart[0];
1909 
1910 }
1911 
1912 VPUBLIC void Vfetk_PDE_mapBoundary(int dim, int dimII, int vertexType,
1913  int chart, double vx[3]) {
1914 
1915 }
1916 
1917 VPUBLIC int Vfetk_PDE_markSimplex(int dim, int dimII, int simplexType,
1918  int faceType[VAPBS_NVS], int vertexType[VAPBS_NVS], int chart[], double vx[][3],
1919  void *simplex) {
1920 
1921  double targetRes, edgeLength, srad, swin, myAcc, refAcc;
1922  int i, natoms;
1923  Vsurf_Meth srfm;
1924  Vhal_PBEType type;
1925  FEMparm *feparm = VNULL;
1926  PBEparm *pbeparm = VNULL;
1927  Vpbe *pbe = VNULL;
1928  Vacc *acc = VNULL;
1929  Vcsm *csm = VNULL;
1930  SS *simp = VNULL;
1931 
1932  VASSERT(var.fetk->feparm != VNULL);
1933  feparm = var.fetk->feparm;
1934  VASSERT(var.fetk->pbeparm != VNULL);;
1935  pbeparm = var.fetk->pbeparm;
1936  pbe = var.fetk->pbe;
1937  csm = Vfetk_getVcsm(var.fetk);
1938  acc = pbe->acc;
1939  targetRes = feparm->targetRes;
1940  srfm = pbeparm->srfm;
1941  srad = pbeparm->srad;
1942  swin = pbeparm->swin;
1943  simp = (SS *)simplex;
1944  type = var.fetk->type;
1945 
1946  /* Check to see if this simplex is smaller than the target size */
1947  /* NAB WARNING: I am providing face=-1 here to conform to the new MC API; however, I'm not sure if this is the correct behavior. */
1948  Gem_longestEdge(var.fetk->gm, simp, -1, &edgeLength);
1949  if (edgeLength < targetRes) return 0;
1950 
1951  /* For non-regularized PBE, check charge-simplex map */
1952  if ((type == PBE_LPBE) || (type == PBE_NPBE) || (type == PBE_SMPBE) /* SMPBE Added */) {
1953  natoms = Vcsm_getNumberAtoms(csm, SS_id(simp));
1954  if (natoms > 0) {
1955  return 1;
1956  }
1957  }
1958 
1959  /* We would like to resolve the mesh between the van der Waals surface the
1960  * max distance from this surface where there could be coefficient
1961  * changes */
1962  switch(srfm) {
1963  case VSM_MOL:
1964  refAcc = Vacc_molAcc(acc, vx[0], srad);
1965  for (i=1; i<(dim+1); i++) {
1966  myAcc = Vacc_molAcc(acc, vx[i], srad);
1967  if (myAcc != refAcc) {
1968  return 1;
1969  }
1970  }
1971  break;
1972  case VSM_MOLSMOOTH:
1973  refAcc = Vacc_molAcc(acc, vx[0], srad);
1974  for (i=1; i<(dim+1); i++) {
1975  myAcc = Vacc_molAcc(acc, vx[i], srad);
1976  if (myAcc != refAcc) {
1977  return 1;
1978  }
1979  }
1980  break;
1981  case VSM_SPLINE:
1982  refAcc = Vacc_splineAcc(acc, vx[0], swin, 0.0);
1983  for (i=1; i<(dim+1); i++) {
1984  myAcc = Vacc_splineAcc(acc, vx[i], swin, 0.0);
1985  if (myAcc != refAcc) {
1986  return 1;
1987  }
1988  }
1989  break;
1990  default:
1991  VASSERT(0);
1992  break;
1993  }
1994 
1995  return 0;
1996 }
1997 
1998 VPUBLIC void Vfetk_PDE_oneChart(int dim, int dimII, int objType, int chart[],
1999  double vx[][3], int dimV) {
2000 
2001 }
2002 
2003 VPUBLIC double Vfetk_PDE_Ju(PDE *thee, int key) {
2004 
2005  int i, ichop;
2006  double dielE, qmE, coef2, u2;
2007  double value = 0.;
2008  Vhal_PBEType type;
2009 
2010  type = var.fetk->type;
2011 
2012  /* interior form case */
2013  if (key == 0) {
2014  dielE = 0;
2015  for (i=0; i<3; i++) dielE += VSQR(var.dU[0][i]);
2016  dielE = dielE*var.diel;
2017 
2018  switch (type) {
2019  case PBE_LPBE:
2020  if ((var.ionacc > VSMALL) && (var.zkappa2 > VSMALL)) {
2021  qmE = var.ionacc*var.zkappa2*VSQR(var.U[0]);
2022  } else qmE = 0;
2023  break;
2024  case PBE_NPBE:
2025  if ((var.ionacc > VSMALL) && (var.zks2 > VSMALL)) {
2026  qmE = 0.;
2027  for (i=0; i<var.nion; i++) {
2028  coef2 = var.ionacc * var.zks2 * var.ionConc[i] * var.ionQ[i];
2029  u2 = -1.0 * (var.U[0]) * var.ionQ[i];
2030  qmE += (coef2 * (Vcap_exp(u2, &ichop) - 1.0));
2031  }
2032  } else qmE = 0;
2033  break;
2034  case PBE_LRPBE:
2035  if ((var.ionacc > VSMALL) && (var.zkappa2 > VSMALL)) {
2036  qmE = var.ionacc*var.zkappa2*VSQR((var.U[0] + var.W));
2037  } else qmE = 0;
2038  break;
2039  case PBE_NRPBE:
2040  if ((var.ionacc > VSMALL) && (var.zks2 > VSMALL)) {
2041  qmE = 0.;
2042  for (i=0; i<var.nion; i++) {
2043  coef2 = var.ionacc * var.zks2 * var.ionConc[i] * var.ionQ[i];
2044  u2 = -1.0 * (var.U[0] + var.W) * var.ionQ[i];
2045  qmE += (coef2 * (Vcap_exp(u2, &ichop) - 1.0));
2046  }
2047  } else qmE = 0;
2048  break;
2049  case PBE_SMPBE: /* SMPBE Temp */
2050  if ((var.ionacc > VSMALL) && (var.zks2 > VSMALL)) {
2051  qmE = 0.;
2052  for (i=0; i<var.nion; i++) {
2053  coef2 = var.ionacc * var.zks2 * var.ionConc[i] * var.ionQ[i];
2054  u2 = -1.0 * (var.U[0]) * var.ionQ[i];
2055  qmE += (coef2 * (Vcap_exp(u2, &ichop) - 1.0));
2056  }
2057  } else qmE = 0;
2058  break;
2059  default:
2060  Vnm_print(2, "Vfetk_PDE_Ju: Invalid PBE type (%d)!\n", type);
2061  VASSERT(0);
2062  break;
2063  }
2064 
2065  value = 0.5*(dielE + qmE)/Vpbe_getZmagic(var.fetk->pbe);
2066 
2067  /* boundary form case */
2068  } else if (key == 1) {
2069  value = 0.0;
2070 
2071  /* how did we get here? */
2072  } else VASSERT(0);
2073 
2074  return value;
2075 
2076 }
2077 
2078 VPUBLIC void Vfetk_externalUpdateFunction(SS **simps, int num) {
2079 
2080  Vcsm *csm = VNULL;
2081  int rc;
2082 
2083  VASSERT(var.fetk != VNULL);
2084  csm = Vfetk_getVcsm(var.fetk);
2085  VASSERT(csm != VNULL);
2086 
2087  rc = Vcsm_update(csm, simps, num);
2088 
2089  if (!rc) {
2090  Vnm_print(2, "Error while updating charge-simplex map!\n");
2091  VASSERT(0);
2092  }
2093 }
2094 
2095 VPRIVATE void polyEval(int numP, double p[], double c[][VMAXP], double xv[]) {
2096  int i;
2097  double x, y, z;
2098 
2099  x = xv[0];
2100  y = xv[1];
2101  z = xv[2];
2102  for (i=0; i<numP; i++) {
2103  p[i] = c[i][0]
2104  + c[i][1] * x
2105  + c[i][2] * y
2106  + c[i][3] * z
2107  + c[i][4] * x*x
2108  + c[i][5] * y*y
2109  + c[i][6] * z*z
2110  + c[i][7] * x*y
2111  + c[i][8] * x*z
2112  + c[i][9] * y*z
2113  + c[i][10] * x*x*x
2114  + c[i][11] * y*y*y
2115  + c[i][12] * z*z*z
2116  + c[i][13] * x*x*y
2117  + c[i][14] * x*x*z
2118  + c[i][15] * x*y*y
2119  + c[i][16] * y*y*z
2120  + c[i][17] * x*z*z
2121  + c[i][18] * y*z*z;
2122  }
2123 }
2124 
2125 VPRIVATE void setCoef(int numP, double c[][VMAXP], double cx[][VMAXP],
2126  double cy[][VMAXP], double cz[][VMAXP], int ic[][VMAXP], int icx[][VMAXP],
2127  int icy[][VMAXP], int icz[][VMAXP]) {
2128 
2129  int i, j;
2130  for (i=0; i<numP; i++) {
2131  for (j=0; j<VMAXP; j++) {
2132  c[i][j] = 0.5 * (double)ic[i][j];
2133  cx[i][j] = 0.5 * (double)icx[i][j];
2134  cy[i][j] = 0.5 * (double)icy[i][j];
2135  cz[i][j] = 0.5 * (double)icz[i][j];
2136  }
2137  }
2138 }
2139 
2140 VPUBLIC int Vfetk_PDE_simplexBasisInit(int key, int dim, int comp, int *ndof,
2141  int dof[]) {
2142 
2143  int qorder, bump, dimIS[VAPBS_NVS];
2144 
2145  /* necessary quadrature order to return at the end */
2146  qorder = P_DEG;
2147 
2148  /* deal with bump function requests */
2149  if ((key == 0) || (key == 1)) {
2150  bump = 0;
2151  } else if ((key == 2) || (key == 3)) {
2152  bump = 1;
2153  } else { VASSERT(0); }
2154 
2155  /* for now use same element for all components, both trial and test */
2156  if (dim==2) {
2157  /* 2D simplex dimensions */
2158  dimIS[0] = 3; /* number of vertices */
2159  dimIS[1] = 3; /* number of edges */
2160  dimIS[2] = 0; /* number of faces (3D only) */
2161  dimIS[3] = 1; /* number of simplices (always=1) */
2162  if (bump==0) {
2163  if (P_DEG==1) {
2164  init_2DP1(dimIS, ndof, dof, c, cx, cy, cz);
2165  } else if (P_DEG==2) {
2166  init_2DP1(dimIS, ndof, dof, c, cx, cy, cz);
2167  } else if (P_DEG==3) {
2168  init_2DP1(dimIS, ndof, dof, c, cx, cy, cz);
2169  } else Vnm_print(2, "..bad order..");
2170  } else if (bump==1) {
2171  if (P_DEG==1) {
2172  init_2DP1(dimIS, ndof, dof, c, cx, cy, cz);
2173  } else Vnm_print(2, "..bad order..");
2174  } else Vnm_print(2, "..bad bump..");
2175  } else if (dim==3) {
2176  /* 3D simplex dimensions */
2177  dimIS[0] = 4; /* number of vertices */
2178  dimIS[1] = 6; /* number of edges */
2179  dimIS[2] = 4; /* number of faces (3D only) */
2180  dimIS[3] = 1; /* number of simplices (always=1) */
2181  if (bump==0) {
2182  if (P_DEG==1) {
2183  init_3DP1(dimIS, ndof, dof, c, cx, cy, cz);
2184  } else if (P_DEG==2) {
2185  init_3DP1(dimIS, ndof, dof, c, cx, cy, cz);
2186  } else if (P_DEG==3) {
2187  init_3DP1(dimIS, ndof, dof, c, cx, cy, cz);
2188  } else Vnm_print(2, "..bad order..");
2189  } else if (bump==1) {
2190  if (P_DEG==1) {
2191  init_3DP1(dimIS, ndof, dof, c, cx, cy, cz);
2192  } else Vnm_print(2, "..bad order..");
2193  } else Vnm_print(2, "..bad bump..");
2194  } else Vnm_print(2, "..bad dimension..");
2195 
2196  /* save number of DF */
2197  numP = *ndof;
2198 
2199  /* return the required quarature order */
2200  return qorder;
2201 }
2202 
2203 VPUBLIC void Vfetk_PDE_simplexBasisForm(int key, int dim, int comp, int pdkey,
2204  double xq[], double basis[]) {
2205 
2206  if (pdkey == 0) {
2207  polyEval(numP, basis, c, xq);
2208  } else if (pdkey == 1) {
2209  polyEval(numP, basis, cx, xq);
2210  } else if (pdkey == 2) {
2211  polyEval(numP, basis, cy, xq);
2212  } else if (pdkey == 3) {
2213  polyEval(numP, basis, cz, xq);
2214  } else { VASSERT(0); }
2215 }
2216 
2217 VPRIVATE void init_2DP1(int dimIS[], int *ndof, int dof[], double c[][VMAXP],
2218  double cx[][VMAXP], double cy[][VMAXP], double cz[][VMAXP]) {
2219 
2220  int i;
2221 
2222  /* dof number and locations */
2223  dof[0] = 1;
2224  dof[1] = 0;
2225  dof[2] = 0;
2226  dof[3] = 0;
2227  *ndof = 0;
2228  for (i=0; i<VAPBS_NVS; i++) *ndof += dimIS[i] * dof[i];
2229  VASSERT( *ndof == dim_2DP1 );
2230  VASSERT( *ndof <= VMAXP );
2231 
2232  /* coefficients of the polynomials */
2233  setCoef( *ndof, c, cx, cy, cz, lgr_2DP1, lgr_2DP1x, lgr_2DP1y, lgr_2DP1z );
2234 }
2235 
2236 VPRIVATE void init_3DP1(int dimIS[], int *ndof, int dof[], double c[][VMAXP],
2237  double cx[][VMAXP], double cy[][VMAXP], double cz[][VMAXP]) {
2238 
2239  int i;
2240 
2241  /* dof number and locations */
2242  dof[0] = 1;
2243  dof[1] = 0;
2244  dof[2] = 0;
2245  dof[3] = 0;
2246  *ndof = 0;
2247  for (i=0; i<VAPBS_NVS; i++) *ndof += dimIS[i] * dof[i];
2248  VASSERT( *ndof == dim_3DP1 );
2249  VASSERT( *ndof <= VMAXP );
2250 
2251  /* coefficients of the polynomials */
2252  setCoef( *ndof, c, cx, cy, cz, lgr_3DP1, lgr_3DP1x, lgr_3DP1y, lgr_3DP1z );
2253 }
2254 
2255 VPUBLIC void Vfetk_dumpLocalVar() {
2256 
2257  int i;
2258 
2259  Vnm_print(1, "DEBUG: nvec = (%g, %g, %g)\n", var.nvec[0], var.nvec[1],
2260  var.nvec[2]);
2261  Vnm_print(1, "DEBUG: nverts = %d\n", var.nverts);
2262  for (i=0; i<var.nverts; i++) {
2263  Vnm_print(1, "DEBUG: verts[%d] ID = %d\n", i, VV_id(var.verts[i]));
2264  Vnm_print(1, "DEBUG: vx[%d] = (%g, %g, %g)\n", i, var.vx[i][0],
2265  var.vx[i][1], var.vx[i][2]);
2266  }
2267  Vnm_print(1, "DEBUG: simp ID = %d\n", SS_id(var.simp));
2268  Vnm_print(1, "DEBUG: sType = %d\n", var.sType);
2269  Vnm_print(1, "DEBUG: fType = %d\n", var.fType);
2270  Vnm_print(1, "DEBUG: xq = (%g, %g, %g)\n", var.xq[0], var.xq[1], var.xq[2]);
2271  Vnm_print(1, "DEBUG: U[0] = %g\n", var.U[0]);
2272  Vnm_print(1, "DEBUG: dU[0] = (%g, %g, %g)\n", var.dU[0][0], var.dU[0][1],
2273  var.dU[0][2]);
2274  Vnm_print(1, "DEBUG: W = %g\n", var.W);
2275  Vnm_print(1, "DEBUG: d2W = %g\n", var.d2W);
2276  Vnm_print(1, "DEBUG: dW = (%g, %g, %g)\n", var.dW[0], var.dW[1], var.dW[2]);
2277  Vnm_print(1, "DEBUG: diel = %g\n", var.diel);
2278  Vnm_print(1, "DEBUG: ionacc = %g\n", var.ionacc);
2279  Vnm_print(1, "DEBUG: A = %g\n", var.A);
2280  Vnm_print(1, "DEBUG: F = %g\n", var.F);
2281  Vnm_print(1, "DEBUG: B = %g\n", var.B);
2282  Vnm_print(1, "DEBUG: DB = %g\n", var.DB);
2283  Vnm_print(1, "DEBUG: nion = %d\n", var.nion);
2284  for (i=0; i<var.nion; i++) {
2285  Vnm_print(1, "DEBUG: ionConc[%d] = %g\n", i, var.ionConc[i]);
2286  Vnm_print(1, "DEBUG: ionQ[%d] = %g\n", i, var.ionQ[i]);
2287  Vnm_print(1, "DEBUG: ionRadii[%d] = %g\n", i, var.ionRadii[i]);
2288  }
2289  Vnm_print(1, "DEBUG: zkappa2 = %g\n", var.zkappa2);
2290  Vnm_print(1, "DEBUG: zks2 = %g\n", var.zks2);
2291  Vnm_print(1, "DEBUG: Fu_v = %g\n", var.Fu_v);
2292  Vnm_print(1, "DEBUG: DFu_wv = %g\n", var.DFu_wv);
2293  Vnm_print(1, "DEBUG: delta = %g\n", var.delta);
2294  Vnm_print(1, "DEBUG: u_D = %g\n", var.u_D);
2295  Vnm_print(1, "DEBUG: u_T = %g\n", var.u_T);
2296 
2297 };
2298 
2299 VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type) {
2300 
2301  int i, j, ichop;
2302  double coord[3], chi, q, conc, val;
2303  VV *vert;
2304  Bvec *u, *u_d;
2305  AM *am;
2306  Gem *gm;
2307  PBEparm *pbeparm;
2308  Vacc *acc;
2309  Vpbe *pbe;
2310 
2311  gm = thee->gm;
2312  am = thee->am;
2313  pbe = thee->pbe;
2314  pbeparm = thee->pbeparm;
2315  acc = pbe->acc;
2316 
2317  /* Make sure vec has enough rows to accomodate the vertex data */
2318  if (Bvec_numRB(vec, 0) != Gem_numVV(gm)) {
2319  Vnm_print(2, "Vfetk_fillArray: insufficient space in Bvec!\n");
2320  Vnm_print(2, "Vfetk_fillArray: Have %d, need %d!\n", Bvec_numRB(vec, 0),
2321  Gem_numVV(gm));
2322  return 0;
2323  }
2324 
2325  switch (type) {
2326 
2327  case VDT_CHARGE:
2328  Vnm_print(2, "Vfetk_fillArray: can't write out charge distribution!\n");
2329  return 0;
2330  break;
2331 
2332  case VDT_POT:
2333  u = am->u;
2334  u_d = am->ud;
2335  /* Copy in solution */
2336  Bvec_copy(vec, u);
2337  /* Add dirichlet condition */
2338  Bvec_axpy(vec, u_d, 1.0);
2339  break;
2340 
2341  case VDT_SMOL:
2342  for (i=0; i<Gem_numVV(gm); i++) {
2343  vert = Gem_VV(gm, i);
2344  for (j=0; j<3; j++) coord[j] = VV_coord(vert, j);
2345  chi = Vacc_molAcc(acc, coord, pbe->solventRadius);
2346  Bvec_set(vec, i, chi);
2347  }
2348  break;
2349 
2350  case VDT_SSPL:
2351  for (i=0; i<Gem_numVV(gm); i++) {
2352  vert = Gem_VV(gm, i);
2353  for (j=0; j<3; j++) coord[j] = VV_coord(vert, j);
2354  chi = Vacc_splineAcc(acc, coord, pbeparm->swin, 0.0);
2355  Bvec_set(vec, i, chi);
2356  }
2357  break;
2358 
2359  case VDT_VDW:
2360  for (i=0; i<Gem_numVV(gm); i++) {
2361  vert = Gem_VV(gm, i);
2362  for (j=0; j<3; j++) coord[j] = VV_coord(vert, j);
2363  chi = Vacc_vdwAcc(acc, coord);
2364  Bvec_set(vec, i, chi);
2365  }
2366  break;
2367 
2368  case VDT_IVDW:
2369  for (i=0; i<Gem_numVV(gm); i++) {
2370  vert = Gem_VV(gm, i);
2371  for (j=0; j<3; j++) coord[j] = VV_coord(vert, j);
2372  chi = Vacc_ivdwAcc(acc, coord, pbe->maxIonRadius);
2373  Bvec_set(vec, i, chi);
2374  }
2375  break;
2376 
2377  case VDT_LAP:
2378  Vnm_print(2, "Vfetk_fillArray: can't write out Laplacian!\n");
2379  return 0;
2380  break;
2381 
2382  case VDT_EDENS:
2383  Vnm_print(2, "Vfetk_fillArray: can't write out energy density!\n");
2384  return 0;
2385  break;
2386 
2387  case VDT_NDENS:
2388  u = am->u;
2389  u_d = am->ud;
2390  /* Copy in solution */
2391  Bvec_copy(vec, u);
2392  /* Add dirichlet condition */
2393  Bvec_axpy(vec, u_d, 1.0);
2394  /* Load up ions */
2395  ichop = 0;
2396  for (i=0; i<Gem_numVV(gm); i++) {
2397  val = 0;
2398  for (j=0; j<pbe->numIon; j++) {
2399  q = pbe->ionQ[j];
2400  conc = pbe->ionConc[j];
2401  if (thee->type == PBE_NPBE || thee->type == PBE_SMPBE /* SMPBE Added */) {
2402  val += (conc*Vcap_exp(-q*Bvec_val(vec, i), &ichop));
2403  } else if (thee->type == PBE_LPBE) {
2404  val += (conc * ( 1 - q*Bvec_val(vec, i)));
2405  }
2406  }
2407  Bvec_set(vec, i, val);
2408  }
2409  break;
2410 
2411  case VDT_QDENS:
2412  u = am->u;
2413  u_d = am->ud;
2414  /* Copy in solution */
2415  Bvec_copy(vec, u);
2416  /* Add dirichlet condition */
2417  Bvec_axpy(vec, u_d, 1.0);
2418  /* Load up ions */
2419  ichop = 0;
2420  for (i=0; i<Gem_numVV(gm); i++) {
2421  val = 0;
2422  for (j=0; j<pbe->numIon; j++) {
2423  q = pbe->ionQ[j];
2424  conc = pbe->ionConc[j];
2425  if (thee->type == PBE_NPBE || thee->type == PBE_SMPBE /* SMPBE Added */) {
2426  val += (q*conc*Vcap_exp(-q*Bvec_val(vec, i), &ichop));
2427  } else if (thee->type == PBE_LPBE) {
2428  val += (q*conc*(1 - q*Bvec_val(vec, i)));
2429  }
2430  }
2431  Bvec_set(vec, i, val);
2432  }
2433  break;
2434 
2435  case VDT_DIELX:
2436  Vnm_print(2, "Vfetk_fillArray: can't write out x-shifted diel!\n");
2437  return 0;
2438  break;
2439 
2440  case VDT_DIELY:
2441  Vnm_print(2, "Vfetk_fillArray: can't write out y-shifted diel!\n");
2442  return 0;
2443  break;
2444 
2445  case VDT_DIELZ:
2446  Vnm_print(2, "Vfetk_fillArray: can't write out z-shifted diel!\n");
2447  return 0;
2448  break;
2449 
2450  case VDT_KAPPA:
2451  Vnm_print(2, "Vfetk_fillArray: can't write out kappa!\n");
2452  return 0;
2453  break;
2454 
2455  default:
2456  Vnm_print(2, "Vfetk_fillArray: invalid data type (%d)!\n", type);
2457  return 0;
2458  break;
2459  }
2460 
2461  return 1;
2462 }
2463 
2464 VPUBLIC int Vfetk_write(Vfetk *thee, const char *iodev, const char *iofmt,
2465  const char *thost, const char *fname, Bvec *vec, Vdata_Format format) {
2466 
2467  int i, j, ichop;
2468  Aprx *aprx;
2469  Gem *gm;
2470  Vio *sock;
2471 
2472  VASSERT(thee != VNULL);
2473  aprx = thee->aprx;
2474  gm = thee->gm;
2475 
2476  sock = Vio_ctor(iodev,iofmt,thost,fname,"w");
2477  if (sock == VNULL) {
2478  Vnm_print(2, "Vfetk_write: Problem opening virtual socket %s\n",
2479  fname);
2480  return 0;
2481  }
2482  if (Vio_connect(sock, 0) < 0) {
2483  Vnm_print(2, "Vfetk_write: Problem connecting to virtual socket %s\n",
2484  fname);
2485  return 0;
2486  }
2487 
2488  /* Make sure vec has enough rows to accomodate the vertex data */
2489  if (Bvec_numRB(vec, 0) != Gem_numVV(gm)) {
2490  Vnm_print(2, "Vfetk_fillArray: insufficient space in Bvec!\n");
2491  Vnm_print(2, "Vfetk_fillArray: Have %d, need %d!\n", Bvec_numRB(vec, 0),
2492  Gem_numVV(gm));
2493  return 0;
2494  }
2495 
2496  switch (format) {
2497 
2498  case VDF_DX:
2499  Aprx_writeSOL(aprx, sock, vec, "DX");
2500  break;
2501  case VDF_AVS:
2502  Aprx_writeSOL(aprx, sock, vec, "UCD");
2503  break;
2504  case VDF_UHBD:
2505  Vnm_print(2, "Vfetk_write: UHBD format not supported!\n");
2506  return 0;
2507  default:
2508  Vnm_print(2, "Vfetk_write: Invalid data format (%d)!\n", format);
2509  return 0;
2510  }
2511 
2512 
2513  Vio_connectFree(sock);
2514  Vio_dtor(&sock);
2515 
2516  return 1;
2517 }
2518 
double ionstr
Definition: vfetk.h:246
Vacc * acc
Definition: vpbe.h:90
Definition: vhal.h:271
double xq[VAPBS_DIM]
Definition: vfetk.h:218
VEXTERNC void Vfetk_PDE_initPoint(PDE *thee, int pointType, int chart, double txq[], double tU[], double tdU[][VAPBS_DIM])
Do once-per-point initialization.
double jumpDiel
Definition: vfetk.h:232
double targetRes
Definition: femparm.h:160
Definition: vfetk.h:122
double DFu_wv
Definition: vfetk.h:249
VPUBLIC void Vfetk_PDE_delta(PDE *thee, int type, int chart, double txq[], void *user, double F[])
Evaluate a (discretized) delta function source term at the given point.
Definition: vfetk.c:1780
double ionQ[MAXION]
Definition: vpbe.h:106
#define VRINGMAX
Maximum number of simplices in a simplex ring.
Definition: vfetk.c:1776
VEXTERNC int Vfetk_PDE_markSimplex(int dim, int dimII, int simplexType, int faceType[VAPBS_NVS], int vertexType[VAPBS_NVS], int chart[], double vx[][VAPBS_DIM], void *simplex)
User-defined error estimator – in our case, a geometry-based refinement method; forcing simplex ref...
FEMparm * feparm
Definition: vfetk.h:198
enum eVdata_Format Vdata_Format
Declaration of the Vdata_Format type as the Vdata_Format enum.
Definition: vhal.h:323
PBEparm * pbeparm
Definition: vfetk.h:197
int lmax
Definition: vfetk.h:188
VPUBLIC Valist * Vpbe_getValist(Vpbe *thee)
Get atom list.
Definition: vpbe.c:69
double u_T
Definition: vfetk.h:252
VPUBLIC void Vfetk_setAtomColors(Vfetk *thee)
Transfer color (partition ID) information frmo a partitioned mesh to the atoms.
Definition: vfetk.c:849
VPUBLIC Vcsm * Vcsm_ctor(Valist *alist, Gem *gm)
Construct Vcsm object.
Definition: vcsm.c:136
VPUBLIC int Vcsm_update(Vcsm *thee, SS **simps, int num)
Update the charge-simplex and simplex-charge maps after refinement.
Definition: vcsm.c:326
enum eVsurf_Meth Vsurf_Meth
Declaration of the Vsurf_Meth type as the Vsurf_Meth enum.
Definition: vhal.h:133
Vfetk LocalVar subclass.
Definition: vfetk.h:215
VPUBLIC Vfetk * Vfetk_ctor(Vpbe *pbe, Vhal_PBEType type)
Constructor for Vfetk object.
Definition: vfetk.c:532
VPUBLIC void Bmat_printHB(Bmat *thee, char *fname)
Writes a Bmat to disk in Harwell-Boeing sparse matrix format.
Definition: vfetk.c:1026
Gem * gm
Definition: vfetk.h:179
VPUBLIC Vcsm * Vfetk_getVcsm(Vfetk *thee)
Get a pointer to the Vcsm (charge-simplex map) object.
Definition: vfetk.c:510
Vfetk_LsolvType lkey
Definition: vfetk.h:187
double vx[4][VAPBS_DIM]
Definition: vfetk.h:217
Valist * alist
Definition: vcsm.h:91
VPUBLIC void Vfetk_PDE_simplexBasisForm(int key, int dim, int comp, int pdkey, double xq[], double basis[])
Evaluate the bases for the trial or test space, for a particular component of the system,...
Definition: vfetk.c:2203
Contains public data members for Vpbe class/module.
Definition: vpbe.h:84
Oracle for solvent- and ion-accessibility around a biomolecule.
Definition: vacc.h:108
VPUBLIC double Vpbe_getSoluteDiel(Vpbe *thee)
Get solute dielectric constant.
Definition: vpbe.c:99
double ionConc[MAXION]
Definition: vfetk.h:241
double A
Definition: vfetk.h:228
VPUBLIC double Vpbe_getZkappa2(Vpbe *thee)
Get modified squared Debye-Huckel parameter.
Definition: vpbe.c:148
VPUBLIC void Vfetk_externalUpdateFunction(SS **simps, int num)
External hook to simplex subdivision routines in Gem. Called each time a simplex is subdivided (we us...
Definition: vfetk.c:2078
VPUBLIC void Vfetk_dtor2(Vfetk *thee)
FORTRAN stub object destructor.
Definition: vfetk.c:633
#define VATOMMAX
Maximum number of atoms associated with a vertex.
Definition: vfetk.c:1779
VPUBLIC int Vfetk_PDE_simplexBasisInit(int key, int dim, int comp, int *ndof, int dof[])
Initialize the bases for the trial or the test space, for a particular component of the system,...
Definition: vfetk.c:2140
VEXTERNC void Vfetk_PDE_mapBoundary(int dim, int dimII, int vertexType, int chart, double vx[VAPBS_DIM])
Map a boundary point to some pre-defined shape.
struct sVfetk Vfetk
Declaration of the Vfetk class as the Vfetk structure.
Definition: vfetk.h:207
VPUBLIC void Vfetk_PDE_dtor(PDE **thee)
Destroys FEtk PDE object.
Definition: vfetk.c:1231
VPUBLIC Vpbe * Vfetk_getVpbe(Vfetk *thee)
Get a pointer to the Vpbe (PBE manager) object.
Definition: vfetk.c:503
double swin
Definition: pbeparm.h:154
VEXTERNC double Vacc_ivdwAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report inflated van der Waals accessibility.
VPUBLIC int Vcsm_getNumberSimplices(Vcsm *thee, int iatom)
Get number of simplices associated with an atom.
Definition: vcsm.c:99
double U[MAXV]
Definition: vfetk.h:219
int numIon
Definition: vpbe.h:103
VPUBLIC void Vfetk_PDE_initAssemble(PDE *thee, int ip[], double rp[])
Do once-per-assembly initialization.
Definition: vfetk.c:1508
int nmax
Definition: vfetk.h:191
VPUBLIC PDE * Vfetk_PDE_ctor(Vfetk *fetk)
Constructs the FEtk PDE object.
Definition: vfetk.c:1176
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition: valist.c:105
double u_D
Definition: vfetk.h:251
VPUBLIC double Vatom_getPartID(Vatom *thee)
Get partition ID.
Definition: vatom.c:70
double d2W
Definition: vfetk.h:223
VEXTERNC double Vfetk_PDE_DFu_wv(PDE *thee, int key, double W[], double dW[][VAPBS_DIM], double V[], double dV[][VAPBS_DIM])
This is the linearization of the weak form of the PBE; e.g., for use in a Newton iteration....
VPUBLIC double Vpbe_getTemperature(Vpbe *thee)
Get temperature.
Definition: vpbe.c:91
VPUBLIC double Vpbe_getBulkIonicStrength(Vpbe *thee)
Get bulk ionic strength.
Definition: vpbe.c:84
Vmem * vmem
Definition: vfetk.h:178
VPUBLIC Vgreen * Vgreen_ctor(Valist *alist)
Construct the Green's function oracle.
Definition: vgreen.c:156
VPUBLIC void Vatom_setPartID(Vatom *thee, int partID)
Set partition ID.
Definition: vatom.c:77
VPUBLIC void Vfetk_dumpLocalVar()
Debugging routine to print out local variables used by PDE object.
Definition: vfetk.c:2255
VPUBLIC AM * Vfetk_getAM(Vfetk *thee)
Get a pointer to the AM (algebra manager) object.
Definition: vfetk.c:497
VPUBLIC int Vfetk_getAtomColor(Vfetk *thee, int iatom)
Get the partition information for a particular atom.
Definition: vfetk.c:517
double nvec[VAPBS_DIM]
Definition: vfetk.h:216
VPUBLIC SS * Vcsm_getSimplex(Vcsm *thee, int isimp, int iatom)
Get particular simplex associated with an atom.
Definition: vcsm.c:109
VPUBLIC void Vfetk_PDE_u_D(PDE *thee, int type, int chart, double txq[], double F[])
Evaluate the Dirichlet boundary condition at the given point.
Definition: vfetk.c:1867
VPUBLIC int Vpbe_getIons(Vpbe *thee, int *nion, double ionConc[MAXION], double ionRadii[MAXION], double ionQ[MAXION])
Get information about the counterion species present.
Definition: vpbe.c:535
VPUBLIC Vrc_Codes Vfetk_loadMesh(Vfetk *thee, double center[3], double length[3], Vfetk_MeshLoad meshType, Vio *sock)
Loads a mesh into the Vfetk (and associated) object(s).
Definition: vfetk.c:980
Vsurf_Meth srfm
Definition: pbeparm.h:150
Vfetk_PrecType lprec
Definition: vfetk.h:194
VPUBLIC double Vcap_exp(double x, int *ichop)
Provide a capped exp() function.
Definition: vcap.c:59
VPUBLIC int Vfetk_ctor2(Vfetk *thee, Vpbe *pbe, Vhal_PBEType type)
FORTRAN stub constructor for Vfetk object.
Definition: vfetk.c:545
VPUBLIC Vatom * Vcsm_getAtom(Vcsm *thee, int iatom, int isimp)
Get particular atom associated with a simplex.
Definition: vcsm.c:77
VEXTERNC void Gem_setExternalUpdateFunction(Gem *thee, void(*externalUpdate)(SS **simps, int num))
External function for FEtk Gem class to use during mesh refinement.
VPUBLIC void Vfetk_dtor(Vfetk **thee)
Object destructor.
Definition: vfetk.c:625
VPUBLIC void Vfetk_setParameters(Vfetk *thee, PBEparm *pbeparm, FEMparm *feparm)
Set the parameter objects.
Definition: vfetk.c:615
VPUBLIC double Vfetk_PDE_Fu_v(PDE *thee, int key, double V[], double dV[][VAPBS_DIM])
This is the weak form of the PBE; i.e. the strong form integrated with a test function to give: wher...
Definition: vfetk.c:1703
Vfetk_NsolvType nkey
Definition: vfetk.h:190
VEXTERNC double Vacc_vdwAcc(Vacc *thee, double center[VAPBS_DIM])
Report van der Waals accessibility.
int level
Definition: vfetk.h:200
VPUBLIC double Vpbe_getXkappa(Vpbe *thee)
Get Debye-Huckel parameter.
Definition: vpbe.c:134
Vfetk_GuessType gues
Definition: vfetk.h:193
double ionConc[MAXION]
Definition: vpbe.h:104
enum eVhal_PBEType Vhal_PBEType
Declaration of the Vhal_PBEType type as the Vhal_PBEType enum.
Definition: vhal.h:151
VPUBLIC void Vcsm_dtor(Vcsm **thee)
Destroy Vcsm object.
Definition: vcsm.c:292
Definition: vhal.h:141
Definition: vhal.h:277
double ionRadii[MAXION]
Definition: vfetk.h:243
Charge-simplex map class.
Definition: vcsm.h:89
VEXTERNC void Vfetk_PDE_initElement(PDE *thee, int elementType, int chart, double tvx[][VAPBS_DIM], void *data)
Do once-per-element initialization.
Vgreen * green
Definition: vfetk.h:234
Definition: vhal.h:279
VPUBLIC double * Vfetk_getSolution(Vfetk *thee, int *length)
Create an array containing the solution (electrostatic potential in units of ) at the finest mesh lev...
Definition: vfetk.c:641
VPUBLIC double Vfetk_PDE_Ju(PDE *thee, int key)
Energy functional. This returns the energy (less delta function terms) in the form: for a 1:1 electr...
Definition: vfetk.c:2003
double F
Definition: vfetk.h:229
Definition: vhal.h:312
Definition: vhal.h:310
VV * verts[4]
Definition: vfetk.h:239
Definition: vhal.h:103
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition: valist.c:115
double solventRadius
Definition: vpbe.h:95
double ionQ[MAXION]
Definition: vfetk.h:242
Contains declarations for class Vfetk.
Parameter structure for FEM-specific variables from input files.
Definition: femparm.h:133
Vhal_PBEType type
Definition: vfetk.h:199
double B
Definition: vfetk.h:230
VPUBLIC double Vfetk_qfEnergy(Vfetk *thee, int color)
Get the "fixed charge" contribution to the electrostatic energy.
Definition: vfetk.c:732
VPUBLIC int Vgreen_coulombD(Vgreen *thee, int npos, double *x, double *y, double *z, double *pot, double *gradx, double *grady, double *gradz)
Get gradient of the Coulomb's Law Green's function (solution to Laplace's equation) integrated over t...
Definition: vgreen.c:362
PDE * pde
Definition: vfetk.h:184
double zks2
Definition: vfetk.h:245
Definition: vhal.h:275
VPUBLIC double Vfetk_dqmEnergy(Vfetk *thee, int color)
Get the "mobile charge" and "polarization" contributions to the electrostatic energy.
Definition: vfetk.c:842
#define Vunit_kb
Boltzmann constant.
Definition: vunit.h:96
Definition: vhal.h:273
Vcsm * csm
Definition: vfetk.h:186
Aprx * aprx
Definition: vfetk.h:183
double delta
Definition: vfetk.h:250
#define VAPBS_DIM
Our dimension.
Definition: vhal.h:402
VPUBLIC void Vfetk_PDE_Fu(PDE *thee, int key, double F[])
Evaluate strong form of PBE. For interior points, this is: where is the (possibly nonlinear) mobile...
Definition: vfetk.c:1695
Definition: vhal.h:311
VEXTERNC void Vfetk_PDE_bisectEdge(int dim, int dimII, int edgeType, int chart[], double vx[][VAPBS_DIM])
Define the way manifold edges are bisected.
double zkappa2
Definition: vfetk.h:244
VPUBLIC int Vcsm_getAtomIndex(Vcsm *thee, int iatom, int isimp)
Get ID of particular atom in a simplex.
Definition: vcsm.c:88
Parameter structure for PBE variables from input files.
Definition: pbeparm.h:117
Definition: vfetk.h:158
Definition: vhal.h:140
double dW[VAPBS_DIM]
Definition: vfetk.h:222
VPUBLIC void Vfetk_PDE_dtor2(PDE *thee)
FORTRAN stub: destroys FEtk PDE object.
Definition: vfetk.c:1246
AM * am
Definition: vfetk.h:182
VPUBLIC Gem * Vfetk_getGem(Vfetk *thee)
Get a pointer to the Gem (grid manager) object.
Definition: vfetk.c:490
Definition: vhal.h:281
VPUBLIC double Vacc_molAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility.
Definition: vacc.c:608
VPUBLIC double Vpbe_getMaxIonRadius(Vpbe *thee)
Get maximum radius of ion species.
Definition: vpbe.c:127
double DB
Definition: vfetk.h:231
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition: vatom.c:63
VPUBLIC void Vfetk_PDE_initFace(PDE *thee, int faceType, int chart, double tnvec[])
Do once-per-face initialization.
Definition: vfetk.c:1559
VPUBLIC unsigned long int Vcsm_memChk(Vcsm *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition: vcsm.c:129
VPUBLIC double Vpbe_getZmagic(Vpbe *thee)
Get charge scaling factor.
Definition: vpbe.c:155
enum eVdata_Type Vdata_Type
Declaration of the Vdata_Type type as the Vdata_Type enum.
Definition: vhal.h:302
double Fu_v
Definition: vfetk.h:248
Contains public data members for Vatom class/module.
Definition: vatom.h:84
VPUBLIC void Vgreen_dtor(Vgreen **thee)
Destruct the Green's function oracle.
Definition: vgreen.c:192
VPUBLIC void Vfetk_PDE_u_T(PDE *thee, int type, int chart, double txq[], double F[])
Evaluate the "true solution" at the given point for comparison with the numerical solution.
Definition: vfetk.c:1886
Container class for list of atom objects.
Definition: valist.h:78
Valist * alist
Definition: vpbe.h:88
#define Vunit_eps0
Vacuum permittivity.
Definition: vunit.h:108
VEXTERNC void Vfetk_PDE_oneChart(int dim, int dimII, int objType, int chart[], double vx[][VAPBS_DIM], int dimV)
Unify the chart for different coordinate systems – a no-op for us.
double srad
Definition: pbeparm.h:152
VPUBLIC Vrc_Codes Vfetk_genCube(Vfetk *thee, double center[3], double length[3], Vfetk_MeshLoad meshType)
Construct a rectangular mesh (in the current Vfetk object)
Definition: vfetk.c:885
Contains public data members for Vfetk class/module.
Definition: vfetk.h:176
Definition: vfetk.h:88
VPUBLIC int Vcsm_getNumberAtoms(Vcsm *thee, int isimp)
Get number of atoms associated with a simplex.
Definition: vcsm.c:69
double ltol
Definition: vfetk.h:189
VPUBLIC int Vfetk_PDE_ctor2(PDE *thee, Vfetk *fetk)
Intializes the FEtk PDE object.
Definition: vfetk.c:1187
VPUBLIC void Vcsm_init(Vcsm *thee)
Initialize charge-simplex map with mesh and atom data.
Definition: vcsm.c:170
Vpbe * pbe
Definition: vfetk.h:185
int initGreen
Definition: vfetk.h:235
VPUBLIC double Vfetk_energy(Vfetk *thee, int color, int nonlin)
Return the total electrostatic energy.
Definition: vfetk.c:693
enum eVfetk_MeshLoad Vfetk_MeshLoad
Declare FEMparm_GuessType type.
Definition: vfetk.h:114
#define VAPBS_NVS
Number of vertices per simplex (hard-coded to 3D)
Definition: vhal.h:397
double dU[MAXV][VAPBS_DIM]
Definition: vfetk.h:220
VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type)
Fill an array with the specified data.
Definition: vfetk.c:2299
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
int pjac
Definition: vfetk.h:195
VPUBLIC double Vacc_splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad)
Report spline-based accessibility.
Definition: vacc.c:528
VPUBLIC int Vfetk_write(Vfetk *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, Bvec *vec, Vdata_Format format)
Write out data.
Definition: vfetk.c:2464
#define Vunit_ec
Charge of an electron in C.
Definition: vunit.h:92
VPUBLIC unsigned long int Vfetk_memChk(Vfetk *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition: vfetk.c:867
double maxIonRadius
Definition: vpbe.h:100
Vfetk * fetk
Definition: vfetk.h:233
Gem * gm
Definition: vcsm.h:94
double ntol
Definition: vfetk.h:192
double W
Definition: vfetk.h:221
VPUBLIC double Vpbe_getSolventDiel(Vpbe *thee)
Get solvent dielectric constant.
Definition: vpbe.c:113