APBS 3.0.0
Loading...
Searching...
No Matches
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 */
74VPRIVATE 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 */
86VPRIVATE 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 */
93VPRIVATE char *diriCubeString =
94"mcsf_begin=1;\n\
95\n\
96dim=3;\n\
97dimii=3;\n\
98vertices=8;\n\
99simplices=6;\n\
100\n\
101vert=[\n\
1020 0 -0.5 -0.5 -0.5\n\
1031 0 0.5 -0.5 -0.5\n\
1042 0 -0.5 0.5 -0.5\n\
1053 0 0.5 0.5 -0.5\n\
1064 0 -0.5 -0.5 0.5\n\
1075 0 0.5 -0.5 0.5\n\
1086 0 -0.5 0.5 0.5\n\
1097 0 0.5 0.5 0.5\n\
110];\n\
111\n\
112simp=[\n\
1130 0 0 0 1 0 1 0 5 1 2\n\
1141 0 0 0 1 1 0 0 5 2 4\n\
1152 0 0 0 1 0 1 1 5 3 2\n\
1163 0 0 0 1 0 1 3 5 7 2\n\
1174 0 0 1 1 0 0 2 5 7 6\n\
1185 0 0 1 1 0 0 2 5 6 4\n\
119];\n\
120\n\
121mcsf_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 */
130VPRIVATE char *neumCubeString =
131"mcsf_begin=1;\n\
132\n\
133dim=3;\n\
134dimii=3;\n\
135vertices=8;\n\
136simplices=6;\n\
137\n\
138vert=[\n\
1390 0 -0.5 -0.5 -0.5\n\
1401 0 0.5 -0.5 -0.5\n\
1412 0 -0.5 0.5 -0.5\n\
1423 0 0.5 0.5 -0.5\n\
1434 0 -0.5 -0.5 0.5\n\
1445 0 0.5 -0.5 0.5\n\
1456 0 -0.5 0.5 0.5\n\
1467 0 0.5 0.5 0.5\n\
147];\n\
148\n\
149simp=[\n\
1500 0 0 0 2 0 2 0 5 1 2\n\
1511 0 0 0 2 2 0 0 5 2 4\n\
1522 0 0 0 2 0 2 1 5 3 2\n\
1533 0 0 0 2 0 2 3 5 7 2\n\
1544 0 0 2 2 0 0 2 5 7 6\n\
1555 0 0 2 2 0 0 2 5 6 4\n\
156];\n\
157\n\
158mcsf_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 */
170VPRIVATE 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 */
180VPRIVATE 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 */
194VPRIVATE 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 */
213VPRIVATE 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 */
230VPRIVATE 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 */
250VPRIVATE 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 */
270VPRIVATE 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 */
292VPRIVATE 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 */
316VPRIVATE 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 */
354VPRIVATE 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 */
366VPRIVATE 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 */
386VPRIVATE 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};
395VPRIVATE 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};
402VPRIVATE 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};
409VPRIVATE 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 */
437VPRIVATE int dim_3DP1 = VAPBS_NVS;
438VPRIVATE 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};
446VPRIVATE 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};
454VPRIVATE 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};
462VPRIVATE 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 */
476VPRIVATE const int P_DEG=1;
477
478/*
479 * @brief Another Holst variable
480 * @ingroup Vfetk
481 * @author Mike Holst */
482VPRIVATE int numP;
483VPRIVATE double c[VMAXP][VMAXP];
484VPRIVATE double cx[VMAXP][VMAXP];
485VPRIVATE double cy[VMAXP][VMAXP];
486VPRIVATE double cz[VMAXP][VMAXP];
487
488#if !defined(VINLINE_VFETK)
489
490VPUBLIC Gem* Vfetk_getGem(Vfetk *thee) {
491
492 VASSERT(thee != VNULL);
493 return thee->gm;
494
495}
496
497VPUBLIC AM* Vfetk_getAM(Vfetk *thee) {
498
499 VASSERT(thee != VNULL);
500 return thee->am;
501}
502
503VPUBLIC Vpbe* Vfetk_getVpbe(Vfetk *thee) {
504
505 VASSERT(thee != VNULL);
506 return thee->pbe;
507
508}
509
510VPUBLIC Vcsm* Vfetk_getVcsm(Vfetk *thee) {
511
512 VASSERT(thee != VNULL);
513 return thee->csm;
514
515}
516
517VPUBLIC int Vfetk_getAtomColor(Vfetk *thee,
518 int iatom
519 ) {
520
521 int natoms;
522
523 VASSERT(thee != VNULL);
524
526 VASSERT(iatom < natoms);
527
528 return Vatom_getPartID(Valist_getAtom(Vpbe_getValist(thee->pbe), iatom));
529}
530#endif /* if !defined(VINLINE_VFETK) */
531
532VPUBLIC 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
545VPUBLIC 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
615VPUBLIC 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
625VPUBLIC 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
633VPUBLIC 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
641VPUBLIC 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
693VPUBLIC 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
732VPUBLIC 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
774VPRIVATE 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
842VPUBLIC double Vfetk_dqmEnergy(Vfetk *thee,
843 int color) {
844
845 return AM_evalJ(thee->am);
846
847}
848
849VPUBLIC 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
867VPUBLIC 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
885VPUBLIC 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
980VPUBLIC 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
1026VPUBLIC 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
1176VPUBLIC 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
1187VPUBLIC 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
1231VPUBLIC 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
1246VPUBLIC void Vfetk_PDE_dtor2(PDE *thee) {
1247 var.fetk = VNULL;
1248}
1249
1250VPRIVATE 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
1279VPRIVATE 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
1326VPRIVATE 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
1369VPRIVATE 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
1404VPRIVATE 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
1440VPRIVATE 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
1508VPUBLIC 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
1528VPUBLIC 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
1559VPUBLIC 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
1571VPUBLIC 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
1695VPUBLIC 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
1703VPUBLIC 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
1741VPUBLIC 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
1779#define VATOMMAX 1000000
1780VPUBLIC 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 \
1842and Gem_pointInSimplex detected misplaced point charge!\n");
1843 Vnm_print(2, "delta: I think you have problems: \
1844phi = {");
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
1867VPUBLIC 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
1886VPUBLIC 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
1902VPUBLIC 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
1912VPUBLIC void Vfetk_PDE_mapBoundary(int dim, int dimII, int vertexType,
1913 int chart, double vx[3]) {
1914
1915}
1916
1917VPUBLIC 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
1998VPUBLIC void Vfetk_PDE_oneChart(int dim, int dimII, int objType, int chart[],
1999 double vx[][3], int dimV) {
2000
2001}
2002
2003VPUBLIC 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
2078VPUBLIC 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
2095VPRIVATE 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
2125VPRIVATE 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
2140VPUBLIC 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
2203VPUBLIC 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
2217VPRIVATE 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
2236VPRIVATE 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
2255VPUBLIC 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
2299VPUBLIC 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
2464VPUBLIC 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}
VPUBLIC double Vacc_molAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility.
Definition vacc.c:608
VPUBLIC double Vacc_splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad)
Report spline-based accessibility.
Definition vacc.c:528
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_getPartID(Vatom *thee)
Get partition ID.
Definition vatom.c:70
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
VPUBLIC void Vatom_setPartID(Vatom *thee, int partID)
Set partition ID.
Definition vatom.c:77
VPUBLIC double Vcap_exp(double x, int *ichop)
Provide a capped exp() function.
Definition vcap.c:59
VPUBLIC int Vcsm_getNumberAtoms(Vcsm *thee, int isimp)
Get number of atoms associated with a simplex.
Definition vcsm.c:69
VPUBLIC int Vcsm_getAtomIndex(Vcsm *thee, int iatom, int isimp)
Get ID of particular atom in a simplex.
Definition vcsm.c:88
VPUBLIC Vcsm * Vcsm_ctor(Valist *alist, Gem *gm)
Construct Vcsm object.
Definition vcsm.c:136
VPUBLIC SS * Vcsm_getSimplex(Vcsm *thee, int isimp, int iatom)
Get particular simplex associated with an atom.
Definition vcsm.c:109
VPUBLIC int Vcsm_update(Vcsm *thee, SS **simps, int num)
Update the charge-simplex and simplex-charge maps after refinement.
Definition vcsm.c:326
VPUBLIC Vatom * Vcsm_getAtom(Vcsm *thee, int iatom, int isimp)
Get particular atom associated with a simplex.
Definition vcsm.c:77
VPUBLIC void Vcsm_init(Vcsm *thee)
Initialize charge-simplex map with mesh and atom data.
Definition vcsm.c:170
VPUBLIC int Vcsm_getNumberSimplices(Vcsm *thee, int iatom)
Get number of simplices associated with an atom.
Definition vcsm.c:99
VEXTERNC void Gem_setExternalUpdateFunction(Gem *thee, void(*externalUpdate)(SS **simps, int num))
External function for FEtk Gem class to use during mesh refinement.
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 void Vcsm_dtor(Vcsm **thee)
Destroy Vcsm object.
Definition vcsm.c:292
enum eVfetk_MeshLoad Vfetk_MeshLoad
Declare FEMparm_GuessType type.
Definition vfetk.h:114
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
VPUBLIC double Vfetk_PDE_Ju(PDE *thee, int key)
Energy functional. This returns the energy (less delta function terms) in the form:
Definition vfetk.c:2003
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
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_PDE_initFace(PDE *thee, int faceType, int chart, double tnvec[])
Do once-per-face initialization.
Definition vfetk.c:1559
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
VPUBLIC void Vfetk_dtor(Vfetk **thee)
Object destructor.
Definition vfetk.c:625
VPUBLIC void Vfetk_PDE_dtor(PDE **thee)
Destroys FEtk PDE object.
Definition vfetk.c:1231
VPUBLIC int Vfetk_getAtomColor(Vfetk *thee, int iatom)
Get the partition information for a particular atom.
Definition vfetk.c:517
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
VPUBLIC void Bmat_printHB(Bmat *thee, char *fname)
Writes a Bmat to disk in Harwell-Boeing sparse matrix format.
Definition vfetk.c:1026
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:
Definition vfetk.c:1703
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
VPUBLIC Gem * Vfetk_getGem(Vfetk *thee)
Get a pointer to the Gem (grid manager) object.
Definition vfetk.c:490
struct sVfetk Vfetk
Declaration of the Vfetk class as the Vfetk structure.
Definition vfetk.h:207
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 Vfetk * Vfetk_ctor(Vpbe *pbe, Vhal_PBEType type)
Constructor for Vfetk object.
Definition vfetk.c:532
VPUBLIC double Vfetk_dqmEnergy(Vfetk *thee, int color)
Get the "mobile charge" and "polarization" contributions to the electrostatic energy.
Definition vfetk.c:842
#define VATOMMAX
Maximum number of atoms associated with a vertex.
Definition vfetk.c:1779
VPUBLIC unsigned long int Vfetk_memChk(Vfetk *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition vfetk.c:867
#define VRINGMAX
Maximum number of simplices in a simplex ring.
Definition vfetk.c:1776
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 void Vfetk_PDE_initAssemble(PDE *thee, int ip[], double rp[])
Do once-per-assembly initialization.
Definition vfetk.c:1508
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
VPUBLIC PDE * Vfetk_PDE_ctor(Vfetk *fetk)
Constructs the FEtk PDE object.
Definition vfetk.c:1176
VPUBLIC void Vfetk_PDE_dtor2(PDE *thee)
FORTRAN stub: destroys FEtk PDE object.
Definition vfetk.c:1246
VPUBLIC int Vfetk_PDE_ctor2(PDE *thee, Vfetk *fetk)
Intializes the FEtk PDE object.
Definition vfetk.c:1187
VPUBLIC void Vfetk_dtor2(Vfetk *thee)
FORTRAN stub object destructor.
Definition vfetk.c:633
VPUBLIC double Vfetk_qfEnergy(Vfetk *thee, int color)
Get the "fixed charge" contribution to the electrostatic energy.
Definition vfetk.c:732
VPUBLIC Vcsm * Vfetk_getVcsm(Vfetk *thee)
Get a pointer to the Vcsm (charge-simplex map) object.
Definition vfetk.c:510
VPUBLIC double Vfetk_energy(Vfetk *thee, int color, int nonlin)
Return the total electrostatic energy.
Definition vfetk.c:693
VPUBLIC Vpbe * Vfetk_getVpbe(Vfetk *thee)
Get a pointer to the Vpbe (PBE manager) object.
Definition vfetk.c:503
VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type)
Fill an array with the specified data.
Definition vfetk.c:2299
VPUBLIC void Vfetk_PDE_Fu(PDE *thee, int key, double F[])
Evaluate strong form of PBE. For interior points, this is:
Definition vfetk.c:1695
VPUBLIC int Vfetk_ctor2(Vfetk *thee, Vpbe *pbe, Vhal_PBEType type)
FORTRAN stub constructor for Vfetk object.
Definition vfetk.c:545
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
VPUBLIC void Vfetk_setAtomColors(Vfetk *thee)
Transfer color (partition ID) information frmo a partitioned mesh to the atoms.
Definition vfetk.c:849
VPUBLIC AM * Vfetk_getAM(Vfetk *thee)
Get a pointer to the AM (algebra manager) object.
Definition vfetk.c:497
VPUBLIC void Vfetk_dumpLocalVar()
Debugging routine to print out local variables used by PDE object.
Definition vfetk.c:2255
VPUBLIC void Vfetk_setParameters(Vfetk *thee, PBEparm *pbeparm, FEMparm *feparm)
Set the parameter objects.
Definition vfetk.c:615
@ VLT_MG
Definition vfetk.h:88
@ VML_DIRICUBE
Definition vfetk.h:105
@ VML_NEUMCUBE
Definition vfetk.h:106
@ VML_EXTERNAL
Definition vfetk.h:107
@ VNT_NEW
Definition vfetk.h:122
@ VPT_MG
Definition vfetk.h:158
@ VGT_ZERO
Definition vfetk.h:139
VPUBLIC void Vgreen_dtor(Vgreen **thee)
Destruct the Green's function oracle.
Definition vgreen.c:192
VPUBLIC Vgreen * Vgreen_ctor(Valist *alist)
Construct the Green's function oracle.
Definition vgreen.c:156
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
enum eVdata_Format Vdata_Format
Declaration of the Vdata_Format type as the Vdata_Format enum.
Definition vhal.h:323
#define VAPBS_NVS
Number of vertices per simplex (hard-coded to 3D)
Definition vhal.h:397
enum eVsurf_Meth Vsurf_Meth
Declaration of the Vsurf_Meth type as the Vsurf_Meth enum.
Definition vhal.h:133
enum eVhal_PBEType Vhal_PBEType
Declaration of the Vhal_PBEType type as the Vhal_PBEType enum.
Definition vhal.h:151
enum eVdata_Type Vdata_Type
Declaration of the Vdata_Type type as the Vdata_Type enum.
Definition vhal.h:302
#define VAPBS_DIM
Our dimension.
Definition vhal.h:402
@ VSM_MOLSMOOTH
Definition vhal.h:107
@ VSM_MOL
Definition vhal.h:103
@ VSM_SPLINE
Definition vhal.h:109
@ VDT_SSPL
Definition vhal.h:275
@ VDT_CHARGE
Definition vhal.h:270
@ VDT_NDENS
Definition vhal.h:284
@ VDT_KAPPA
Definition vhal.h:294
@ VDT_IVDW
Definition vhal.h:279
@ VDT_SMOL
Definition vhal.h:273
@ VDT_DIELX
Definition vhal.h:288
@ VDT_EDENS
Definition vhal.h:282
@ VDT_POT
Definition vhal.h:271
@ VDT_DIELZ
Definition vhal.h:292
@ VDT_DIELY
Definition vhal.h:290
@ VDT_QDENS
Definition vhal.h:286
@ VDT_VDW
Definition vhal.h:277
@ VDT_LAP
Definition vhal.h:281
@ VDF_AVS
Definition vhal.h:312
@ VDF_UHBD
Definition vhal.h:311
@ VDF_DX
Definition vhal.h:310
@ VRC_FAILURE
Definition vhal.h:69
@ VRC_SUCCESS
Definition vhal.h:70
@ PBE_LRPBE
Definition vhal.h:142
@ PBE_SMPBE
Definition vhal.h:144
@ PBE_LPBE
Definition vhal.h:140
@ PBE_NPBE
Definition vhal.h:141
VPUBLIC double Vpbe_getSoluteDiel(Vpbe *thee)
Get solute dielectric constant.
Definition vpbe.c:99
VPUBLIC double Vpbe_getZkappa2(Vpbe *thee)
Get modified squared Debye-Huckel parameter.
Definition vpbe.c:148
VPUBLIC double Vpbe_getXkappa(Vpbe *thee)
Get Debye-Huckel parameter.
Definition vpbe.c:134
VPUBLIC double Vpbe_getZmagic(Vpbe *thee)
Get charge scaling factor.
Definition vpbe.c:155
VPUBLIC double Vpbe_getSolventDiel(Vpbe *thee)
Get solvent dielectric constant.
Definition vpbe.c:113
VPUBLIC double Vpbe_getBulkIonicStrength(Vpbe *thee)
Get bulk ionic strength.
Definition vpbe.c:84
VPUBLIC double Vpbe_getMaxIonRadius(Vpbe *thee)
Get maximum radius of ion species.
Definition vpbe.c:127
VPUBLIC Valist * Vpbe_getValist(Vpbe *thee)
Get atom list.
Definition vpbe.c:69
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 double Vpbe_getTemperature(Vpbe *thee)
Get temperature.
Definition vpbe.c:91
#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
Parameter structure for FEM-specific variables from input files.
Definition femparm.h:133
double targetRes
Definition femparm.h:160
Parameter structure for PBE variables from input files.
Definition pbeparm.h:117
double swin
Definition pbeparm.h:154
Vsurf_Meth srfm
Definition pbeparm.h:150
double srad
Definition pbeparm.h:152
Oracle for solvent- and ion-accessibility around a biomolecule.
Definition vacc.h:108
Container class for list of atom objects.
Definition valist.h:78
Contains public data members for Vatom class/module.
Definition vatom.h:84
Charge-simplex map class.
Definition vcsm.h:89
Vfetk LocalVar subclass.
Definition vfetk.h:215
double u_D
Definition vfetk.h:251
double vx[4][VAPBS_DIM]
Definition vfetk.h:217
double ionRadii[MAXION]
Definition vfetk.h:243
double U[MAXV]
Definition vfetk.h:219
double xq[VAPBS_DIM]
Definition vfetk.h:218
double zks2
Definition vfetk.h:245
double nvec[VAPBS_DIM]
Definition vfetk.h:216
double zkappa2
Definition vfetk.h:244
double Fu_v
Definition vfetk.h:248
double u_T
Definition vfetk.h:252
double delta
Definition vfetk.h:250
double d2W
Definition vfetk.h:223
double DB
Definition vfetk.h:231
VV * verts[4]
Definition vfetk.h:239
double ionstr
Definition vfetk.h:246
double dW[VAPBS_DIM]
Definition vfetk.h:222
Vgreen * green
Definition vfetk.h:234
double jumpDiel
Definition vfetk.h:232
Vfetk * fetk
Definition vfetk.h:233
double ionConc[MAXION]
Definition vfetk.h:241
double dU[MAXV][VAPBS_DIM]
Definition vfetk.h:220
double DFu_wv
Definition vfetk.h:249
double ionQ[MAXION]
Definition vfetk.h:242
Contains public data members for Vfetk class/module.
Definition vfetk.h:176
double ltol
Definition vfetk.h:189
AM * am
Definition vfetk.h:182
Vpbe * pbe
Definition vfetk.h:185
Vhal_PBEType type
Definition vfetk.h:199
PDE * pde
Definition vfetk.h:184
Gem * gm
Definition vfetk.h:179
int lmax
Definition vfetk.h:188
double ntol
Definition vfetk.h:192
Vfetk_LsolvType lkey
Definition vfetk.h:187
Vfetk_GuessType gues
Definition vfetk.h:193
int nmax
Definition vfetk.h:191
Vfetk_PrecType lprec
Definition vfetk.h:194
Aprx * aprx
Definition vfetk.h:183
int level
Definition vfetk.h:200
Vcsm * csm
Definition vfetk.h:186
Vfetk_NsolvType nkey
Definition vfetk.h:190
Vmem * vmem
Definition vfetk.h:178
int pjac
Definition vfetk.h:195
FEMparm * feparm
Definition vfetk.h:198
PBEparm * pbeparm
Definition vfetk.h:197
Contains public data members for Vpbe class/module.
Definition vpbe.h:84
Vacc * acc
Definition vpbe.h:90
int numIon
Definition vpbe.h:103
double maxIonRadius
Definition vpbe.h:100
double solventRadius
Definition vpbe.h:95
Valist * alist
Definition vpbe.h:88
double ionConc[MAXION]
Definition vpbe.h:104
double ionQ[MAXION]
Definition vpbe.h:106
Contains declarations for class Vfetk.