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;
599 if (var.ionstr > 0.0) var.zks2 = 0.5*var.zkappa2/var.ionstr;
601 Vpbe_getIons(var.fetk->pbe, &(var.nion), var.ionConc, var.ionRadii,
603 for (i=0; i<var.nion; i++) {
604 var.ionConc[i] = var.zks2 * var.ionConc[i] * var.ionQ[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;
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;
1298 if (VABS(epsp - epsw) < VSMALL)
return epsp;
1301 eps = ((epsw-epsp)*
Vacc_molAcc(acc, var.xq, srad) + 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);
1316 eps = ((epsw-epsp)*
Vacc_splineAcc(acc, var.xq, swin, 0.0) + epsp);
1319 Vnm_print(2,
"Undefined surface method (%d)!\n", srfm);
1326VPRIVATE
double ionacc() {
1329 double dist[5], coeff[5], irad, swin, *vx, accval;
1334 VASSERT(var.fetk->pbeparm != VNULL);
1335 pbeparm = var.fetk->pbeparm;
1336 srfm = pbeparm->
srfm;
1338 swin = pbeparm->
swin;
1339 acc = var.fetk->pbe->
acc;
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];
1456 if (!
Vgreen_coulombD(var.green, 1, &x, &y, &z, &pot, &fx, &fy, &fz)) {
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]);
1513 if (var.initGreen) {
1520 if (!var.initGreen) {
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;
1580 pdetype = var.fetk->type;
1581 pbe = var.fetk->pbe;
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);
1608 var.DB = var.ionacc*var.zkappa2*var.ionstr;
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;
1625 var.DB += (coef2 *
Vcap_exp(u2, &ichop));
1631 var.DB = var.ionacc*var.zkappa2*var.ionstr;
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;
1649 var.DB += (coef2 *
Vcap_exp(u2, &ichop));
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;
1667 var.DB += (coef2 *
Vcap_exp(u2, &ichop));
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");
1714 type = var.fetk->type;
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(
1754 type = var.fetk->type;
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;
1791 pdetype = var.fetk->type;
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");
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];
1859 }
else if ((pdetype ==
PBE_LRPBE) || (pdetype == PBE_NRPBE)) {
1861 }
else { VASSERT(0); }
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);
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;
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;
1940 srfm = pbeparm->
srfm;
1941 srad = pbeparm->
srad;
1942 swin = pbeparm->
swin;
1943 simp = (SS *)simplex;
1944 type = var.fetk->type;
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;
2010 type = var.fetk->type;
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++) {
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));
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++) {
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));
2050 if ((var.ionacc > VSMALL) && (var.zks2 > VSMALL)) {
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));
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);
struct sFEMparm FEMparm
Declaration of the FEMparm class as the FEMparm structure.
struct sPBEparm PBEparm
Declaration of the PBEparm class as the PBEparm structure.
struct sVacc Vacc
Declaration of the Vacc class as the Vacc structure.
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.
struct sValist Valist
Declaration of the Valist class as the Valist structure.
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.
struct sVatom Vatom
Declaration of the Vatom class as the Vatom structure.
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.
struct sVcsm Vcsm
Declaration of the Vcsm class as the Vcsm structure.
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.
struct sVfetk_LocalVar Vfetk_LocalVar
Declaration of the Vfetk_LocalVar subclass as the Vfetk_LocalVar structure.
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.
struct sVpbe Vpbe
Declaration of the Vpbe class as the Vpbe structure.
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.
Contains declarations for class Vfetk.