74VPRIVATE
double Vfetk_qfEnergyAtom(
93VPRIVATE
char *diriCubeString =
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\
130VPRIVATE
char *neumCubeString =
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\
170VPRIVATE
double diel();
180VPRIVATE
double ionacc();
194VPRIVATE
double smooth(
213VPRIVATE
double debye_U(
230VPRIVATE
double debye_Udiff(
250VPRIVATE
void coulomb(
270VPRIVATE
void init_2DP1(
292VPRIVATE
void init_3DP1(
316VPRIVATE
void setCoef(
354VPRIVATE
void polyEval(
366VPRIVATE
int dim_2DP1 = 3;
386VPRIVATE
int lgr_2DP1[3][VMAXP] = {
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 }
395VPRIVATE
int lgr_2DP1x[3][VMAXP] = {
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 }
402VPRIVATE
int lgr_2DP1y[3][VMAXP] = {
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 }
409VPRIVATE
int lgr_2DP1z[3][VMAXP] = {
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 }
438VPRIVATE
int lgr_3DP1[
VAPBS_NVS][VMAXP] = {
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 }
446VPRIVATE
int lgr_3DP1x[
VAPBS_NVS][VMAXP] = {
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 }
454VPRIVATE
int lgr_3DP1y[
VAPBS_NVS][VMAXP] = {
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 }
462VPRIVATE
int lgr_3DP1z[
VAPBS_NVS][VMAXP] = {
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 },
476VPRIVATE
const int P_DEG=1;
483VPRIVATE
double c[VMAXP][VMAXP];
484VPRIVATE
double cx[VMAXP][VMAXP];
485VPRIVATE
double cy[VMAXP][VMAXP];
486VPRIVATE
double cz[VMAXP][VMAXP];
488#if !defined(VINLINE_VFETK)
492 VASSERT(thee != VNULL);
499 VASSERT(thee != VNULL);
505 VASSERT(thee != VNULL);
512 VASSERT(thee != VNULL);
523 VASSERT(thee != VNULL);
526 VASSERT(iatom < natoms);
538 thee = (
Vfetk*)Vmem_malloc(VNULL, 1,
sizeof(
Vfetk) );
539 VASSERT(thee != VNULL);
554 VASSERT(pbe != VNULL);
556 VASSERT(pbe->
alist != VNULL);
557 VASSERT(pbe->
acc != VNULL);
563 thee->
vmem = Vmem_ctor(
"APBS::VFETK");
566 Vnm_print(0,
"Vfetk_ctor2: Constructing PDE...\n");
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");
572 Vnm_print(0,
"Vfetk_ctor2: Constructing Aprx...\n");
580 thee->
lmax = 1000000;
584 thee->
nmax = 1000000;
603 for (i=0; i<var.
nion; i++) {
620 VASSERT(thee != VNULL);
626 if ((*thee) != VNULL) {
635 AM_dtor(&(thee->
am));
636 Aprx_dtor(&(thee->
aprx));
638 Vmem_dtor(&(thee->
vmem));
650 VASSERT(thee != VNULL);
655 Bvec_copy(am->w0, am->u);
657 Bvec_axpy(am->w0, am->ud, 1.);
659 solution = Bvec_addr(am->w0);
661 *length = Bvec_numRT(am->w0);
664 VASSERT(1 == Bvec_numB(am->w0));
667 theAnswer = (
double*)Vmem_malloc(VNULL, *length,
sizeof(
double));
668 VASSERT(theAnswer != VNULL);
669 for (i=0; i<(*length); i++) theAnswer[i] = solution[i];
702 double totEnergy = 0.0,
706 VASSERT(thee != VNULL);
709 Vnm_print(0,
"Vfetk_energy: calculating full PBE energy\n");
710 Vnm_print(0,
"Vfetk_energy: bulk ionic strength = %g M\n",
713 Vnm_print(0,
"Vfetk_energy: dqmEnergy = %g kT\n", dqmEnergy);
715 Vnm_print(0,
"Vfetk_energy: qfEnergy = %g kT\n", qfEnergy);
717 totEnergy = qfEnergy - dqmEnergy;
719 Vnm_print(0,
"Vfetk_energy: calculating only q-phi energy\n");
721 Vnm_print(0,
"Vfetk_energy: dqmEnergy = %g kT (NOT USED)\n", dqmEnergy);
723 Vnm_print(0,
"Vfetk_energy: qfEnergy = %g kT\n", qfEnergy);
724 totEnergy = 0.5*qfEnergy;
743 VASSERT(thee != VNULL);
749 VASSERT(sol != VNULL);
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");
761 for (iatom=0; iatom<natoms; iatom++) {
763 energy = energy + Vfetk_qfEnergyAtom(thee, iatom, color, sol);
768 Vmem_free(VNULL, nsol,
sizeof(
double), (
void **)&sol);
774VPRIVATE
double Vfetk_qfEnergyAtom(
802 usingColor = (color >= 0);
804 if (usingColor && (icolor<0)) {
805 Vnm_print(2,
"Vfetk_qfEnergy: Atom colors not set!\n");
810 if ((icolor==color) || (!usingColor)) {
817 for (isimp=0; isimp<nsimps; isimp++) {
823 if ((SS_chart(simp)==color)||(color<0)) {
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);
845 return AM_evalJ(thee->
am);
856 VASSERT(thee != VNULL);
859 for (i=0; i<natoms; i++) {
871 if (thee == VNULL)
return 0;
873 memUse = memUse +
sizeof(
Vfetk);
891 VASSERT(thee != VNULL);
903 *iohost =
"localhost",
911 VASSERT(am != VNULL);
913 VASSERT(gm != VNULL);
920 bufsize = strlen(diriCubeString);
921 VASSERT( bufsize <= VMAX_BUFSIZE );
922 strncpy(buf, diriCubeString, VMAX_BUFSIZE);
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);
932 Vnm_print(2,
"Vfetk_genCube: Got request for external mesh!\n");
933 Vnm_print(2,
"Vfetk_genCube: How did we get here?\n");
936 Vnm_print(2,
"Vfetk_genCube: Unknown mesh type (%d)\n", meshType);
940 VASSERT( VNULL != (sock=Vio_socketOpen(key,iodev,iofmt,iohost,iofile)) );
941 Vio_bufTake(sock, buf, bufsize);
942 AM_read(am, skey, sock);
945 Vio_connectFree(sock);
952 for (i=0; i<Gem_numVV(gm); i++) {
954 for (j=0; j<3; j++) {
957 VV_setCoord(vx, j, x);
962 for (i=0; i<Gem_numVV(gm); i++) {
964 for (j=0; j<3; j++) {
967 VV_setCoord(vx, j, x);
994 Vnm_print(2,
"Vfetk_loadMesh: Got NULL socket!\n");
997 AM_read(thee->
am, skey, sock);
998 Vio_connectFree(sock);
1009 Vnm_print(2,
"Vfetk_loadMesh: unrecognized mesh type (%d)!\n",
1015 Vnm_print(0,
"Vfetk_ctor2: Constructing Vcsm...\n");
1019 VASSERT(thee->
csm != VNULL);
1038 char mmkey[] = {
"8charkey"};
1039 int totc = 0, ptrc = 0, indc = 0, valc = 0;
1040 char mxtyp[] = {
"RUA"};
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"};
1048 VASSERT( thee->numB == 1 );
1049 Ablock = thee->AD[0][0];
1051 VASSERT( Mat_format( Ablock ) == DRC_FORMAT );
1053 pqsym = Mat_sym( Ablock );
1055 if ( pqsym == IS_SYM ) {
1057 }
else if ( pqsym == ISNOT_SYM ) {
1063 nrow = Bmat_numRT( thee );
1064 ncol = Bmat_numCT( thee );
1065 numZ = Bmat_numZT( thee );
1067 nrowdigits = (int) (log( nrow )/log( 10 )) + 1;
1068 numZdigits = (int) (log( numZ )/log( 10 )) + 1;
1070 nptrline = (int) ( 80 / (numZdigits + 1) );
1071 nindline = (int) ( 80 / (nrowdigits + 1) );
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);
1078 ptrc = (int) ( ( (ncol + 1) - 1 ) / nptrline ) + 1;
1079 indc = (int) ( (numZ - 1) / nindline ) + 1;
1080 valc = (int) ( (numZ - 1) / nvalline ) + 1;
1082 totc = ptrc + indc + valc;
1084 sprintf( mmtitle,
"Sparse '%s' Matrix - Harwell-Boeing Format - '%s'",
1085 thee->name, fname );
1089 fp = fopen( fname,
"w" );
1091 Vnm_print(2,
"Bmat_printHB: Ouch couldn't open file <%s>\n",fname);
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" );
1108 if ( pqsym == IS_SYM ) {
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" );
1119 if ( ( (ncol+1) % nptrline ) != 0 ) {
1120 fprintf( fp,
"\n" );
1126 for (i=0; i<ncol; i++) {
1127 fprintf( fp, indfmtstr, i+1);
1128 if ( ( (j+1) % nindline ) == 0 ) {
1129 fprintf( fp,
"\n" );
1132 for (jj=IA[i]; jj<IA[i+1]; jj++) {
1133 fprintf( fp, indfmtstr, JA[jj] + 1 );
1134 if ( ( (j+1) % nindline ) == 0 ) {
1135 fprintf( fp,
"\n" );
1141 if ( ( j % nindline ) != 0 ) {
1142 fprintf( fp,
"\n" );
1148 for (i=0; i<ncol; i++) {
1149 fprintf( fp, valfmtstr, D[i] );
1150 if ( ( (j+1) % nvalline ) == 0 ) {
1151 fprintf( fp,
"\n" );
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" );
1163 if ( ( j % nvalline ) != 0 ) {
1164 fprintf( fp,
"\n" );
1180 thee = (PDE*)Vmem_malloc(fetk->
vmem, 1,
sizeof(PDE));
1181 VASSERT(thee != VNULL);
1193 if (thee == VNULL) {
1194 Vnm_print(2,
"Vfetk_PDE_ctor2: Got NULL thee!\n");
1203 thee->initElement = Vfetk_PDE_initElement;
1205 thee->initPoint = Vfetk_PDE_initPoint;
1208 thee->DFu_wv = Vfetk_PDE_DFu_wv;
1214 thee->sym[0][0] = 1;
1216 for (i=0; i<VMAX_BDTYPE; i++) thee->bmap[0][i] = i;
1219 thee->bisectEdge = Vfetk_PDE_bisectEdge;
1220 thee->mapBoundary = Vfetk_PDE_mapBoundary;
1221 thee->markSimplex = Vfetk_PDE_markSimplex;
1222 thee->oneChart = Vfetk_PDE_oneChart;
1233 if ((*thee) != VNULL) {
1250VPRIVATE
double smooth(
int nverts,
double dist[
VAPBS_NVS],
double coeff[
VAPBS_NVS],
int meth) {
1257 for (i=0; i<nverts; i++) {
1258 if (dist[i] < VSMALL)
return coeff[i];
1259 weight = 1.0/dist[i];
1261 num += (weight * coeff[i]);
1263 }
else if (meth == 1) {
1266 if (coeff[i] < VSMALL) {
1270 num += weight; den += (weight/coeff[i]);
1279VPRIVATE
double diel() {
1282 double eps, epsp, epsw, dist[5], coeff[5], srad, swin, *vx;
1291 srfm = pbeparm->
srfm;
1292 srad = pbeparm->
srad;
1293 swin = pbeparm->
swin;
1298 if (VABS(epsp - epsw) < VSMALL)
return epsp;
1304 for (i=0; i<var.
nverts; i++) {
1307 for (j=0; j<3; j++) {
1308 dist[i] += VSQR(var.
xq[j] - vx[j]);
1310 dist[i] = VSQRT(dist[i]);
1311 coeff[i] = (epsw-epsp)*
Vacc_molAcc(acc, var.
xq, srad) + epsp;
1313 eps = smooth(var.
nverts, dist, coeff, 1);
1319 Vnm_print(2,
"Undefined surface method (%d)!\n", srfm);
1326VPRIVATE
double ionacc() {
1329 double dist[5], coeff[5], irad, swin, *vx, accval;
1336 srfm = pbeparm->
srfm;
1338 swin = pbeparm->
swin;
1341 if (var.
zks2 < VSMALL)
return 0.0;
1344 accval = Vacc_ivdwAcc(acc, var.
xq, irad);
1347 for (i=0; i<var.
nverts; i++) {
1350 for (j=0; j<3; j++) {
1351 dist[i] += VSQR(var.
xq[j] - vx[j]);
1353 dist[i] = VSQRT(dist[i]);
1354 coeff[i] = Vacc_ivdwAcc(acc, var.
xq, irad);
1356 accval = smooth(var.
nverts, dist, coeff, 1);
1362 Vnm_print(2,
"Undefined surface method (%d)!\n", srfm);
1369VPRIVATE
double debye_U(
Vpbe *pbe,
int d,
double x[]) {
1371 double size, *position, charge, xkappa, eps_w, dist, T, pot, val;
1389 for (i=0; i<d; i++) {
1390 dist += VSQR(position[i] - x[i]);
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));
1404VPRIVATE
double debye_Udiff(
Vpbe *pbe,
int d,
double x[]) {
1406 double size, *position, charge, eps_p, dist, T, pot, val;
1412 Ufull = debye_U(pbe, d, x);
1426 for (i=0; i<d; i++) {
1427 dist += VSQR(position[i] - x[i]);
1429 dist = (1.0e-10)*VSQRT(dist);
1430 val = (charge)/(4*VPI*
Vunit_eps0*eps_p*dist);
1440VPRIVATE
void coulomb(
Vpbe *pbe,
int d,
double pt[],
double eps,
double *U,
1441 double dU[],
double *d2U) {
1444 double T, pot, fx, fy, fz, x, y, z, scale;
1445 double *position, charge, dist, dist2, val, vec[3], dUold[3], Uold;
1452 pot = 0; fx = 0; fy = 0; fz = 0;
1453 x = pt[0]; y = pt[1]; z = pt[2];
1457 Vnm_print(2,
"Error calculating Green's function!\n");
1473 Uold = 0.0; dUold[0] = 0.0; dUold[1] = 0.0; dUold[2] = 0.0;
1479 for (i=0; i<d; i++) {
1480 vec[i] = (position[i] - pt[i]);
1481 dist2 += VSQR(vec[i]);
1483 dist = VSQRT(dist2);
1486 Uold = Uold + charge/dist;
1489 for (i=0; i<d; i++) dUold[i] = dUold[i] + vec[i]*charge/(dist2*dist);
1493 for (i=0; i<d; i++) {
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]);
1528VPUBLIC
void Vfetk_PDE_initElement(PDE *thee,
int elementType,
int chart,
1529 double tvx[][3],
void *data) {
1536 VASSERT(data != NULL);
1537 var.
simp = (SS *)data;
1540 var.
sType = elementType;
1543 var.
nverts = thee->dim+1;
1544 for (i=0; i<thee->dim+1; i++) var.
verts[i] = SS_vertex(var.
simp, i);
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];
1565 for (i=0; i<thee->dim; i++) var.
nvec[i] = tnvec[i];
1568 var.
fType = faceType;
1571VPUBLIC
void Vfetk_PDE_initPoint(PDE *thee,
int pointType,
int chart,
1572 double txq[],
double tU[],
double tdU[][3]) {
1575 double u2, coef2, eps_p;
1585 if ((pdetype ==
PBE_LRPBE) || (pdetype == PBE_NRPBE)) {
1586 coulomb(pbe, thee->dim, txq, eps_p, &(var.
W), var.
dW, &(var.
d2W));
1588 for (i=0; i<thee->vec; i++) {
1590 for (j=0; j<thee->dim; j++) {
1592 var.
dU[i][j] = tdU[i][j];
1597 if (pointType == 0) {
1601 var.ionacc = ionacc();
1603 var.
F = (var.diel - eps_p);
1609 var.
B = var.
DB*var.
U[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];
1621 coef2 = -1.0 * var.ionacc * var.
zks2 * var.
ionConc[i];
1622 var.
B += (coef2 *
Vcap_exp(u2, &ichop));
1624 coef2 = -1.0 * var.
ionQ[i] * coef2;
1632 var.
B = var.
DB*(var.
U[0]+var.
W);
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];
1644 coef2 = -1.0 * var.ionacc * var.
zks2 * var.
ionConc[i];
1645 var.
B += (coef2 *
Vcap_exp(u2, &ichop));
1648 coef2 = -1.0 * var.
ionQ[i] * coef2;
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];
1663 coef2 = -1.0 * var.ionacc * var.
zks2 * var.
ionConc[i];
1664 var.
B += (coef2 *
Vcap_exp(u2, &ichop));
1666 coef2 = -1.0 * var.
ionQ[i] * coef2;
1672 Vnm_print(2,
"Vfetk_PDE_initPoint: Unknown PBE type (%d)!\n",
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");
1719 for (i=0; i<thee->dim; i++) value += ( var.
A * var.
dU[0][i] * dV[0][i] );
1720 value += var.
B * V[0];
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]);
1733 Vnm_print(2,
"Vfetk: Whoa! I was just asked to evaluate a boundary weak form for point type %d!\n", key);
1741VPUBLIC
double Vfetk_PDE_DFu_wv(
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] );
1766 Vnm_print(2,
"Vfetk: Whoa! I was just asked to evaluate a boundary weak form for point type %d!\n", key);
1776#define VRINGMAX 1000
1779#define VATOMMAX 1000000
1781 void *user,
double F[]) {
1783 int iatom, jatom, natoms, atomIndex, atomList[
VATOMMAX], nAtomList;
1784 int gotAtom, numSring, isimp, ivert, sid;
1789 VV *vertex = (VV *)user;
1796 VASSERT( vertex != VNULL);
1798 sring[numSring] = VV_firstSS(vertex);
1799 while (sring[numSring] != VNULL) {
1801 sring[numSring] = SS_link(sring[numSring-1], vertex);
1803 VASSERT( numSring > 0 );
1810 for (isimp=0; isimp<numSring; isimp++) {
1811 sid = SS_id(sring[isimp]);
1813 for (iatom=0; iatom<natoms; iatom++) {
1818 for (jatom=0; jatom<nAtomList; jatom++) {
1819 if (atomList[jatom] == atomIndex) {
1826 atomList[nAtomList] = atomIndex;
1838 sring[isimp], position, phi, phix)) {
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: \
1845 for (ivert=0; ivert<Gem_dimVV(
Vfetk_getGem(var.
fetk)); ivert++) Vnm_print(2,
"%e ", phi[ivert]);
1846 Vnm_print(2,
"}\n");
1851 if (VV_id(SS_vertex(sring[isimp], ivert)) == VV_id(vertex)) value += phi[ivert];
1859 }
else if ((pdetype ==
PBE_LRPBE) || (pdetype == PBE_NRPBE)) {
1861 }
else { VASSERT(0); }
1871 F[0] = debye_U(var.
fetk->
pbe, thee->dim, txq);
1873 F[0] = debye_Udiff(var.
fetk->
pbe, thee->dim, txq);
1902VPUBLIC
void Vfetk_PDE_bisectEdge(
int dim,
int dimII,
int edgeType,
1903 int chart[],
double vx[][3]) {
1907 for (i=0; i<dimII; i++) vx[2][i] = .5 * (vx[0][i] + vx[1][i]);
1908 chart[2] = chart[0];
1912VPUBLIC
void Vfetk_PDE_mapBoundary(
int dim,
int dimII,
int vertexType,
1913 int chart,
double vx[3]) {
1917VPUBLIC
int Vfetk_PDE_markSimplex(
int dim,
int dimII,
int simplexType,
1921 double targetRes, edgeLength, srad, swin, myAcc, refAcc;
1940 srfm = pbeparm->
srfm;
1941 srad = pbeparm->
srad;
1942 swin = pbeparm->
swin;
1943 simp = (SS *)simplex;
1948 Gem_longestEdge(var.
fetk->
gm, simp, -1, &edgeLength);
1949 if (edgeLength < targetRes)
return 0;
1965 for (i=1; i<(dim+1); i++) {
1967 if (myAcc != refAcc) {
1974 for (i=1; i<(dim+1); i++) {
1976 if (myAcc != refAcc) {
1983 for (i=1; i<(dim+1); i++) {
1985 if (myAcc != refAcc) {
1998VPUBLIC
void Vfetk_PDE_oneChart(
int dim,
int dimII,
int objType,
int chart[],
1999 double vx[][3],
int dimV) {
2006 double dielE, qmE, coef2, u2;
2015 for (i=0; i<3; i++) dielE += VSQR(var.
dU[0][i]);
2016 dielE = dielE*var.diel;
2020 if ((var.ionacc > VSMALL) && (var.
zkappa2 > VSMALL)) {
2021 qmE = var.ionacc*var.
zkappa2*VSQR(var.
U[0]);
2025 if ((var.ionacc > VSMALL) && (var.
zks2 > VSMALL)) {
2027 for (i=0; i<var.
nion; i++) {
2029 u2 = -1.0 * (var.
U[0]) * var.
ionQ[i];
2030 qmE += (coef2 * (
Vcap_exp(u2, &ichop) - 1.0));
2035 if ((var.ionacc > VSMALL) && (var.
zkappa2 > VSMALL)) {
2036 qmE = var.ionacc*var.
zkappa2*VSQR((var.
U[0] + var.
W));
2040 if ((var.ionacc > VSMALL) && (var.
zks2 > VSMALL)) {
2042 for (i=0; i<var.
nion; i++) {
2044 u2 = -1.0 * (var.
U[0] + var.
W) * var.
ionQ[i];
2045 qmE += (coef2 * (
Vcap_exp(u2, &ichop) - 1.0));
2050 if ((var.ionacc > VSMALL) && (var.
zks2 > VSMALL)) {
2052 for (i=0; i<var.
nion; i++) {
2054 u2 = -1.0 * (var.
U[0]) * var.
ionQ[i];
2055 qmE += (coef2 * (
Vcap_exp(u2, &ichop) - 1.0));
2060 Vnm_print(2,
"Vfetk_PDE_Ju: Invalid PBE type (%d)!\n", type);
2068 }
else if (key == 1) {
2083 VASSERT(var.
fetk != VNULL);
2085 VASSERT(csm != VNULL);
2090 Vnm_print(2,
"Error while updating charge-simplex map!\n");
2095VPRIVATE
void polyEval(
int numP,
double p[],
double c[][VMAXP],
double xv[]) {
2102 for (i=0; i<numP; i++) {
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]) {
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];
2149 if ((key == 0) || (key == 1)) {
2151 }
else if ((key == 2) || (key == 3)) {
2153 }
else { VASSERT(0); }
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) {
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) {
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) {
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..");
2204 double xq[],
double basis[]) {
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); }
2217VPRIVATE
void init_2DP1(
int dimIS[],
int *ndof,
int dof[],
double c[][VMAXP],
2218 double cx[][VMAXP],
double cy[][VMAXP],
double cz[][VMAXP]) {
2228 for (i=0; i<
VAPBS_NVS; i++) *ndof += dimIS[i] * dof[i];
2229 VASSERT( *ndof == dim_2DP1 );
2230 VASSERT( *ndof <= VMAXP );
2233 setCoef( *ndof, c, cx, cy, cz, lgr_2DP1, lgr_2DP1x, lgr_2DP1y, lgr_2DP1z );
2236VPRIVATE
void init_3DP1(
int dimIS[],
int *ndof,
int dof[],
double c[][VMAXP],
2237 double cx[][VMAXP],
double cy[][VMAXP],
double cz[][VMAXP]) {
2247 for (i=0; i<
VAPBS_NVS; i++) *ndof += dimIS[i] * dof[i];
2248 VASSERT( *ndof == dim_3DP1 );
2249 VASSERT( *ndof <= VMAXP );
2252 setCoef( *ndof, c, cx, cy, cz, lgr_3DP1, lgr_3DP1x, lgr_3DP1y, lgr_3DP1z );
2259 Vnm_print(1,
"DEBUG: nvec = (%g, %g, %g)\n", var.
nvec[0], var.
nvec[1],
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]);
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],
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]);
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);
2302 double coord[3], chi, q, conc, val;
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),
2328 Vnm_print(2,
"Vfetk_fillArray: can't write out charge distribution!\n");
2338 Bvec_axpy(vec, u_d, 1.0);
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);
2346 Bvec_set(vec, i, chi);
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);
2355 Bvec_set(vec, i, chi);
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);
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);
2373 Bvec_set(vec, i, chi);
2378 Vnm_print(2,
"Vfetk_fillArray: can't write out Laplacian!\n");
2383 Vnm_print(2,
"Vfetk_fillArray: can't write out energy density!\n");
2393 Bvec_axpy(vec, u_d, 1.0);
2396 for (i=0; i<Gem_numVV(gm); i++) {
2398 for (j=0; j<pbe->
numIon; j++) {
2402 val += (conc*
Vcap_exp(-q*Bvec_val(vec, i), &ichop));
2404 val += (conc * ( 1 - q*Bvec_val(vec, i)));
2407 Bvec_set(vec, i, val);
2417 Bvec_axpy(vec, u_d, 1.0);
2420 for (i=0; i<Gem_numVV(gm); i++) {
2422 for (j=0; j<pbe->
numIon; j++) {
2426 val += (q*conc*
Vcap_exp(-q*Bvec_val(vec, i), &ichop));
2428 val += (q*conc*(1 - q*Bvec_val(vec, i)));
2431 Bvec_set(vec, i, val);
2436 Vnm_print(2,
"Vfetk_fillArray: can't write out x-shifted diel!\n");
2441 Vnm_print(2,
"Vfetk_fillArray: can't write out y-shifted diel!\n");
2446 Vnm_print(2,
"Vfetk_fillArray: can't write out z-shifted diel!\n");
2451 Vnm_print(2,
"Vfetk_fillArray: can't write out kappa!\n");
2456 Vnm_print(2,
"Vfetk_fillArray: invalid data type (%d)!\n", type);
2465 const char *thost,
const char *fname, Bvec *vec,
Vdata_Format format) {
2472 VASSERT(thee != VNULL);
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",
2482 if (Vio_connect(sock, 0) < 0) {
2483 Vnm_print(2,
"Vfetk_write: Problem connecting to virtual socket %s\n",
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),
2499 Aprx_writeSOL(aprx, sock, vec,
"DX");
2502 Aprx_writeSOL(aprx, sock, vec,
"UCD");
2505 Vnm_print(2,
"Vfetk_write: UHBD format not supported!\n");
2508 Vnm_print(2,
"Vfetk_write: Invalid data format (%d)!\n", format);
2513 Vio_connectFree(sock);
VPUBLIC double Vacc_molAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility.
VPUBLIC double Vacc_splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad)
Report spline-based accessibility.
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
VPUBLIC double Vatom_getPartID(Vatom *thee)
Get partition ID.
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
VPUBLIC void Vatom_setPartID(Vatom *thee, int partID)
Set partition ID.
VPUBLIC double Vcap_exp(double x, int *ichop)
Provide a capped exp() function.
VPUBLIC int Vcsm_getNumberAtoms(Vcsm *thee, int isimp)
Get number of atoms associated with a simplex.
VPUBLIC int Vcsm_getAtomIndex(Vcsm *thee, int iatom, int isimp)
Get ID of particular atom in a simplex.
VPUBLIC Vcsm * Vcsm_ctor(Valist *alist, Gem *gm)
Construct Vcsm object.
VPUBLIC SS * Vcsm_getSimplex(Vcsm *thee, int isimp, int iatom)
Get particular simplex associated with an atom.
VPUBLIC int Vcsm_update(Vcsm *thee, SS **simps, int num)
Update the charge-simplex and simplex-charge maps after refinement.
VPUBLIC Vatom * Vcsm_getAtom(Vcsm *thee, int iatom, int isimp)
Get particular atom associated with a simplex.
VPUBLIC void Vcsm_init(Vcsm *thee)
Initialize charge-simplex map with mesh and atom data.
VPUBLIC int Vcsm_getNumberSimplices(Vcsm *thee, int iatom)
Get number of simplices associated with an atom.
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.
VPUBLIC void Vcsm_dtor(Vcsm **thee)
Destroy Vcsm object.
enum eVfetk_MeshLoad Vfetk_MeshLoad
Declare FEMparm_GuessType type.
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.
VPUBLIC double Vfetk_PDE_Ju(PDE *thee, int key)
Energy functional. This returns the energy (less delta function terms) in the form:
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).
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...
VPUBLIC void Vfetk_PDE_initFace(PDE *thee, int faceType, int chart, double tnvec[])
Do once-per-face initialization.
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,...
VPUBLIC void Vfetk_dtor(Vfetk **thee)
Object destructor.
VPUBLIC void Vfetk_PDE_dtor(PDE **thee)
Destroys FEtk PDE object.
VPUBLIC int Vfetk_getAtomColor(Vfetk *thee, int iatom)
Get the partition information for a particular atom.
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.
VPUBLIC void Bmat_printHB(Bmat *thee, char *fname)
Writes a Bmat to disk in Harwell-Boeing sparse matrix format.
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:
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,...
VPUBLIC Gem * Vfetk_getGem(Vfetk *thee)
Get a pointer to the Gem (grid manager) object.
struct sVfetk Vfetk
Declaration of the Vfetk class as the Vfetk structure.
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.
VPUBLIC Vfetk * Vfetk_ctor(Vpbe *pbe, Vhal_PBEType type)
Constructor for Vfetk object.
VPUBLIC double Vfetk_dqmEnergy(Vfetk *thee, int color)
Get the "mobile charge" and "polarization" contributions to the electrostatic energy.
#define VATOMMAX
Maximum number of atoms associated with a vertex.
VPUBLIC unsigned long int Vfetk_memChk(Vfetk *thee)
Return the memory used by this structure (and its contents) in bytes.
#define VRINGMAX
Maximum number of simplices in a simplex ring.
VPUBLIC double * Vfetk_getSolution(Vfetk *thee, int *length)
Create an array containing the solution (electrostatic potential in units of ) at the finest mesh lev...
VPUBLIC void Vfetk_PDE_initAssemble(PDE *thee, int ip[], double rp[])
Do once-per-assembly initialization.
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)
VPUBLIC PDE * Vfetk_PDE_ctor(Vfetk *fetk)
Constructs the FEtk PDE object.
VPUBLIC void Vfetk_PDE_dtor2(PDE *thee)
FORTRAN stub: destroys FEtk PDE object.
VPUBLIC int Vfetk_PDE_ctor2(PDE *thee, Vfetk *fetk)
Intializes the FEtk PDE object.
VPUBLIC void Vfetk_dtor2(Vfetk *thee)
FORTRAN stub object destructor.
VPUBLIC double Vfetk_qfEnergy(Vfetk *thee, int color)
Get the "fixed charge" contribution to the electrostatic energy.
VPUBLIC Vcsm * Vfetk_getVcsm(Vfetk *thee)
Get a pointer to the Vcsm (charge-simplex map) object.
VPUBLIC double Vfetk_energy(Vfetk *thee, int color, int nonlin)
Return the total electrostatic energy.
VPUBLIC Vpbe * Vfetk_getVpbe(Vfetk *thee)
Get a pointer to the Vpbe (PBE manager) object.
VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type)
Fill an array with the specified data.
VPUBLIC void Vfetk_PDE_Fu(PDE *thee, int key, double F[])
Evaluate strong form of PBE. For interior points, this is:
VPUBLIC int Vfetk_ctor2(Vfetk *thee, Vpbe *pbe, Vhal_PBEType type)
FORTRAN stub constructor for Vfetk object.
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.
VPUBLIC void Vfetk_setAtomColors(Vfetk *thee)
Transfer color (partition ID) information frmo a partitioned mesh to the atoms.
VPUBLIC AM * Vfetk_getAM(Vfetk *thee)
Get a pointer to the AM (algebra manager) object.
VPUBLIC void Vfetk_dumpLocalVar()
Debugging routine to print out local variables used by PDE object.
VPUBLIC void Vfetk_setParameters(Vfetk *thee, PBEparm *pbeparm, FEMparm *feparm)
Set the parameter objects.
VPUBLIC void Vgreen_dtor(Vgreen **thee)
Destruct the Green's function oracle.
VPUBLIC Vgreen * Vgreen_ctor(Valist *alist)
Construct the Green's function oracle.
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...
enum eVdata_Format Vdata_Format
Declaration of the Vdata_Format type as the Vdata_Format enum.
#define VAPBS_NVS
Number of vertices per simplex (hard-coded to 3D)
enum eVsurf_Meth Vsurf_Meth
Declaration of the Vsurf_Meth type as the Vsurf_Meth enum.
enum eVhal_PBEType Vhal_PBEType
Declaration of the Vhal_PBEType type as the Vhal_PBEType enum.
enum eVdata_Type Vdata_Type
Declaration of the Vdata_Type type as the Vdata_Type enum.
#define VAPBS_DIM
Our dimension.
VPUBLIC double Vpbe_getSoluteDiel(Vpbe *thee)
Get solute dielectric constant.
VPUBLIC double Vpbe_getZkappa2(Vpbe *thee)
Get modified squared Debye-Huckel parameter.
VPUBLIC double Vpbe_getXkappa(Vpbe *thee)
Get Debye-Huckel parameter.
VPUBLIC double Vpbe_getZmagic(Vpbe *thee)
Get charge scaling factor.
VPUBLIC double Vpbe_getSolventDiel(Vpbe *thee)
Get solvent dielectric constant.
VPUBLIC double Vpbe_getBulkIonicStrength(Vpbe *thee)
Get bulk ionic strength.
VPUBLIC double Vpbe_getMaxIonRadius(Vpbe *thee)
Get maximum radius of ion species.
VPUBLIC Valist * Vpbe_getValist(Vpbe *thee)
Get atom list.
VPUBLIC int Vpbe_getIons(Vpbe *thee, int *nion, double ionConc[MAXION], double ionRadii[MAXION], double ionQ[MAXION])
Get information about the counterion species present.
VPUBLIC double Vpbe_getTemperature(Vpbe *thee)
Get temperature.
#define Vunit_ec
Charge of an electron in C.
#define Vunit_kb
Boltzmann constant.
#define Vunit_eps0
Vacuum permittivity.
Parameter structure for FEM-specific variables from input files.
Parameter structure for PBE variables from input files.
Oracle for solvent- and ion-accessibility around a biomolecule.
Container class for list of atom objects.
Contains public data members for Vatom class/module.
Charge-simplex map class.
double dU[MAXV][VAPBS_DIM]
Contains public data members for Vfetk class/module.
Contains public data members for Vpbe class/module.
Contains declarations for class Vfetk.