68 Vnm_tprint( 1,
"Reading parameter data from %s.\n",
72 Vnm_tprint(2,
"Error reading parameter file (%s)!\n", nosh->
parmpath);
77 Vnm_tprint( 1,
"Reading parameter data from %s.\n",
81 Vnm_tprint(2,
"Error reading parameter file (%s)!\n", nosh->
parmpath);
86 Vnm_tprint(2,
"Error! Undefined parameter file type (%d)!\n", nosh->
parmfmt);
103 Vnm_tprint( 1,
"Got paths for %d molecules\n", nosh->
nmol);
104 if (nosh->
nmol <= 0) {
105 Vnm_tprint(2,
"You didn't specify any molecules (correctly)!\n");
106 Vnm_tprint(2,
"Bailing out!\n");
111 if (param == VNULL) {
112 Vnm_tprint(2,
"Error! You don't have a valid parameter object!\n");
118 for (i=0; i<nosh->
nmol; i++) {
119 if(alist[i] == VNULL){
126 switch (nosh->
molfmt[i]) {
131 Vnm_print(2,
"\nWARNING!! Radius/charge information from PQR file %s\n", nosh->
molpath[i]);
132 Vnm_print(2,
"will be replaced with data from parameter file (%s)!\n", nosh->
parmpath);
134 Vnm_tprint( 1,
"Reading PQR-format atom data from %s.\n",
136 sock = Vio_ctor(
"FILE",
"ASC", VNULL, nosh->
molpath[i],
"r");
138 Vnm_print(2,
"Problem opening virtual socket %s!\n",
142 if (Vio_accept(sock, 0) < 0) {
143 Vnm_print(2,
"Problem accepting virtual socket %s!\n",
152 if(rc == 0)
return 0;
154 Vio_acceptFree(sock);
160 Vnm_tprint(2,
"NOsh: Error! Can't read PDB without specifying PARM file!\n");
163 Vnm_tprint( 1,
"Reading PDB-format atom data from %s.\n",
165 sock = Vio_ctor(
"FILE",
"ASC", VNULL, nosh->
molpath[i],
"r");
167 Vnm_print(2,
"Problem opening virtual socket %s!\n",
171 if (Vio_accept(sock, 0) < 0) {
172 Vnm_print(2,
"Problem accepting virtual socket %s!\n", nosh->
molpath[i]);
181 Vio_acceptFree(sock);
185 Vnm_tprint( 1,
"Reading XML-format atom data from %s.\n",
187 sock = Vio_ctor(
"FILE",
"ASC", VNULL, nosh->
molpath[i],
"r");
189 Vnm_print(2,
"Problem opening virtual socket %s!\n",
193 if (Vio_accept(sock, 0) < 0) {
194 Vnm_print(2,
"Problem accepting virtual socket %s!\n",
206 Vio_acceptFree(sock);
210 Vnm_tprint(2,
"NOsh: Error! Undefined molecule file type \ 211 (%d)!\n", nosh->
molfmt[i]);
216 Vnm_tprint( 2,
"Error while reading molecule from %s\n",
222 Vnm_tprint( 1,
" Centered at (%4.3e, %4.3e, %4.3e)\n",
223 alist[i]->center[0], alist[i]->center[1],
224 alist[i]->center[2]);
225 Vnm_tprint( 1,
" Net charge %3.2e e\n", alist[i]->charge);
238 Vnm_tprint( 1,
"Destroying %d molecules\n", nosh->
nmol);
241 for (i=0; i<nosh->
nmol; i++)
252 Vgrid *dielYMap[NOSH_MAXMOL],
253 Vgrid *dielZMap[NOSH_MAXMOL]
271 Vnm_tprint( 1,
"Got paths for %d dielectric map sets\n",
277 for (i=0; i<nosh->
ndiel; i++) {
278 Vnm_tprint( 1,
"Reading x-shifted dielectric map data from \ 280 dielXMap[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
288 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
294 nx = dielXMap[i]->nx;
295 ny = dielXMap[i]->ny;
296 nz = dielXMap[i]->nz;
299 hx = dielXMap[i]->hx;
300 hy = dielXMap[i]->hy;
301 hzed = dielXMap[i]->hzed;
304 xmin = dielXMap[i]->xmin;
305 ymin = dielXMap[i]->ymin;
306 zmin = dielXMap[i]->zmin;
307 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
308 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
309 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
312 for (ii=0; ii<(nx*ny*
nz); ii++)
313 sum += (dielXMap[i]->data[ii]);
314 sum = sum*hx*hy*
hzed;
315 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
320 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
326 nx = dielXMap[i]->nx;
327 ny = dielXMap[i]->ny;
328 nz = dielXMap[i]->nz;
331 hx = dielXMap[i]->hx;
332 hy = dielXMap[i]->hy;
333 hzed = dielXMap[i]->hzed;
336 xmin = dielXMap[i]->xmin;
337 ymin = dielXMap[i]->ymin;
338 zmin = dielXMap[i]->zmin;
339 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
340 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
341 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
344 for (ii=0; ii<(nx*ny*
nz); ii++)
345 sum += (dielXMap[i]->data[ii]);
346 sum = sum*hx*hy*
hzed;
347 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
351 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
355 Vnm_tprint( 2,
"AVS input not supported yet!\n");
359 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
362 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
366 Vnm_tprint( 1,
"Reading y-shifted dielectric map data from \ 368 dielYMap[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
376 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
382 nx = dielYMap[i]->nx;
383 ny = dielYMap[i]->ny;
384 nz = dielYMap[i]->nz;
387 hx = dielYMap[i]->hx;
388 hy = dielYMap[i]->hy;
389 hzed = dielYMap[i]->hzed;
392 xmin = dielYMap[i]->xmin;
393 ymin = dielYMap[i]->ymin;
394 zmin = dielYMap[i]->zmin;
395 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
396 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
397 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
400 for (ii=0; ii<(nx*ny*
nz); ii++)
401 sum += (dielYMap[i]->data[ii]);
402 sum = sum*hx*hy*
hzed;
403 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
408 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
414 nx = dielYMap[i]->nx;
415 ny = dielYMap[i]->ny;
416 nz = dielYMap[i]->nz;
419 hx = dielYMap[i]->hx;
420 hy = dielYMap[i]->hy;
421 hzed = dielYMap[i]->hzed;
424 xmin = dielYMap[i]->xmin;
425 ymin = dielYMap[i]->ymin;
426 zmin = dielYMap[i]->zmin;
427 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
428 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
429 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
432 for (ii=0; ii<(nx*ny*
nz); ii++)
433 sum += (dielYMap[i]->data[ii]);
434 sum = sum*hx*hy*
hzed;
435 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
439 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
443 Vnm_tprint( 2,
"AVS input not supported yet!\n");
447 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
450 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
455 Vnm_tprint( 1,
"Reading z-shifted dielectric map data from \ 457 dielZMap[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
465 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
471 nx = dielZMap[i]->nx;
472 ny = dielZMap[i]->ny;
473 nz = dielZMap[i]->nz;
476 hx = dielZMap[i]->hx;
477 hy = dielZMap[i]->hy;
478 hzed = dielZMap[i]->hzed;
481 xmin = dielZMap[i]->xmin;
482 ymin = dielZMap[i]->ymin;
483 zmin = dielZMap[i]->zmin;
484 Vnm_tprint(1,
" %d x %d x %d grid\n",
486 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
488 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
491 for (ii=0; ii<(nx*ny*
nz); ii++) sum += (dielZMap[i]->data[ii]);
492 sum = sum*hx*hy*
hzed;
493 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
498 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
504 nx = dielZMap[i]->nx;
505 ny = dielZMap[i]->ny;
506 nz = dielZMap[i]->nz;
509 hx = dielZMap[i]->hx;
510 hy = dielZMap[i]->hy;
511 hzed = dielZMap[i]->hzed;
514 xmin = dielZMap[i]->xmin;
515 ymin = dielZMap[i]->ymin;
516 zmin = dielZMap[i]->zmin;
517 Vnm_tprint(1,
" %d x %d x %d grid\n",
519 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
521 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
524 for (ii=0; ii<(nx*ny*
nz); ii++) sum += (dielZMap[i]->data[ii]);
525 sum = sum*hx*hy*
hzed;
526 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
530 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
534 Vnm_tprint( 2,
"AVS input not supported yet!\n");
538 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
541 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
552 Vgrid *dielYMap[NOSH_MAXMOL],
553 Vgrid *dielZMap[NOSH_MAXMOL]) {
557 if (nosh->
ndiel > 0) {
559 Vnm_tprint( 1,
"Destroying %d dielectric map sets\n",
562 for (i=0; i<nosh->
ndiel; i++) {
585 Vnm_tprint( 1,
"Got paths for %d kappa maps\n", nosh->
nkappa);
588 for (i=0; i<nosh->
nkappa; i++) {
589 Vnm_tprint( 1,
"Reading kappa map data from %s:\n",
591 map[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
599 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
603 Vnm_tprint(1,
" %d x %d x %d grid\n",
604 map[i]->
nx, map[i]->
ny, map[i]->
nz);
605 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
606 map[i]->
hx, map[i]->
hy, map[i]->
hzed);
607 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
610 for (ii = 0, len = map[i]->
nx * map[i]->
ny * map[i]->
nz;
614 sum += (map[i]->data[ii]);
616 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
617 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
621 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
625 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
629 Vnm_tprint( 2,
"AVS input not supported yet!\n");
634 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
638 Vnm_tprint(1,
" %d x %d x %d grid\n",
639 map[i]->
nx, map[i]->
ny, map[i]->
nz);
640 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
641 map[i]->
hx, map[i]->
hy, map[i]->
hzed);
642 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
645 for (ii=0, len=map[i]->
nx*map[i]->
ny*map[i]->
nz; ii<len; ii++) {
646 sum += (map[i]->data[ii]);
648 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
649 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
652 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
668 Vnm_tprint( 1,
"Destroying %d kappa maps\n", nosh->
nkappa);
689 Vnm_tprint( 1,
"Got paths for %d potential maps\n", nosh->
npot);
692 for (i=0; i<nosh->
npot; i++) {
693 Vnm_tprint( 1,
"Reading potential map data from %s:\n",
695 map[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
696 switch (nosh->
potfmt[i]) {
704 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
710 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
715 Vnm_tprint(1,
" %d x %d x %d grid\n",
716 map[i]->
nx, map[i]->
ny, map[i]->
nz);
717 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
718 map[i]->
hx, map[i]->
hy, map[i]->
hzed);
719 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
722 for (ii=0,len=map[i]->
nx*map[i]->
ny*map[i]->
nz; ii<len; ii++) {
723 sum += (map[i]->data[ii]);
725 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
726 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
730 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
734 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
738 Vnm_tprint( 2,
"AVS input not supported yet!\n");
741 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
757 if (nosh->
npot > 0) {
759 Vnm_tprint( 1,
"Destroying %d potential maps\n", nosh->
npot);
780 Vnm_tprint( 1,
"Got paths for %d charge maps\n", nosh->
ncharge);
783 for (i=0; i<nosh->
ncharge; i++) {
784 Vnm_tprint( 1,
"Reading charge map data from %s:\n",
786 map[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
793 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
797 Vnm_tprint(1,
" %d x %d x %d grid\n",
798 map[i]->
nx, map[i]->
ny, map[i]->
nz);
799 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
800 map[i]->
hx, map[i]->
hy, map[i]->
hzed);
801 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
804 for (ii=0,len=map[i]->
nx*map[i]->
ny*map[i]->
nz; ii<len; ii++) {
805 sum += (map[i]->data[ii]);
807 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
808 Vnm_tprint(1,
" Charge map integral = %3.2e e\n", sum);
811 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
814 Vnm_tprint( 2,
"AVS input not supported yet!\n");
817 Vnm_tprint(2,
"MCSF input not supported yet!\n");
821 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
825 Vnm_tprint(1,
" %d x %d x %d grid\n",
826 map[i]->
nx, map[i]->
ny, map[i]->
nz);
827 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
828 map[i]->
hx, map[i]->
hy, map[i]->
hzed);
829 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
832 for (ii=0,len=map[i]->
nx*map[i]->
ny*map[i]->
nz; ii<len; ii++) {
833 sum += (map[i]->data[ii]);
835 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
836 Vnm_tprint(1,
" Charge map integral = %3.2e e\n", sum);
839 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
857 Vnm_tprint( 1,
"Destroying %d charge maps\n", nosh->
ncharge);
872 for (i=0; i<pbeparm->
nion; i++)
873 ionstr += 0.5*(VSQR(pbeparm->
ionq[i])*pbeparm->
ionc[i]);
875 Vnm_tprint( 1,
" Molecule ID: %d\n", pbeparm->
molid);
878 Vnm_tprint( 1,
" Nonlinear traditional PBE\n");
881 Vnm_tprint( 1,
" Linearized traditional PBE\n");
884 Vnm_tprint( 1,
" Nonlinear regularized PBE\n");
885 Vnm_tprint( 2,
" ** Sorry, but Nathan broke the nonlinear regularized PBE implementation. **\n");
886 Vnm_tprint( 2,
" ** Please let us know if you are interested in using it. **\n");
890 Vnm_tprint( 1,
" Linearized regularized PBE\n");
893 Vnm_tprint( 1,
" Nonlinear Size-Modified PBE\n");
896 Vnm_tprint(2,
" Unknown PBE type (%d)!\n", pbeparm->
pbetype);
900 Vnm_tprint( 1,
" Zero boundary conditions\n");
902 Vnm_tprint( 1,
" Single Debye-Huckel sphere boundary \ 905 Vnm_tprint( 1,
" Multiple Debye-Huckel sphere boundary \ 908 Vnm_tprint( 1,
" Boundary conditions from focusing\n");
910 Vnm_tprint( 1,
" Boundary conditions from potential map\n");
912 Vnm_tprint( 1,
" Membrane potential boundary conditions.\n");
914 Vnm_tprint( 1,
" %d ion species (%4.3f M ionic strength):\n",
915 pbeparm->
nion, ionstr);
916 for (i=0; i<pbeparm->
nion; i++) {
917 Vnm_tprint( 1,
" %4.3f A-radius, %4.3f e-charge, \ 918 %4.3f M concentration\n",
923 Vnm_tprint( 1,
" Lattice spacing: %4.3f A (SMPBE) \n", pbeparm->
smvolume);
924 Vnm_tprint( 1,
" Relative size parameter: %4.3f (SMPBE) \n", pbeparm->
smsize);
927 Vnm_tprint( 1,
" Solute dielectric: %4.3f\n", pbeparm->
pdie);
928 Vnm_tprint( 1,
" Solvent dielectric: %4.3f\n", pbeparm->
sdie);
929 switch (pbeparm->
srfm) {
931 Vnm_tprint( 1,
" Using \"molecular\" surface \ 932 definition; no smoothing\n");
933 Vnm_tprint( 1,
" Solvent probe radius: %4.3f A\n",
937 Vnm_tprint( 1,
" Using \"molecular\" surface definition;\ 938 harmonic average smoothing\n");
939 Vnm_tprint( 1,
" Solvent probe radius: %4.3f A\n",
943 Vnm_tprint( 1,
" Using spline-based surface definition;\ 944 window = %4.3f\n", pbeparm->
swin);
949 Vnm_tprint( 1,
" Temperature: %4.3f K\n", pbeparm->
temp);
951 energies will be calculated\n");
953 forces will be calculated \n");
955 solvent forces will be calculated\n");
956 for (i=0; i<pbeparm->
numwrite; i++) {
959 Vnm_tprint(1,
" Charge distribution to be written to ");
962 Vnm_tprint(1,
" Potential to be written to ");
965 Vnm_tprint(1,
" Molecular solvent accessibility \ 969 Vnm_tprint(1,
" Spline-based solvent accessibility \ 973 Vnm_tprint(1,
" van der Waals solvent accessibility \ 977 Vnm_tprint(1,
" Ion accessibility to be written to ");
980 Vnm_tprint(1,
" Potential Laplacian to be written to ");
983 Vnm_tprint(1,
" Energy density to be written to ");
986 Vnm_tprint(1,
" Ion number density to be written to ");
989 Vnm_tprint(1,
" Ion charge density to be written to ");
992 Vnm_tprint(1,
" X-shifted dielectric map to be written \ 996 Vnm_tprint(1,
" Y-shifted dielectric map to be written \ 1000 Vnm_tprint(1,
" Z-shifted dielectric map to be written \ 1004 Vnm_tprint(1,
" Kappa map to be written to ");
1007 Vnm_tprint(1,
" Atom potentials to be written to ");
1010 Vnm_tprint(2,
" Invalid data type for writing!\n");
1015 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"dx");
1018 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"dx.gz");
1021 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"grd");
1024 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"ucd");
1027 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"mcsf");
1030 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"txt");
1033 Vnm_tprint(2,
" Invalid format for writing!\n");
1043 switch (mgparm->
chgm) {
1045 Vnm_tprint(1,
" Using linear spline charge discretization.\n");
1048 Vnm_tprint(1,
" Using cubic spline charge discretization.\n");
1054 Vnm_tprint( 1,
" Partition overlap fraction = %g\n",
1056 Vnm_tprint( 1,
" Processor array = %d x %d x %d\n",
1059 Vnm_tprint( 1,
" Grid dimensions: %d x %d x %d\n",
1061 Vnm_tprint( 1,
" Grid spacings: %4.3f x %4.3f x %4.3f\n",
1063 Vnm_tprint( 1,
" Grid lengths: %4.3f x %4.3f x %4.3f\n",
1065 Vnm_tprint( 1,
" Grid center: (%4.3f, %4.3f, %4.3f)\n",
1066 realCenter[0], realCenter[1], realCenter[2]);
1067 Vnm_tprint( 1,
" Multigrid levels: %d\n", mgparm->
nlev);
1077 double realCenter[3],
1080 Vgrid *dielXMap[NOSH_MAXMOL],
1081 Vgrid *dielYMap[NOSH_MAXMOL],
1082 Vgrid *dielZMap[NOSH_MAXMOL],
1083 Vgrid *kappaMap[NOSH_MAXMOL],
1084 Vgrid *chargeMap[NOSH_MAXMOL],
1085 Vpmgp *pmgp[NOSH_MAXCALC],
1086 Vpmg *pmg[NOSH_MAXCALC],
1087 Vgrid *potMap[NOSH_MAXMOL]
1098 Vatom *atom = VNULL;
1099 Vgrid *theDielXMap = VNULL,
1100 *theDielYMap = VNULL,
1101 *theDielZMap = VNULL;
1102 Vgrid *theKappaMap = VNULL,
1104 *theChargeMap = VNULL;
1110 for (j=0; j<3; j++) realCenter[j] = mgparm->
center[j];
1114 myalist = alist[pbeparm->
molid-1];
1128 Vnm_tprint(0,
"Setting up PBE object...\n");
1130 sparm = pbeparm->
swin;
1132 sparm = pbeparm->
srad;
1134 if (pbeparm->
nion > 0) {
1135 iparm = pbeparm->
ionr[0];
1141 Vnm_tprint( 2,
"Can't focus first calculation!\n");
1153 pbeparm->
sdie, sparm, focusFlag, pbeparm->
sdens,
1158 Vnm_tprint(0,
"Setting up PDE object...\n");
1163 mgparm->
method = (mgparm->
useAqua == 1) ? VSOL_NewtonAqua : VSOL_Newton;
1169 mgparm->
method = (mgparm->
useAqua == 1) ? VSOL_CGMGAqua : VSOL_MG;
1173 Vnm_tprint(2,
"Sorry, LRPBE isn't supported with the MG solver!\n");
1176 Vnm_tprint(2,
"Sorry, NRPBE isn't supported with the MG solver!\n");
1183 pbe[icalc]->smsize = pbeparm->
smsize;
1184 pbe[icalc]->smvolume = pbeparm->
smvolume;
1185 pbe[icalc]->ipkey = pmgp[icalc]->ipkey;
1189 Vnm_tprint(2,
"Error! Unknown PBE type (%d)!\n", pbeparm->
pbetype);
1192 Vnm_tprint(0,
"Setting PDE center to local center...\n");
1193 pmgp[icalc]->bcfl = pbeparm->
bcfl;
1194 pmgp[icalc]->xcent = realCenter[0];
1195 pmgp[icalc]->ycent = realCenter[1];
1196 pmgp[icalc]->zcent = realCenter[2];
1200 Vnm_tprint( 2,
"Can't focus first calculation!\n");
1205 pmg[icalc] =
Vpmg_ctor(pmgp[icalc], pbe[icalc], 1, pmg[icalc-1],
1211 if (icalc>0)
Vpmg_dtor(&(pmg[icalc-1]));
1212 pmg[icalc] =
Vpmg_ctor(pmgp[icalc], pbe[icalc], 0, VNULL, mgparm,
PCE_NO);
1220 theDielXMap = dielXMap[pbeparm->
dielMapID-1];
1222 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1229 theDielYMap = dielYMap[pbeparm->
dielMapID-1];
1231 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1238 theDielZMap = dielZMap[pbeparm->
dielMapID-1];
1240 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1247 theKappaMap = kappaMap[pbeparm->
kappaMapID-1];
1249 Vnm_print(2,
"Error! %d is not a valid kappa map ID!\n",
1256 thePotMap = potMap[pbeparm->
potMapID-1];
1258 Vnm_print(2,
"Error! %d is not a valid potential map ID!\n",
1267 Vnm_print(2,
"Error! %d is not a valid charge map ID!\n",
1274 Vnm_print(2,
"Warning: You specified 'bcfl map' in the input file, but no potential map was found.\n");
1275 Vnm_print(2,
" You must specify 'usemap pot' statement in the APBS input file!\n");
1276 Vnm_print(2,
"Bailing out ...\n");
1289 Vnm_print(2,
"initMG: problems setting up coefficients (fillco)!\n");
1295 Vnm_tprint(1,
" Debye length: %g A\n",
Vpbe_getDeblen(pbe[icalc]));
1302 bytesTotal = Vmem_bytesTotal();
1303 highWater = Vmem_highWaterTotal();
1306 Vnm_tprint( 1,
" Current memory usage: %4.3f MB total, \ 1307 %4.3f MB high water\n", (
double)(bytesTotal)/(1024.*1024.),
1308 (
double)(highWater)/(1024.*1024.));
1316 Vpmgp *pmgp[NOSH_MAXCALC],
Vpmg *pmg[NOSH_MAXCALC]) {
1321 Vnm_tprint(1,
"Destroying multigrid structures.\n");
1334 for(i=0;i<nosh->
ncalc;i++){
1351 if (nosh != VNULL) {
1352 if (nosh->
bogus)
return 1;
1360 Vnm_print(2,
" Error during PDE solution!\n");
1364 Vnm_tprint( 1,
" Skipping solve for mg-dummy run; zeroing \ 1369 for (i=0; i<nx*ny*
nz; i++) pmg->
u[i] = 0.0;
1386 if (nosh->
bogus)
return 1;
1389 for (j=0; j<3; j++) {
1394 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part center = (%g, %g, %g)\n",
1400 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part lower corner = (%g, %g, %g)\n",
1401 __FILE__, __LINE__, partMin[0], partMin[1], partMin[2]);
1402 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part upper corner = (%g, %g, %g)\n",
1404 partMax[0], partMax[1], partMax[2]);
1407 for (j=0; j<3; j++) {
1408 partMin[j] = mgparm->
center[j] - 0.5*mgparm->
glen[j];
1409 partMax[j] = mgparm->
center[j] + 0.5*mgparm->
glen[j];
1450 if (nosh->
bogus == 0) {
1453 Vnm_tprint( 1,
" Total electrostatic energy = %1.12E kJ/mol\n",
1456 }
else *totEnergy = 0;
1464 Vnm_tprint( 1,
" Total electrostatic energy = %1.12E \ 1466 Vnm_tprint( 1,
" Fixed charge energy = %g kJ/mol\n",
1468 Vnm_tprint( 1,
" Mobile charge energy = %g kJ/mol\n",
1470 Vnm_tprint( 1,
" Dielectric energy = %g kJ/mol\n",
1472 Vnm_tprint( 1,
" Per-atom energies:\n");
1479 Vnm_tprint( 1,
" Atom %d: %1.12E kJ/mol\n", i,
1483 }
else *nenergy = 0;
1509 Vnm_tprint( 1,
" Calculating forces...\n");
1516 for (j=0; j<3; j++) {
1517 (*atomForce)[0].
qfForce[j] = 0;
1518 (*atomForce)[0].ibForce[j] = 0;
1519 (*atomForce)[0].dbForce[j] = 0;
1522 if (nosh->
bogus == 0) {
1527 for (k=0; k<3; k++) {
1533 for (k=0; k<3; k++) {
1534 (*atomForce)[0].qfForce[k] += qfForce[k];
1535 (*atomForce)[0].ibForce[k] += ibForce[k];
1536 (*atomForce)[0].dbForce[k] += dbForce[k];
1540 Vnm_tprint( 1,
" Printing net forces for molecule %d (kJ/mol/A)\n",
1542 Vnm_tprint( 1,
" Legend:\n");
1543 Vnm_tprint( 1,
" qf -- fixed charge force\n");
1544 Vnm_tprint( 1,
" db -- dielectric boundary force\n");
1545 Vnm_tprint( 1,
" ib -- ionic boundary force\n");
1546 Vnm_tprint( 1,
" qf %4.3e %4.3e %4.3e\n",
1550 Vnm_tprint( 1,
" ib %4.3e %4.3e %4.3e\n",
1554 Vnm_tprint( 1,
" db %4.3e %4.3e %4.3e\n",
1561 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
1564 Vnm_tprint( 1,
" Printing per-atom forces for molecule %d (kJ/mol/A)\n",
1566 Vnm_tprint( 1,
" Legend:\n");
1567 Vnm_tprint( 1,
" tot n -- total force for atom n\n");
1568 Vnm_tprint( 1,
" qf n -- fixed charge force for atom n\n");
1569 Vnm_tprint( 1,
" db n -- dielectric boundary force for atom n\n");
1570 Vnm_tprint( 1,
" ib n -- ionic boundary force for atom n\n");
1573 if (nosh->
bogus == 0) {
1581 for (k=0; k<3; k++) {
1582 (*atomForce)[j].qfForce[k] = 0;
1583 (*atomForce)[j].ibForce[k] = 0;
1584 (*atomForce)[j].dbForce[k] = 0;
1588 Vnm_tprint( 1,
"mgF tot %d %4.3e %4.3e %4.3e\n", j,
1590 *((*atomForce)[j].qfForce[0]+(*atomForce)[j].ibForce[0]+
1591 (*atomForce)[j].dbForce[0]),
1593 *((*atomForce)[j].qfForce[1]+(*atomForce)[j].ibForce[1]+
1594 (*atomForce)[j].dbForce[1]),
1596 *((*atomForce)[j].qfForce[2]+(*atomForce)[j].ibForce[2]+
1597 (*atomForce)[j].dbForce[2]));
1598 Vnm_tprint( 1,
"mgF qf %d %4.3e %4.3e %4.3e\n", j,
1600 *(*atomForce)[j].qfForce[0],
1602 *(*atomForce)[j].qfForce[1],
1604 *(*atomForce)[j].qfForce[2]);
1605 Vnm_tprint( 1,
"mgF ib %d %4.3e %4.3e %4.3e\n", j,
1607 *(*atomForce)[j].ibForce[0],
1609 *(*atomForce)[j].ibForce[1],
1611 *(*atomForce)[j].ibForce[2]);
1612 Vnm_tprint( 1,
"mgF db %d %4.3e %4.3e %4.3e\n", j,
1614 *(*atomForce)[j].dbForce[0],
1616 *(*atomForce)[j].dbForce[1],
1618 *(*atomForce)[j].dbForce[2]);
1631 Vnm_tprint(1,
"No energy arrays to destroy.\n");
1642 Vnm_tprint(1,
"Destroying force arrays.\n");
1645 for (i=0; i<nosh->
ncalc; i++) {
1647 if (nforce[i] > 0) Vmem_free(mem, nforce[i],
sizeof(
AtomForce),
1648 (
void **)&(atomForce[i]));
1655 char writematstem[VMAX_ARGLEN];
1656 char outpath[VMAX_ARGLEN];
1660 if (nosh->
bogus)
return 1;
1663 strlenmax = VMAX_ARGLEN-14;
1665 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1667 Vnm_tprint(2,
" Not writing matrix!\n");
1670 sprintf(writematstem,
"%s-PE%d", pbeparm->
writematstem, rank);
1672 strlenmax = (int)(VMAX_ARGLEN)-1;
1674 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1676 Vnm_tprint(2,
" Not writing matrix!\n");
1687 strlenmax = VMAX_ARGLEN-5;
1689 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1691 Vnm_tprint(2,
" Not writing matrix!\n");
1694 sprintf(outpath,
"%s.%s", writematstem,
"mat");
1700 Vnm_tprint( 1,
" Writing Poisson operator matrix \ 1701 to %s...\n", outpath);
1705 Vnm_tprint( 1,
" Writing linearization of full \ 1706 Poisson-Boltzmann operator matrix to %s...\n", outpath);
1709 Vnm_tprint( 2,
" Bogus matrix specification\ 1714 Vnm_tprint(0,
" Printing operator...\n");
1733 *atomEnergy = (
double *)Vmem_malloc(pmg->
vmem, *nenergy,
sizeof(
double));
1735 for (i=0; i<*nenergy; i++) {
1746 double qfEnergy[NOSH_MAXCALC],
1747 double qmEnergy[NOSH_MAXCALC],
1748 double dielEnergy[NOSH_MAXCALC],
1749 int nenergy[NOSH_MAXCALC],
1750 double *atomEnergy[NOSH_MAXCALC],
1751 int nforce[NOSH_MAXCALC],
1756 int ielec, icalc, i, j;
1757 char *timestring = VNULL;
1760 double conversion, ltenergy, gtenergy, scalar;
1762 if (nosh->
bogus)
return 1;
1768 file = fopen(fname,
"w");
1769 if (file == VNULL) {
1770 Vnm_print(2,
"writedataFlat: Problem opening virtual socket %s\n",
1778 timestring = ctime(&now);
1779 fprintf(file,
"%s\n", timestring);
1781 for (ielec=0; ielec<nosh->
nelec;ielec++) {
1791 fprintf(file,
"elec");
1793 fprintf(file,
" name %s\n", nosh->
elecname[ielec]);
1794 }
else fprintf(file,
"\n");
1796 switch (mgparm->
type) {
1798 fprintf(file,
" mg-dummy\n");
1801 fprintf(file,
" mg-manual\n");
1804 fprintf(file,
" mg-auto\n");
1807 fprintf(file,
" mg-para\n");
1813 fprintf(file,
" mol %d\n", pbeparm->
molid);
1814 fprintf(file,
" dime %d %d %d\n", mgparm->
dime[0], mgparm->
dime[1],\
1819 fprintf(file,
" npbe\n");
1822 fprintf(file,
" lpbe\n");
1828 if (pbeparm->
nion > 0) {
1829 for (i=0; i<pbeparm->
nion; i++) {
1830 fprintf(file,
" ion %4.3f %4.3f %4.3f\n",
1835 fprintf(file,
" pdie %4.3f\n", pbeparm->
pdie);
1836 fprintf(file,
" sdie %4.3f\n", pbeparm->
sdie);
1838 switch (pbeparm->
srfm) {
1840 fprintf(file,
" srfm mol\n");
1841 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1844 fprintf(file,
" srfm smol\n");
1845 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1848 fprintf(file,
" srfm spl2\n");
1849 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1855 switch (pbeparm->
bcfl) {
1857 fprintf(file,
" bcfl zero\n");
1860 fprintf(file,
" bcfl sdh\n");
1863 fprintf(file,
" bcfl mdh\n");
1866 fprintf(file,
" bcfl focus\n");
1869 fprintf(file,
" bcfl map\n");
1872 fprintf(file,
" bcfl mem\n");
1878 fprintf(file,
" temp %4.3f\n", pbeparm->
temp);
1880 for (;icalc<=nosh->
elec2calc[ielec];icalc++){
1886 fprintf(file,
" calc\n");
1887 fprintf(file,
" id %i\n", (icalc+1));
1888 fprintf(file,
" grid %4.3f %4.3f %4.3f\n",
1890 fprintf(file,
" glen %4.3f %4.3f %4.3f\n",
1894 fprintf(file,
" totEnergy %1.12E kJ/mol\n",
1895 (totEnergy[icalc]*conversion));
1897 fprintf(file,
" totEnergy %1.12E kJ/mol\n",
1898 (totEnergy[icalc]*conversion));
1899 fprintf(file,
" qfEnergy %1.12E kJ/mol\n",
1900 (0.5*qfEnergy[icalc]*conversion));
1901 fprintf(file,
" qmEnergy %1.12E kJ/mol\n",
1902 (qmEnergy[icalc]*conversion));
1903 fprintf(file,
" dielEnergy %1.12E kJ/mol\n",
1904 (dielEnergy[icalc]*conversion));
1905 for (i=0; i<nenergy[icalc]; i++){
1906 fprintf(file,
" atom %i %1.12E kJ/mol\n", i,
1907 (0.5*atomEnergy[icalc][i]*conversion));
1913 fprintf(file,
" qfForce %1.12E %1.12E %1.12E kJ/mol/A\n",
1914 (atomForce[icalc][0].qfForce[0]*conversion),
1915 (atomForce[icalc][0].qfForce[1]*conversion),
1916 (atomForce[icalc][0].qfForce[2]*conversion));
1917 fprintf(file,
" ibForce %1.12E %1.12E %1.12E kJ/mol/A\n",
1918 (atomForce[icalc][0].ibForce[0]*conversion),
1919 (atomForce[icalc][0].ibForce[1]*conversion),
1920 (atomForce[icalc][0].ibForce[2]*conversion));
1921 fprintf(file,
" dbForce %1.12E %1.12E %1.12E kJ/mol/A\n",
1922 (atomForce[icalc][0].dbForce[0]*conversion),
1923 (atomForce[icalc][0].dbForce[1]*conversion),
1924 (atomForce[icalc][0].dbForce[2]*conversion));
1926 fprintf(file,
" end\n");
1929 fprintf(file,
"end\n");
1934 for (i=0; i<nosh->
nprint; i++) {
1938 fprintf(file,
"print energy");
1939 fprintf(file,
" %d", nosh->
printcalc[i][0]+1);
1942 if (nosh->
printop[i][j-1] == 0) fprintf(file,
" +");
1943 else if (nosh->
printop[i][j-1] == 1) fprintf(file,
" -");
1944 fprintf(file,
" %d", nosh->
printcalc[i][j]+1);
1947 fprintf(file,
"\n");
1956 if (nosh->
printop[i][j-1] == 0) scalar = 1.0;
1957 else if (nosh->
printop[i][j-1] == 1) scalar = -1.0;
1962 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
1965 fprintf(file,
" localEnergy %1.12E kJ/mol\n", \
1967 fprintf(file,
" globalEnergy %1.12E kJ/mol\nend\n", \
1979 double qfEnergy[NOSH_MAXCALC],
1980 double qmEnergy[NOSH_MAXCALC],
1981 double dielEnergy[NOSH_MAXCALC],
1982 int nenergy[NOSH_MAXCALC],
1983 double *atomEnergy[NOSH_MAXCALC],
1984 int nforce[NOSH_MAXCALC],
1989 int ielec, icalc, i, j;
1990 char *timestring = VNULL;
1994 double conversion, ltenergy, gtenergy, scalar;
1996 if (nosh->
bogus)
return 1;
2002 file = fopen(fname,
"w");
2003 if (file == VNULL) {
2004 Vnm_print(2,
"writedataXML: Problem opening virtual socket %s\n",
2009 fprintf(file,
"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
2010 fprintf(file,
"<APBS>\n");
2015 timestring = ctime(&now);
2016 for(c = timestring; *c !=
'\n'; c++);
2018 fprintf(file,
" <date>%s</date>\n", timestring);
2020 for (ielec=0; ielec<nosh->
nelec;ielec++){
2030 fprintf(file,
" <elec>\n");
2032 fprintf(file,
" <name>%s</name>\n", nosh->
elecname[ielec]);
2035 switch (mgparm->
type) {
2037 fprintf(file,
" <type>mg-dummy</type>\n");
2040 fprintf(file,
" <type>mg-manual</type>\n");
2043 fprintf(file,
" <type>mg-auto</type>\n");
2046 fprintf(file,
" <type>mg-para</type>\n");
2052 fprintf(file,
" <molid>%d</molid>\n", pbeparm->
molid);
2053 fprintf(file,
" <nx>%d</nx>\n", mgparm->
dime[0]);
2054 fprintf(file,
" <ny>%d</ny>\n", mgparm->
dime[1]);
2055 fprintf(file,
" <nz>%d</nz>\n", mgparm->
dime[2]);
2059 fprintf(file,
" <pbe>npbe</pbe>\n");
2062 fprintf(file,
" <pbe>lpbe</pbe>\n");
2068 if (pbeparm->
nion > 0) {
2069 for (i=0; i<pbeparm->
nion; i++) {
2070 fprintf(file,
" <ion>\n");
2071 fprintf(file,
" <radius>%4.3f A</radius>\n",
2073 fprintf(file,
" <charge>%4.3f A</charge>\n",
2075 fprintf(file,
" <concentration>%4.3f M</concentration>\n",
2077 fprintf(file,
" </ion>\n");
2082 fprintf(file,
" <pdie>%4.3f</pdie>\n", pbeparm->
pdie);
2083 fprintf(file,
" <sdie>%4.3f</sdie>\n", pbeparm->
sdie);
2085 switch (pbeparm->
srfm) {
2087 fprintf(file,
" <srfm>mol</srfm>\n");
2088 fprintf(file,
" <srad>%4.3f</srad>\n", pbeparm->
srad);
2091 fprintf(file,
" <srfm>smol</srfm>\n");
2092 fprintf(file,
" <srad>%4.3f</srad>\n", pbeparm->
srad);
2095 fprintf(file,
" <srfm>spl2</srfm>\n");
2101 switch (pbeparm->
bcfl) {
2103 fprintf(file,
" <bcfl>zero</bcfl>\n");
2106 fprintf(file,
" <bcfl>sdh</bcfl>\n");
2109 fprintf(file,
" <bcfl>mdh</bcfl>\n");
2112 fprintf(file,
" <bcfl>focus</bcfl>\n");
2115 fprintf(file,
" <bcfl>map</bcfl>\n");
2118 fprintf(file,
" <bcfl>mem</bcfl>\n");
2124 fprintf(file,
" <temp>%4.3f K</temp>\n", pbeparm->
temp);
2126 for (;icalc<=nosh->
elec2calc[ielec];icalc++){
2132 fprintf(file,
" <calc>\n");
2133 fprintf(file,
" <id>%i</id>\n", (icalc+1));
2134 fprintf(file,
" <hx>%4.3f A</hx>\n", mgparm->
grid[0]);
2135 fprintf(file,
" <hy>%4.3f A</hy>\n", mgparm->
grid[1]);
2136 fprintf(file,
" <hz>%4.3f A</hz>\n", mgparm->
grid[2]);
2137 fprintf(file,
" <xlen>%4.3f A</xlen>\n", mgparm->
glen[0]);
2138 fprintf(file,
" <ylen>%4.3f A</ylen>\n", mgparm->
glen[1]);
2139 fprintf(file,
" <zlen>%4.3f A</zlen>\n", mgparm->
glen[2]);
2142 fprintf(file,
" <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2143 (totEnergy[icalc]*conversion));
2145 fprintf(file,
" <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2146 (totEnergy[icalc]*conversion));
2147 fprintf(file,
" <qfEnergy>%1.12E kJ/mol</qfEnergy>\n",
2148 (0.5*qfEnergy[icalc]*conversion));
2149 fprintf(file,
" <qmEnergy>%1.12E kJ/mol</qmEnergy>\n",
2150 (qmEnergy[icalc]*conversion));
2151 fprintf(file,
" <dielEnergy>%1.12E kJ/mol</dielEnergy>\n",
2152 (dielEnergy[icalc]*conversion));
2153 for (i=0; i<nenergy[icalc]; i++){
2154 fprintf(file,
" <atom>\n");
2155 fprintf(file,
" <id>%i</id>\n", i+1);
2156 fprintf(file,
" <energy>%1.12E kJ/mol</energy>\n",
2157 (0.5*atomEnergy[icalc][i]*conversion));
2158 fprintf(file,
" </atom>\n");
2164 fprintf(file,
" <qfforce_x>%1.12E</qfforce_x>\n",
2165 atomForce[icalc][0].qfForce[0]*conversion);
2166 fprintf(file,
" <qfforce_y>%1.12E</qfforce_y>\n",
2167 atomForce[icalc][0].qfForce[1]*conversion);
2168 fprintf(file,
" <qfforce_z>%1.12E</qfforce_z>\n",
2169 atomForce[icalc][0].qfForce[2]*conversion);
2170 fprintf(file,
" <ibforce_x>%1.12E</ibforce_x>\n",
2171 atomForce[icalc][0].ibForce[0]*conversion);
2172 fprintf(file,
" <ibforce_y>%1.12E</ibforce_y>\n",
2173 atomForce[icalc][0].ibForce[1]*conversion);
2174 fprintf(file,
" <ibforce_z>%1.12E</ibforce_z>\n",
2175 atomForce[icalc][0].ibForce[2]*conversion);
2176 fprintf(file,
" <dbforce_x>%1.12E</dbforce_x>\n",
2177 atomForce[icalc][0].dbForce[0]*conversion);
2178 fprintf(file,
" <dbforce_y>%1.12E</dbforce_y>\n",
2179 atomForce[icalc][0].dbForce[1]*conversion);
2180 fprintf(file,
" <dbforce_z>%1.12E</dbforce_z>\n",
2181 atomForce[icalc][0].dbForce[2]*conversion);
2184 fprintf(file,
" </calc>\n");
2187 fprintf(file,
" </elec>\n");
2192 for (i=0; i<nosh->
nprint; i++) {
2196 fprintf(file,
" <printEnergy>\n");
2197 fprintf(file,
" <equation>%d", nosh->
printcalc[i][0]+1);
2200 if (nosh->
printop[i][j-1] == 0) fprintf(file,
" +");
2201 else if (nosh->
printop[i][j-1] == 1) fprintf(file,
" -");
2202 fprintf(file,
" %d", nosh->
printcalc[i][j] +1);
2205 fprintf(file,
"</equation>\n");
2214 if (nosh->
printop[i][j-1] == 0) scalar = 1.0;
2215 else if (nosh->
printop[i][j-1] == 1) scalar = -1.0;
2220 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2221 fprintf(file,
" <localEnergy>%1.12E kJ/mol</localEnergy>\n", \
2223 fprintf(file,
" <globalEnergy>%1.12E kJ/mol</globalEnergy>\n", \
2226 fprintf(file,
" </printEnergy>\n");
2231 fprintf(file,
"</APBS>\n");
2243 char writestem[VMAX_ARGLEN];
2244 char outpath[VMAX_ARGLEN];
2264 if (nosh->
bogus)
return 1;
2266 for (i=0; i<pbeparm->
numwrite; i++) {
2279 Vnm_tprint(1,
" Writing charge distribution to ");
2283 xmin = xcent - 0.5*(nx-1)*hx;
2284 ymin = ycent - 0.5*(ny-1)*hy;
2285 zmin = zcent - 0.5*(nz-1)*hzed;
2288 sprintf(title,
"CHARGE DISTRIBUTION (e)");
2293 Vnm_tprint(1,
" Writing potential to ");
2297 xmin = xcent - 0.5*(nx-1)*hx;
2298 ymin = ycent - 0.5*(ny-1)*hy;
2299 zmin = zcent - 0.5*(nz-1)*hzed;
2302 sprintf(title,
"POTENTIAL (kT/e)");
2307 Vnm_tprint(1,
" Writing molecular accessibility to ");
2311 xmin = xcent - 0.5*(nx-1)*hx;
2312 ymin = ycent - 0.5*(ny-1)*hy;
2313 zmin = zcent - 0.5*(nz-1)*hzed;
2317 "SOLVENT ACCESSIBILITY -- MOLECULAR (%4.3f PROBE)",
2323 Vnm_tprint(1,
" Writing spline-based accessibility to ");
2327 xmin = xcent - 0.5*(nx-1)*hx;
2328 ymin = ycent - 0.5*(ny-1)*hy;
2329 zmin = zcent - 0.5*(nz-1)*hzed;
2333 "SOLVENT ACCESSIBILITY -- SPLINE (%4.3f WINDOW)",
2339 Vnm_tprint(1,
" Writing van der Waals accessibility to ");
2343 xmin = xcent - 0.5*(nx-1)*hx;
2344 ymin = ycent - 0.5*(ny-1)*hy;
2345 zmin = zcent - 0.5*(nz-1)*hzed;
2348 sprintf(title,
"SOLVENT ACCESSIBILITY -- VAN DER WAALS");
2353 Vnm_tprint(1,
" Writing ion accessibility to ");
2357 xmin = xcent - 0.5*(nx-1)*hx;
2358 ymin = ycent - 0.5*(ny-1)*hy;
2359 zmin = zcent - 0.5*(nz-1)*hzed;
2363 "ION ACCESSIBILITY -- SPLINE (%4.3f RADIUS)",
2369 Vnm_tprint(1,
" Writing potential Laplacian to ");
2373 xmin = xcent - 0.5*(nx-1)*hx;
2374 ymin = ycent - 0.5*(ny-1)*hy;
2375 zmin = zcent - 0.5*(nz-1)*hzed;
2379 "POTENTIAL LAPLACIAN (kT/e/A^2)");
2384 Vnm_tprint(1,
" Writing energy density to ");
2388 xmin = xcent - 0.5*(nx-1)*hx;
2389 ymin = ycent - 0.5*(ny-1)*hy;
2390 zmin = zcent - 0.5*(nz-1)*hzed;
2393 sprintf(title,
"ENERGY DENSITY (kT/e/A)^2");
2398 Vnm_tprint(1,
" Writing number density to ");
2402 xmin = xcent - 0.5*(nx-1)*hx;
2403 ymin = ycent - 0.5*(ny-1)*hy;
2404 zmin = zcent - 0.5*(nz-1)*hzed;
2408 "ION NUMBER DENSITY (M)");
2413 Vnm_tprint(1,
" Writing charge density to ");
2417 xmin = xcent - 0.5*(nx-1)*hx;
2418 ymin = ycent - 0.5*(ny-1)*hy;
2419 zmin = zcent - 0.5*(nz-1)*hzed;
2423 "ION CHARGE DENSITY (e_c * M)");
2428 Vnm_tprint(1,
" Writing x-shifted dielectric map to ");
2432 xmin = xcent - 0.5*(nx-1)*hx;
2433 ymin = ycent - 0.5*(ny-1)*hy;
2434 zmin = zcent - 0.5*(nz-1)*hzed;
2438 "X-SHIFTED DIELECTRIC MAP");
2443 Vnm_tprint(1,
" Writing y-shifted dielectric map to ");
2447 xmin = xcent - 0.5*(nx-1)*hx;
2448 ymin = ycent - 0.5*(ny-1)*hy;
2449 zmin = zcent - 0.5*(nz-1)*hzed;
2453 "Y-SHIFTED DIELECTRIC MAP");
2458 Vnm_tprint(1,
" Writing z-shifted dielectric map to ");
2462 xmin = xcent - 0.5*(nx-1)*hx;
2463 ymin = ycent - 0.5*(ny-1)*hy;
2464 zmin = zcent - 0.5*(nz-1)*hzed;
2468 "Z-SHIFTED DIELECTRIC MAP");
2473 Vnm_tprint(1,
" Writing kappa map to ");
2477 xmin = xcent - 0.5*(nx-1)*hx;
2478 ymin = ycent - 0.5*(ny-1)*hy;
2479 zmin = zcent - 0.5*(nz-1)*hzed;
2488 Vnm_tprint(1,
" Writing atom potentials to ");
2492 xmin = xcent - 0.5*(nx-1)*hx;
2493 ymin = ycent - 0.5*(ny-1)*hy;
2494 zmin = zcent - 0.5*(nz-1)*hzed;
2502 Vnm_tprint(2,
"Invalid data type for writing!\n");
2508 sprintf(writestem,
"%s-PE%d", pbeparm->
writestem[i], rank);
2513 sprintf(writestem,
"%s", pbeparm->
writestem[i]);
2520 sprintf(outpath,
"%s.%s", writestem,
"dx");
2521 Vnm_tprint(1,
"%s\n", outpath);
2522 grid =
Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2530 sprintf(outpath,
"%s.%s", writestem,
"ucd");
2531 Vnm_tprint(1,
"%s\n", outpath);
2532 Vnm_tprint(2,
"Sorry, AVS format isn't supported for \ 2533 uniform meshes yet!\n");
2537 sprintf(outpath,
"%s.%s", writestem,
"mcsf");
2538 Vnm_tprint(1,
"%s\n", outpath);
2539 Vnm_tprint(2,
"Sorry, MCSF format isn't supported for \ 2540 uniform meshes yet!\n");
2544 sprintf(outpath,
"%s.%s", writestem,
"grd");
2545 Vnm_tprint(1,
"%s\n", outpath);
2546 grid =
Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2554 sprintf(outpath,
"%s.%s", writestem,
"dx.gz");
2555 Vnm_tprint(1,
"%s\n", outpath);
2556 grid =
Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2563 sprintf(outpath,
"%s.%s", writestem,
"txt");
2564 Vnm_tprint(1,
"%s\n", outpath);
2565 Vnm_print(0,
"routines: Opening virtual socket...\n");
2566 sock = Vio_ctor(
"FILE",
"ASC",VNULL,outpath,
"w");
2567 if (sock == VNULL) {
2568 Vnm_print(2,
"routines: Problem opening virtual socket %s\n",
2572 if (Vio_connect(sock, 0) < 0) {
2573 Vnm_print(2,
"routines: Problem connecting virtual socket %s\n",
2577 Vio_printf(sock,
"# Data from %s\n", PACKAGE_STRING);
2578 Vio_printf(sock,
"# \n");
2579 Vio_printf(sock,
"# %s\n", title);
2580 Vio_printf(sock,
"# \n");
2582 for(i=0;i<natoms;i++)
2583 Vio_printf(sock,
"%12.6e\n", pmg->
rwork[i]);
2586 Vnm_tprint(2,
"Bogus data format (%d)!\n",
2612 Vnm_tprint( 2,
" No energy available in Calculation %d\n", calcid+1);
2615 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++){
2618 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2619 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2640 Vnm_tprint( 2,
"Warning: The 'energy' print keyword is deprecated.\n" \
2641 " Use elecEnergy for electrostatics energy calcs.\n\n");
2644 Vnm_tprint( 1,
"print energy %d ", nosh->
printcalc[iprint][0]+1);
2646 Vnm_tprint( 1,
"print energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2649 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2650 if (nosh->
printop[iprint][iarg-1] == 0)
2651 Vnm_tprint(1,
"+ ");
2652 else if (nosh->
printop[iprint][iarg-1] == 1)
2653 Vnm_tprint(1,
"- ");
2655 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2660 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2662 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2666 Vnm_tprint(1,
"end\n");
2672 Vnm_tprint( 2,
" Didn't calculate energy in Calculation \ 2676 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2679 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2680 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2686 Vnm_tprint( 1,
" Local net energy (PE %d) = %1.12E kJ/mol\n",
2687 Vcom_rank(com), ltenergy);
2688 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
2689 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2690 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2709 Vnm_tprint( 1,
"\nprint energy %d ", nosh->
printcalc[iprint][0]+1);
2711 Vnm_tprint( 1,
"\nprint energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2714 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2715 if (nosh->
printop[iprint][iarg-1] == 0)
2716 Vnm_tprint(1,
"+ ");
2717 else if (nosh->
printop[iprint][iarg-1] == 1)
2718 Vnm_tprint(1,
"- ");
2720 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2725 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2727 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2731 Vnm_tprint(1,
"end\n");
2737 Vnm_tprint( 2,
" Didn't calculate energy in Calculation \ 2741 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2744 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2745 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2751 Vnm_tprint( 1,
" Local net energy (PE %d) = %1.12E kJ/mol\n",
2752 Vcom_rank(com), ltenergy);
2753 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
2754 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2755 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2772 Vnm_tprint( 1,
"\nprint APOL energy %d ", nosh->
printcalc[iprint][0]+1);
2774 Vnm_tprint( 1,
"\nprint APOL energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2777 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2778 if (nosh->
printop[iprint][iarg-1] == 0)
2779 Vnm_tprint(1,
"+ ");
2780 else if (nosh->
printop[iprint][iarg-1] == 1)
2781 Vnm_tprint(1,
"- ");
2783 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2788 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2790 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2794 Vnm_tprint(1,
"end\n");
2802 Vnm_tprint( 2,
" Didn't calculate energy in Calculation #%d\n", calcid+1);
2805 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2810 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2811 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2813 gtenergy += (scalar * ((apolparm->
gamma*apolparm->
sasa) +
2819 Vnm_tprint( 1,
" Global net APOL energy = %1.12E kJ/mol\n", gtenergy);
2843 Vnm_tprint( 2,
"Warning: The 'force' print keyword is deprecated.\n" \
2844 " Use elecForce for electrostatics force calcs.\n\n");
2847 Vnm_tprint( 1,
"print force %d ", nosh->
printcalc[iprint][0]+1);
2849 Vnm_tprint( 1,
"print force %d (%s) ", nosh->
printcalc[iprint][0]+1,
2852 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2853 if (nosh->
printop[iprint][iarg-1] == 0)
2854 Vnm_tprint(1,
"+ ");
2855 else if (nosh->
printop[iprint][iarg-1] == 1)
2856 Vnm_tprint(1,
"- ");
2858 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2863 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2865 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2869 Vnm_tprint(1,
"end\n");
2874 refnforce = nforce[calcid];
2876 if (refcalcforce ==
PCF_NO) {
2877 Vnm_tprint( 2,
" Didn't calculate force in calculation \ 2881 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2884 Vnm_tprint(2,
" Inconsistent calcforce declarations in \ 2889 if (nforce[calcid] != refnforce) {
2890 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \ 2903 aforce = atomForce[calcid];
2909 for (ivc=0; ivc<3; ivc++) {
2918 for (ifr=0; ifr<refnforce; ifr++) {
2919 for (ivc=0; ivc<3; ivc++) {
2931 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2934 aforce = atomForce[calcid];
2936 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
2937 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2942 for (ivc=0; ivc<3; ivc++) {
2951 for (ifr=0; ifr<refnforce; ifr++) {
2952 for (ivc=0; ivc<3; ivc++) {
2964 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
2965 for (ifr=0; ifr<refnforce; ifr++) {
2966 Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
2967 Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
2968 Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
2973 Vnm_tprint( 1,
" Local net fixed charge force = \ 2974 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
2975 lforce[0].qfForce[1], lforce[0].qfForce[2]);
2976 Vnm_tprint( 1,
" Local net ionic boundary force = \ 2977 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
2978 lforce[0].ibForce[1], lforce[0].ibForce[2]);
2979 Vnm_tprint( 1,
" Local net dielectric boundary force = \ 2980 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
2981 lforce[0].dbForce[1], lforce[0].dbForce[2]);
2983 for (ifr=0; ifr<refnforce; ifr++) {
2984 Vnm_tprint( 1,
" Local fixed charge force \ 2985 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
2986 lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
2987 Vnm_tprint( 1,
" Local ionic boundary force \ 2988 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
2989 lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
2990 Vnm_tprint( 1,
" Local dielectric boundary force \ 2991 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
2992 lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
2998 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A).\n");
2999 Vnm_tprint( 1,
" Legend:\n");
3000 Vnm_tprint( 1,
" tot -- Total force\n");
3001 Vnm_tprint( 1,
" qf -- Fixed charge force\n");
3002 Vnm_tprint( 1,
" db -- Dielectric boundary force\n");
3003 Vnm_tprint( 1,
" ib -- Ionic boundary force\n");
3005 for (ivc=0; ivc<3; ivc++) {
3011 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3012 totforce[1], totforce[2]);
3013 Vnm_tprint( 1,
" qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3014 gforce[0].qfForce[1], gforce[0].qfForce[2]);
3015 Vnm_tprint( 1,
" ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3016 gforce[0].ibForce[1], gforce[0].ibForce[2]);
3017 Vnm_tprint( 1,
" db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3018 gforce[0].dbForce[1], gforce[0].dbForce[2]);
3022 Vnm_tprint( 1,
" Printing per-atom forces (kJ/mol/A).\n");
3023 Vnm_tprint( 1,
" Legend:\n");
3024 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3025 Vnm_tprint( 1,
" qf n -- Fixed charge force for atom n\n");
3026 Vnm_tprint( 1,
" db n -- Dielectric boundary force for atom n\n");
3027 Vnm_tprint( 1,
" ib n -- Ionic boundary force for atom n\n");
3028 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3034 for (ifr=0; ifr<refnforce; ifr++) {
3035 Vnm_tprint( 1,
" qf %d %1.12E %1.12E %1.12E\n", ifr,
3036 gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3037 gforce[ifr].qfForce[2]);
3038 Vnm_tprint( 1,
" ib %d %1.12E %1.12E %1.12E\n", ifr,
3039 gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3040 gforce[ifr].ibForce[2]);
3041 Vnm_tprint( 1,
" db %d %1.12E %1.12E %1.12E\n", ifr,
3042 gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3043 gforce[ifr].dbForce[2]);
3044 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3045 (gforce[ifr].dbForce[0] \
3046 + gforce[ifr].ibForce[0] +
3047 gforce[ifr].qfForce[0]),
3048 (gforce[ifr].dbForce[1] \
3049 + gforce[ifr].ibForce[1] +
3050 gforce[ifr].qfForce[1]),
3051 (gforce[ifr].dbForce[2] \
3052 + gforce[ifr].ibForce[2] +
3053 gforce[ifr].qfForce[2]));
3054 for (ivc=0; ivc<3; ivc++) {
3055 totforce[ivc] += (gforce[ifr].
dbForce[ivc] \
3060 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3061 totforce[1], totforce[2]);
3064 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3065 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3090 Vnm_tprint( 1,
"print force %d ", nosh->
printcalc[iprint][0]+1);
3092 Vnm_tprint( 1,
"print force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3095 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3096 if (nosh->
printop[iprint][iarg-1] == 0)
3097 Vnm_tprint(1,
"+ ");
3098 else if (nosh->
printop[iprint][iarg-1] == 1)
3099 Vnm_tprint(1,
"- ");
3101 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3106 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3108 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3112 Vnm_tprint(1,
"end\n");
3117 refnforce = nforce[calcid];
3119 if (refcalcforce ==
PCF_NO) {
3120 Vnm_tprint( 2,
" Didn't calculate force in calculation \ 3124 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3127 Vnm_tprint(2,
" Inconsistent calcforce declarations in \ 3132 if (nforce[calcid] != refnforce) {
3133 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \ 3146 aforce = atomForce[calcid];
3152 for (ivc=0; ivc<3; ivc++) {
3161 for (ifr=0; ifr<refnforce; ifr++) {
3162 for (ivc=0; ivc<3; ivc++) {
3174 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3177 aforce = atomForce[calcid];
3179 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3180 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3185 for (ivc=0; ivc<3; ivc++) {
3194 for (ifr=0; ifr<refnforce; ifr++) {
3195 for (ivc=0; ivc<3; ivc++) {
3207 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
3208 for (ifr=0; ifr<refnforce; ifr++) {
3209 Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
3210 Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
3211 Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
3216 Vnm_tprint( 1,
" Local net fixed charge force = \ 3217 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
3218 lforce[0].qfForce[1], lforce[0].qfForce[2]);
3219 Vnm_tprint( 1,
" Local net ionic boundary force = \ 3220 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
3221 lforce[0].ibForce[1], lforce[0].ibForce[2]);
3222 Vnm_tprint( 1,
" Local net dielectric boundary force = \ 3223 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
3224 lforce[0].dbForce[1], lforce[0].dbForce[2]);
3226 for (ifr=0; ifr<refnforce; ifr++) {
3227 Vnm_tprint( 1,
" Local fixed charge force \ 3228 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
3229 lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
3230 Vnm_tprint( 1,
" Local ionic boundary force \ 3231 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
3232 lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
3233 Vnm_tprint( 1,
" Local dielectric boundary force \ 3234 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
3235 lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
3241 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A).\n");
3242 Vnm_tprint( 1,
" Legend:\n");
3243 Vnm_tprint( 1,
" tot -- Total force\n");
3244 Vnm_tprint( 1,
" qf -- Fixed charge force\n");
3245 Vnm_tprint( 1,
" db -- Dielectric boundary force\n");
3246 Vnm_tprint( 1,
" ib -- Ionic boundary force\n");
3248 for (ivc=0; ivc<3; ivc++) {
3254 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3255 totforce[1], totforce[2]);
3256 Vnm_tprint( 1,
" qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3257 gforce[0].qfForce[1], gforce[0].qfForce[2]);
3258 Vnm_tprint( 1,
" ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3259 gforce[0].ibForce[1], gforce[0].ibForce[2]);
3260 Vnm_tprint( 1,
" db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3261 gforce[0].dbForce[1], gforce[0].dbForce[2]);
3265 Vnm_tprint( 1,
" Printing per-atom forces (kJ/mol/A).\n");
3266 Vnm_tprint( 1,
" Legend:\n");
3267 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3268 Vnm_tprint( 1,
" qf n -- Fixed charge force for atom n\n");
3269 Vnm_tprint( 1,
" db n -- Dielectric boundary force for atom n\n");
3270 Vnm_tprint( 1,
" ib n -- Ionic boundary force for atom n\n");
3271 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3277 for (ifr=0; ifr<refnforce; ifr++) {
3278 Vnm_tprint( 1,
" qf %d %1.12E %1.12E %1.12E\n", ifr,
3279 gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3280 gforce[ifr].qfForce[2]);
3281 Vnm_tprint( 1,
" ib %d %1.12E %1.12E %1.12E\n", ifr,
3282 gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3283 gforce[ifr].ibForce[2]);
3284 Vnm_tprint( 1,
" db %d %1.12E %1.12E %1.12E\n", ifr,
3285 gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3286 gforce[ifr].dbForce[2]);
3287 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3288 (gforce[ifr].dbForce[0] \
3289 + gforce[ifr].ibForce[0] +
3290 gforce[ifr].qfForce[0]),
3291 (gforce[ifr].dbForce[1] \
3292 + gforce[ifr].ibForce[1] +
3293 gforce[ifr].qfForce[1]),
3294 (gforce[ifr].dbForce[2] \
3295 + gforce[ifr].ibForce[2] +
3296 gforce[ifr].qfForce[2]));
3297 for (ivc=0; ivc<3; ivc++) {
3298 totforce[ivc] += (gforce[ifr].
dbForce[ivc] \
3303 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3304 totforce[1], totforce[2]);
3307 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3308 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3335 Vnm_tprint( 1,
"\nprint APOL force %d ", nosh->
printcalc[iprint][0]+1);
3337 Vnm_tprint( 1,
"\nprint APOL force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3340 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3341 if (nosh->
printop[iprint][iarg-1] == 0)
3342 Vnm_tprint(1,
"+ ");
3343 else if (nosh->
printop[iprint][iarg-1] == 1)
3344 Vnm_tprint(1,
"- ");
3346 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3351 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3353 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3357 Vnm_tprint(1,
"end\n");
3362 refnforce = nforce[calcid];
3364 if (refcalcforce ==
ACF_NO) {
3365 Vnm_tprint( 2,
" Didn't calculate force in calculation \ 3369 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3372 Vnm_tprint(2,
" Inconsistent calcforce declarations in \ 3377 if (nforce[calcid] != refnforce) {
3378 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \ 3391 aforce = atomForce[calcid];
3397 for (ivc=0; ivc<3; ivc++) {
3403 for (ifr=0; ifr<refnforce; ifr++) {
3404 for (ivc=0; ivc<3; ivc++) {
3413 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3416 aforce = atomForce[calcid];
3418 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3419 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3424 for (ivc=0; ivc<3; ivc++) {
3430 for (ifr=0; ifr<refnforce; ifr++) {
3431 for (ivc=0; ivc<3; ivc++) {
3440 Vnm_tprint( 0,
"printForce: Performing global reduction (sum)\n");
3441 for (ifr=0; ifr<refnforce; ifr++) {
3442 Vcom_reduce(com, lforce[ifr].sasaForce, gforce[ifr].sasaForce, 3, 2, 0);
3443 Vcom_reduce(com, lforce[ifr].savForce, gforce[ifr].savForce, 3, 2, 0);
3444 Vcom_reduce(com, lforce[ifr].wcaForce, gforce[ifr].wcaForce, 3, 2, 0);
3448 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A)\n");
3449 Vnm_tprint( 1,
" Legend:\n");
3450 Vnm_tprint( 1,
" tot -- Total force\n");
3451 Vnm_tprint( 1,
" sasa -- SASA force\n");
3452 Vnm_tprint( 1,
" sav -- SAV force\n");
3453 Vnm_tprint( 1,
" wca -- WCA force\n\n");
3455 for (ivc=0; ivc<3; ivc++) {
3461 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3462 totforce[1], totforce[2]);
3463 Vnm_tprint( 1,
" sasa %1.12E %1.12E %1.12E\n", gforce[0].sasaForce[0],
3464 gforce[0].sasaForce[1], gforce[0].sasaForce[2]);
3465 Vnm_tprint( 1,
" sav %1.12E %1.12E %1.12E\n", gforce[0].savForce[0],
3466 gforce[0].savForce[1], gforce[0].savForce[2]);
3467 Vnm_tprint( 1,
" wca %1.12E %1.12E %1.12E\n", gforce[0].wcaForce[0],
3468 gforce[0].wcaForce[1], gforce[0].wcaForce[2]);
3472 Vnm_tprint( 1,
" Printing per atom forces (kJ/mol/A)\n");
3473 Vnm_tprint( 1,
" Legend:\n");
3474 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3475 Vnm_tprint( 1,
" sasa n -- SASA force for atom n\n");
3476 Vnm_tprint( 1,
" sav n -- SAV force for atom n\n");
3477 Vnm_tprint( 1,
" wca n -- WCA force for atom n\n");
3478 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3487 for (ifr=0; ifr<refnforce; ifr++) {
3488 Vnm_tprint( 1,
" sasa %d %1.12E %1.12E %1.12E\n", ifr,
3489 gforce[ifr].sasaForce[0], gforce[ifr].sasaForce[1],
3490 gforce[ifr].sasaForce[2]);
3491 Vnm_tprint( 1,
" sav %d %1.12E %1.12E %1.12E\n", ifr,
3492 gforce[ifr].savForce[0], gforce[ifr].savForce[1],
3493 gforce[ifr].savForce[2]);
3494 Vnm_tprint( 1,
" wca %d %1.12E %1.12E %1.12E\n", ifr,
3495 gforce[ifr].wcaForce[0], gforce[ifr].wcaForce[1],
3496 gforce[ifr].wcaForce[2]);
3497 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3498 (gforce[ifr].wcaForce[0] \
3499 + gforce[ifr].savForce[0] +
3500 gforce[ifr].sasaForce[0]),
3501 (gforce[ifr].wcaForce[1] \
3502 + gforce[ifr].savForce[1] +
3503 gforce[ifr].sasaForce[1]),
3504 (gforce[ifr].wcaForce[2] \
3505 + gforce[ifr].savForce[2] +
3506 gforce[ifr].sasaForce[2]));
3507 for (ivc=0; ivc<3; ivc++) {
3508 totforce[ivc] += (gforce[ifr].
wcaForce[ivc] \
3513 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3514 totforce[1], totforce[2]);
3517 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3518 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3528 Vfetk *fetk[NOSH_MAXCALC],
3535 Vnm_tprint(1,
"Destroying finite element structures.\n");
3538 for(i=0;i<nosh->
ncalc;i++){
3542 for (i=0; i<nosh->
nmesh; i++) {
3560 Vfetk *fetk[NOSH_MAXCALC]
3580 Vatom *atom = VNULL;
3582 Vnm_tstart(27,
"Setup timer");
3585 if (pbeparm->
useDielMap) Vnm_tprint(2,
"FEM ignoring dielectric map!\n");
3586 if (pbeparm->
useKappaMap) Vnm_tprint(2,
"FEM ignoring kappa map!\n");
3587 if (pbeparm->
useChargeMap) Vnm_tprint(2,
"FEM ignoring charge map!\n");
3590 Vnm_tprint(0,
"Re-centering mesh...\n");
3591 theMol = pbeparm->
molid-1;
3592 myalist = alist[theMol];
3593 for (j=0; j<3; j++) {
3594 if (theMol < nosh->nmol) {
3595 center[j] = (myalist)->center[j];
3597 Vnm_tprint(2,
"ERROR! Bogus molecule number (%d)!\n",
3625 Vnm_tprint(0,
"Setting up PBE object...\n");
3627 sparm = pbeparm->
swin;
3630 sparm = pbeparm->
srad;
3632 if (pbeparm->
nion > 0) {
3633 iparm = pbeparm->
ionr[0];
3639 pbeparm->
sdie, sparm, focusFlag, pbeparm->
sdens,
3644 Vnm_tprint(1,
" Debye length: %g A\n",
Vpbe_getDeblen(pbe[icalc]));
3647 Vnm_tprint(0,
"Setting up FEtk object...\n");
3654 Vnm_tprint(0,
"Setting up mesh...\n");
3657 imesh = feparm->
meshID-1;
3659 for (i=0; i<3; i++) {
3663 Vnm_print(0,
"Using mesh %d (%s) in calculation.\n", imesh+1,
3665 switch (nosh->
meshfmt[imesh]) {
3667 Vnm_tprint(2,
"DX finite element mesh input not supported yet!\n");
3670 Vnm_tprint( 2,
"UHBD finite element mesh input not supported!\n");
3673 Vnm_tprint( 2,
"AVS finite element mesh input not supported!\n");
3676 Vnm_tprint(1,
"Reading MCSF-format input finite element mesh from %s.\n",
3678 sock = Vio_ctor(
"FILE",
"ASC", VNULL, nosh->
meshpath[imesh],
"r");
3679 if (sock == VNULL) {
3680 Vnm_print(2,
"Problem opening virtual socket %s!\n",
3684 if (Vio_accept(sock, 0) < 0) {
3685 Vnm_print(2,
"Problem accepting virtual socket %s!\n",
3692 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
3698 for (i=0; i<3; i++) {
3699 center[i] = alist[theMol]->center[i];
3700 length[i] = feparm->
glen[i];
3705 vrc =
Vfetk_loadMesh(fetk[icalc], center, length, meshType, sock);
3707 Vnm_print(2,
"Error constructing finite element mesh!\n");
3711 Gem_shapeChk(fetk[icalc]->gm);
3715 for (j=0; j<2; j++) {
3719 AM_markRefine(fetk[icalc]->am, 0, -1, 0, 0.0);
3721 AM_refine(fetk[icalc]->am, 2, 0, feparm->
pkey);
3723 Gem_shapeChk(fetk[icalc]->gm);
3728 Vnm_tstop(27,
"Setup timer");
3731 bytesTotal = Vmem_bytesTotal();
3732 highWater = Vmem_highWaterTotal();
3735 Vnm_tprint( 1,
" Current memory usage: %4.3f MB total, \ 3736 %4.3f MB high water\n", (
double)(bytesTotal)/(1024.*1024.),
3737 (
double)(highWater)/(1024.*1024.));
3749 Vnm_tprint(1,
" Domain size: %g A x %g A x %g A\n",
3752 switch(feparm->
ekey) {
3754 Vnm_tprint(1,
" Per-simplex error tolerance: %g\n", feparm->
etol);
3757 Vnm_tprint(1,
" Global error tolerance: %g\n", feparm->
etol);
3760 Vnm_tprint(1,
" Fraction of simps to refine: %g\n", feparm->
etol);
3763 Vnm_tprint(2,
"Invalid ekey (%d)!\n", feparm->
ekey);
3769 Vnm_tprint(1,
" Uniform pre-solve refinement.\n");
3772 Vnm_tprint(1,
" Geometry-based pre-solve refinement.\n");
3775 Vnm_tprint(2,
"Invalid akeyPRE (%d)!\n", feparm->
akeyPRE);
3781 Vnm_tprint(1,
" Residual-based a posteriori refinement.\n");
3784 Vnm_tprint(1,
" Dual-based a posteriori refinement.\n");
3787 Vnm_tprint(1,
" Local-based a posteriori refinement.\n");
3790 Vnm_tprint(2,
"Invalid akeySOLVE (%d)!\n", feparm->
akeySOLVE);
3793 Vnm_tprint(1,
" Refinement of initial mesh to ~%d vertices\n",
3795 Vnm_tprint(1,
" Geometry-based refinment lower bound: %g A\n",
3797 Vnm_tprint(1,
" Maximum number of solve-estimate-refine cycles: %d\n",
3799 Vnm_tprint(1,
" Maximum number of vertices in mesh: %d\n",
3803 if (nosh->
bogus)
return;
3805 Vnm_tprint(1,
" HB linear solver: AM_hPcg\n");
3807 Vnm_tprint(1,
" Non-HB linear solver: ");
3808 switch (fetk[icalc]->lkey) {
3810 Vnm_print(1,
"SLU direct\n");
3813 Vnm_print(1,
"multigrid\n");
3816 Vnm_print(1,
"conjugate gradient\n");
3819 Vnm_print(1,
"BiCGStab\n");
3822 Vnm_print(1,
"???\n");
3827 Vnm_tprint(1,
" Linear solver tol.: %g\n", fetk[icalc]->ltol);
3828 Vnm_tprint(1,
" Linear solver max. iters.: %d\n", fetk[icalc]->lmax);
3829 Vnm_tprint(1,
" Linear solver preconditioner: ");
3830 switch (fetk[icalc]->lprec) {
3832 Vnm_print(1,
"identity\n");
3835 Vnm_print(1,
"diagonal\n");
3838 Vnm_print(1,
"multigrid\n");
3841 Vnm_print(1,
"???\n");
3844 Vnm_tprint(1,
" Nonlinear solver: ");
3845 switch (fetk[icalc]->nkey) {
3847 Vnm_print(1,
"newton\n");
3850 Vnm_print(1,
"incremental\n");
3853 Vnm_print(1,
"pseudo-arclength\n");
3856 Vnm_print(1,
"??? ");
3859 Vnm_tprint(1,
" Nonlinear solver tol.: %g\n", fetk[icalc]->ntol);
3860 Vnm_tprint(1,
" Nonlinear solver max. iters.: %d\n", fetk[icalc]->nmax);
3861 Vnm_tprint(1,
" Initial guess: ");
3862 switch (fetk[icalc]->gues) {
3864 Vnm_tprint(1,
"zero\n");
3867 Vnm_tprint(1,
"boundary function\n");
3870 Vnm_tprint(1,
"interpolated previous solution\n");
3873 Vnm_tprint(1,
"???\n");
3898 Vnm_tprint(1,
" Commencing uniform refinement to %d verts.\n",
3902 Vnm_tprint(1,
" Commencing geometry-based refinement to %d \ 3906 Vnm_tprint(2,
"What? You can't do a posteriori error estimation \ 3907 before you solve!\n");
3922 Vnm_tprint(1,
" Initial mesh has %d vertices\n",
3923 Gem_numVV(fetk[icalc]->gm));
3929 nverts = Gem_numVV(fetk[icalc]->gm);
3931 Vnm_tprint(1,
" Hit vertex number limit.\n");
3934 marked = AM_markRefine(fetk[icalc]->am, feparm->
akeyPRE, -1,
3937 Vnm_tprint(1,
" Marked 0 simps; hit error/size tolerance.\n");
3940 Vnm_tprint(1,
" Have %d verts, marked %d. Refining...\n", nverts,
3942 AM_refine(fetk[icalc]->am, 0, 0, feparm->
pkey);
3945 nverts = Gem_numVV(fetk[icalc]->gm);
3946 Vnm_tprint(1,
" Done refining; have %d verts.\n", nverts);
3968 (pbeparm->
pbetype == PBE_NRPBE) ||
3991 Vnm_print(2,
"SORRY! DON'T USE HB!!!\n");
3995 AM_hlSolve(fetk[icalc]->am,
meth, lkeyHB, fetk[icalc]->lmax,
3996 fetk[icalc]->ltol, fetk[icalc]->gues, fetk[icalc]->pjac);
4041 if (nosh->
bogus == 0) {
4043 (pbeparm->
pbetype==PBE_NRPBE) ||
4054 Vnm_tprint(1,
" Total electrostatic energy = %1.12E kJ/mol\n",
4061 Vnm_tprint(2,
"Error! Verbose energy evaluation not available for FEM yet!\n");
4062 Vnm_tprint(2,
"E-mail nathan.baker@pnl.gov if you want this.\n");
4087 nverts = Gem_numVV(fetk[icalc]->gm);
4088 if (nverts > feparm->
maxvert) {
4089 Vnm_tprint(1,
" Current number of vertices (%d) exceeds max (%d)!\n",
4093 Vnm_tprint(1,
" Mesh currently has %d vertices\n", nverts);
4097 Vnm_tprint(1,
" Commencing uniform refinement.\n");
4100 Vnm_tprint(1,
" Commencing geometry-based refinement.\n");
4103 Vnm_tprint(1,
" Commencing residual-based refinement.\n");
4106 Vnm_tprint(1,
" Commencing dual problem-based refinement.\n");
4109 Vnm_tprint(1,
" Commencing local-based refinement.\n.");
4112 Vnm_tprint(2,
" Error -- unknown refinement type (%d)!\n",
4119 marked = AM_markRefine(fetk[icalc]->am, feparm->
akeySOLVE, -1,
4122 Vnm_tprint(1,
" Marked 0 simps; hit error/size tolerance.\n");
4125 Vnm_tprint(1,
" Have %d verts, marked %d. Refining...\n", nverts,
4127 AM_refine(fetk[icalc]->am, 0, 0, feparm->
pkey);
4128 nverts = Gem_numVV(fetk[icalc]->gm);
4129 Vnm_tprint(1,
" Done refining; have %d verts.\n", nverts);
4131 Gem_shapeChk(fetk[icalc]->gm);
4146 char writestem[VMAX_ARGLEN];
4147 char outpath[VMAX_ARGLEN];
4153 if (nosh->
bogus)
return 1;
4158 for (i=0; i<pbeparm->
numwrite; i++) {
4166 Vnm_tprint(2,
" Sorry; can't write charge distribution for FEM!\n");
4172 Vnm_tprint(1,
" Writing potential to ");
4178 Vnm_tprint(1,
" Writing molecular accessibility to ");
4184 Vnm_tprint(1,
" Writing spline-based accessibility to ");
4190 Vnm_tprint(1,
" Writing van der Waals accessibility to ");
4196 Vnm_tprint(1,
" Writing ion accessibility to ");
4202 Vnm_tprint(2,
" Sorry; can't write charge distribution for FEM!\n");
4208 Vnm_tprint(2,
" Sorry; can't write energy density for FEM!\n");
4214 Vnm_tprint(1,
" Writing number density to ");
4220 Vnm_tprint(1,
" Writing charge density to ");
4226 Vnm_tprint(2,
" Sorry; can't write x-shifted dielectric map for FEM!\n");
4232 Vnm_tprint(2,
" Sorry; can't write y-shifted dielectric map for FEM!\n");
4238 Vnm_tprint(2,
" Sorry; can't write z-shifted dielectric map for FEM!\n");
4244 Vnm_tprint(1,
" Sorry; can't write kappa map for FEM!\n");
4250 Vnm_tprint(1,
" Sorry; can't write atom potentials for FEM!\n");
4256 Vnm_tprint(2,
"Invalid data type for writing!\n");
4261 if (!writeit)
return 0;
4265 sprintf(writestem,
"%s-PE%d", pbeparm->
writestem[i], rank);
4270 sprintf(writestem,
"%s", pbeparm->
writestem[i]);
4277 sprintf(outpath,
"%s.%s", writestem,
"dx");
4278 Vnm_tprint(1,
"%s\n", outpath);
4283 sprintf(outpath,
"%s.%s", writestem,
"ucd");
4284 Vnm_tprint(1,
"%s\n", outpath);
4289 Vnm_tprint(2,
"UHBD format not supported for FEM!\n");
4293 Vnm_tprint(2,
"MCSF format not supported yet!\n");
4297 Vnm_tprint(2,
"Bogus data format (%d)!\n",
4325 Vatom *atom = VNULL;
4368 for (i=0; i < natoms; i++) {
4374 if ((x+atomRadius) > xmax) xmax = x + atomRadius;
4375 if ((x-atomRadius) < xmin) xmin = x - atomRadius;
4376 if ((y+atomRadius) > ymax) ymax = y + atomRadius;
4377 if ((y-atomRadius) < ymin) ymin = y - atomRadius;
4378 if ((z+atomRadius) > zmax) zmax = z + atomRadius;
4379 if ((z-atomRadius) < zmin) zmin = z - atomRadius;
4380 disp[0] = (x - center[0]);
4381 disp[1] = (y - center[1]);
4382 disp[2] = (z - center[2]);
4383 dist = (disp[0]*disp[0]) + (disp[1]*disp[1]) + (disp[2]*disp[2]);
4384 dist = VSQRT(dist) + atomRadius;
4387 soluteXlen = xmax -
xmin;
4388 soluteYlen = ymax -
ymin;
4389 soluteZlen = zmax -
zmin;
4392 Vnm_print(0,
"APOL: Setting up hash table and accessibility object...\n");
4393 nhash[0] = soluteXlen/0.5;
4394 nhash[1] = soluteYlen/0.5;
4395 nhash[2] = soluteZlen/0.5;
4396 for (i=0; i<3; i++) inhash[i] = (
int)(nhash[i]);
4399 if (inhash[i] < 3) inhash[i] = 3;
4400 if (inhash[i] > MAX_HASH_DIM) inhash[i] = MAX_HASH_DIM;
4404 srad = apolparm->
srad;
4405 sradPad = srad + (2*apolparm->
dpos);
4411 if (param == VNULL && (apolparm->
bconc != 0.0)) {
4412 Vnm_tprint(2,
"initAPOL: Got NULL Vparam object!\n");
4413 Vnm_tprint(2,
"initAPOL: You are performing an apolar calculation with the van der Waals integral term,\n");
4414 Vnm_tprint(2,
"initAPOL: this term requires van der Waals parameters which are not available from the \n");
4415 Vnm_tprint(2,
"initAPOL: PQR file. Therefore, you need to supply a parameter file with the parm keyword,\n");
4416 Vnm_tprint(2,
"initAPOL: for example,\n");
4417 Vnm_tprint(2,
"initAPOL: read parm flat amber94.dat end\n");
4418 Vnm_tprint(2,
"initAPOL: where the relevant parameter files can be found in apbs/tools/conversion/param/vparam.\n");
4422 if (apolparm->
bconc != 0.0){
4425 if (atomData == VNULL) {
4426 Vnm_tprint(2,
"initAPOL: Couldn't find parameters for WAT OW or WAT O!\n");
4427 Vnm_tprint(2,
"initAPOL: These parameters must be present in your file\n");
4428 Vnm_tprint(2,
"initAPOL: for apolar calculations.\n");
4438 rc =
forceAPOL(acc, mem, apolparm, nforce, atomForce, alist, clist);
4440 Vnm_print(2,
"Error in apolar force calculation!\n");
4452 if (VABS(apolparm->
gamma) > VSMALL) {
4456 for (i = 0; i < len; i++) {
4462 apolparm->
sasa = 0.0;
4464 for (i = 0; i < len; i++) {
4470 if (VABS(apolparm->
press) > VSMALL){
4473 apolparm->
sav = 0.0;
4477 if (VABS(apolparm->
bconc) > VSMALL) {
4479 for (i = 0; i < len; i++) {
4482 Vnm_print(2,
"Error in apolar energy calculation!\n");
4485 atomwcaEnergy[i] = energy;
4490 Vnm_print(2,
"Error in apolar energy calculation!\n");
4511 double atomwcaEnergy[],
4515 double energy = 0.0;
4519 Vnm_print(1,
"\nSolvent Accessible Surface Area (SASA) for each atom:\n");
4520 for (i = 0; i < numatoms; i++) {
4521 Vnm_print(1,
" SASA for atom %i: %1.12E\n", i, atomsasa[i]);
4524 Vnm_print(1,
"\nTotal solvent accessible surface area: %g A^2\n",sasa);
4531 Vnm_print(1,
"energyAPOL: Cannot calculate component energy, skipping.\n");
4534 energy = (apolparm->
gamma*sasa) + (apolparm->
press*sav)
4537 Vnm_print(1,
"\nSurface tension*area energies (gamma * SASA) for each atom:\n");
4538 for (i = 0; i < numatoms; i++) {
4539 Vnm_print(1,
" Surface tension*area energy for atom %i: %1.12E\n", i, apolparm->
gamma*atomsasa[i]);
4542 Vnm_print(1,
"\nTotal surface tension energy: %g kJ/mol\n", apolparm->
gamma*sasa);
4543 Vnm_print(1,
"\nTotal solvent accessible volume: %g A^3\n", sav);
4544 Vnm_print(1,
"\nTotal pressure*volume energy: %g kJ/mol\n", apolparm->
press*sav);
4545 Vnm_print(1,
"\nWCA dispersion Energies for each atom:\n");
4546 for (i = 0; i < numatoms; i++) {
4547 Vnm_print(1,
" WCA energy for atom %i: %1.12E\n", i, atomwcaEnergy[i]);
4550 Vnm_print(1,
"\nTotal WCA energy: %g kJ/mol\n", (apolparm->
wcaEnergy));
4551 Vnm_print(1,
"\nTotal non-polar energy = %1.12E kJ/mol\n", energy);
4555 Vnm_print(2,
"energyAPOL: Error in energyAPOL. Unknown option.\n");
4570 time_t ts, ts_main, ts_sub;
4588 Vatom *atom = VNULL;
4591 srad = apolparm->
srad;
4592 press = apolparm->
press;
4593 gamma = apolparm->
gamma;
4594 offset = apolparm->
dpos;
4595 bconc = apolparm->
bconc;
4600 Vnm_print(0,
"forceAPOL: Trying atom surf...\n");
4602 if (acc->
surf == VNULL) {
4604 for (i=0; i<natom; i++) {
4612 Vnm_print(0,
"forceAPOL: atom surf: Time elapsed: %f\n", ((
double)clock() - ts) / CLOCKS_PER_SEC);
4615 Vnm_print(0,
"forceAPOL: calcforce == ACF_TOTAL\n");
4619 if(*atomForce != VNULL){
4620 Vmem_free(mem,*nforce,
sizeof(
AtomForce), (
void **)atomForce);
4623 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4627 for (j=0; j<3; j++) {
4629 (*atomForce)[0].savForce[j] = 0.0;
4630 (*atomForce)[0].wcaForce[j] = 0.0;
4634 for (i=0; i<natom; i++) {
4643 if(VABS(gamma) > VSMALL) {
4646 if(VABS(press) > VSMALL) {
4649 if(VABS(bconc) > VSMALL) {
4654 (*atomForce)[0].sasaForce[j] += dSASA[j];
4655 (*atomForce)[0].savForce[j] += dSAV[j];
4656 (*atomForce)[0].wcaForce[j] += force[j];
4661 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A)\n");
4662 Vnm_tprint( 1,
" Legend:\n");
4663 Vnm_tprint( 1,
" sasa -- SASA force\n");
4664 Vnm_tprint( 1,
" sav -- SAV force\n");
4665 Vnm_tprint( 1,
" wca -- WCA force\n\n");
4667 Vnm_tprint( 1,
" sasa %4.3e %4.3e %4.3e\n",
4668 (*atomForce)[0].sasaForce[0],
4669 (*atomForce)[0].sasaForce[1],
4670 (*atomForce)[0].sasaForce[2]);
4671 Vnm_tprint( 1,
" sav %4.3e %4.3e %4.3e\n",
4672 (*atomForce)[0].savForce[0],
4673 (*atomForce)[0].savForce[1],
4674 (*atomForce)[0].savForce[2]);
4675 Vnm_tprint( 1,
" wca %4.3e %4.3e %4.3e\n",
4676 (*atomForce)[0].wcaForce[0],
4677 (*atomForce)[0].wcaForce[1],
4678 (*atomForce)[0].wcaForce[2]);
4680 Vnm_print(0,
"forceAPOL: calcforce == ACF_TOTAL: %f\n", ((
double)clock() - ts) / CLOCKS_PER_SEC);
4683 if(*atomForce == VNULL){
4684 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4687 Vmem_free(mem,*nforce,
sizeof(
AtomForce), (
void **)atomForce);
4688 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4693 Vnm_tprint( 1,
" Printing per atom forces (kJ/mol/A)\n");
4694 Vnm_tprint( 1,
" Legend:\n");
4695 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
4696 Vnm_tprint( 1,
" sasa n -- SASA force for atom n\n");
4697 Vnm_tprint( 1,
" sav n -- SAV force for atom n\n");
4698 Vnm_tprint( 1,
" wca n -- WCA force for atom n\n\n");
4700 Vnm_tprint( 1,
" gamma %f\n" \
4706 for (i=0; i<natom; i++) {
4716 for (j=0; j<3; j++) {
4717 (*atomForce)[i].sasaForce[j] = 0.0;
4718 (*atomForce)[i].savForce[j] = 0.0;
4719 (*atomForce)[i].wcaForce[j] = 0.0;
4722 if(VABS(gamma) > VSMALL)
Vacc_atomdSASA(acc, offset, srad, atom, dSASA);
4723 if(VABS(press) > VSMALL)
Vacc_atomdSAV(acc, srad, atom, dSAV);
4726 xF = -((gamma*dSASA[0]) + (press*dSAV[0]) + (bconc*force[0]));
4727 yF = -((gamma*dSASA[1]) + (press*dSAV[1]) + (bconc*force[1]));
4728 zF = -((gamma*dSASA[2]) + (press*dSAV[2]) + (bconc*force[2]));
4731 (*atomForce)[i].sasaForce[j] += dSASA[j];
4732 (*atomForce)[i].savForce[j] += dSAV[j];
4733 (*atomForce)[i].wcaForce[j] += force[j];
4737 Vnm_print( 1,
" tot %i %4.3e %4.3e %4.3e\n",
4742 Vnm_print( 1,
" sasa %i %4.3e %4.3e %4.3e\n",
4744 (*atomForce)[i].sasaForce[0],
4745 (*atomForce)[i].sasaForce[1],
4746 (*atomForce)[i].sasaForce[2]);
4747 Vnm_print( 1,
" sav %i %4.3e %4.3e %4.3e\n",
4749 (*atomForce)[i].savForce[0],
4750 (*atomForce)[i].savForce[1],
4751 (*atomForce)[i].savForce[2]);
4752 Vnm_print( 1,
" wca %i %4.3e %4.3e %4.3e\n",
4754 (*atomForce)[i].wcaForce[0],
4755 (*atomForce)[i].wcaForce[1],
4756 (*atomForce)[i].wcaForce[2]);
4766 Vnm_print(0,
"forceAPOL: Time elapsed: %f\n", ((
double)clock() - ts_main) / CLOCKS_PER_SEC);
VPUBLIC Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
VPUBLIC double Vpmg_qmEnergy(Vpmg *thee, int extFlag)
Get the "mobile charge" contribution to the electrostatic energy.
VPUBLIC double Vpmg_qfEnergy(Vpmg *thee, int extFlag)
Get the "fixed charge" contribution to the electrostatic energy.
Vdata_Format meshfmt[NOSH_MAXMOL]
VPUBLIC void Vgrid_writeDX(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out the data in OpenDX grid format.
int apol2calc[NOSH_MAXCALC]
int Vacc_wcaEnergyAtom(Vacc *thee, APOLparm *apolparm, Valist *alist, Vclist *clist, int iatom, double *value)
Calculate the WCA energy for an atom.
VPUBLIC Vparam * loadParameter(NOsh *nosh)
Loads and returns parameter object.
VPUBLIC void Vfetk_setAtomColors(Vfetk *thee)
Transfer color (partition ID) information frmo a partitioned mesh to the atoms.
VPUBLIC void Vgrid_writeUHBD(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out the data in UHBD grid format.
VPUBLIC Vfetk * Vfetk_ctor(Vpbe *pbe, Vhal_PBEType type)
Constructor for Vfetk object.
VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom)
Return the atomic solvent accessible surface area (SASA)
#define NOSH_MAXMOL
Maximum number of molecules in a run.
Contains public data members for Vpbe class/module.
Oracle for solvent- and ion-accessibility around a biomolecule.
Vdata_Format dielfmt[NOSH_MAXMOL]
VPUBLIC Vparam * Vparam_ctor()
Construct the object.
VPUBLIC int loadChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the charge maps given in NOsh into grid objects.
VPUBLIC void Vpmgp_dtor(Vpmgp **thee)
Object destructor.
VPUBLIC void Vgrid_dtor(Vgrid **thee)
Object destructor.
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
char potpath[NOSH_MAXMOL][VMAX_ARGLEN]
VPUBLIC int energyAPOL(APOLparm *apolparm, double sasa, double sav, double atomsasa[], double atomwcaEnergy[], int numatoms)
Calculate non-polar energies.
NOsh_PrintType printwhat[NOSH_MAXPRINT]
Electrostatic potential oracle for Cartesian mesh data.
VPUBLIC int Vgrid_readGZ(Vgrid *thee, const char *fname)
Read in OpenDX data in GZIP format.
VPUBLIC int writedataFlat(NOsh *nosh, Vcom *com, const char *fname, double totEnergy[NOSH_MAXCALC], double qfEnergy[NOSH_MAXCALC], double qmEnergy[NOSH_MAXCALC], double dielEnergy[NOSH_MAXCALC], int nenergy[NOSH_MAXCALC], double *atomEnergy[NOSH_MAXCALC], int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Write out information to a flat file.
Vdata_Format writefmt[PBEPARM_MAXWRITE]
char elecname[NOSH_MAXCALC][VMAX_ARGLEN]
VPUBLIC int postRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Estimate error, mark mesh, and refine mesh after solve.
NOsh_calc * calc[NOSH_MAXCALC]
Contains public data members for Vpmg class/module.
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
VPUBLIC int solveMG(NOsh *nosh, Vpmg *pmg, MGparm_CalcType type)
Solve the PBE with MG.
VPUBLIC void Vclist_dtor(Vclist **thee)
Destroy object.
VPUBLIC void Vpmg_setPart(Vpmg *thee, double lowerCorner[3], double upperCorner[3], int bflags[6])
Set partition information which restricts the calculation of observables to a (rectangular) subset of...
VPUBLIC Vpbe * Vpbe_ctor(Valist *alist, int ionNum, double *ionConc, double *ionRadii, double *ionQ, double T, double soluteDiel, double solventDiel, double solventRadius, int focusFlag, double sdens, double z_mem, double L, double membraneDiel, double V)
Construct Vpbe object.
VPUBLIC int writematMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out operator matrix from MG calculation to file.
VPUBLIC int setPartMG(NOsh *nosh, MGparm *mgparm, Vpmg *pmg)
Set MG partitions for calculating observables and performing I/O.
Header file for front end auxiliary routines.
VPUBLIC int Vpmg_dbForce(Vpmg *thee, double *dbForce, int atomID, Vsurf_Meth srfm)
Calculate the dielectric boundary forces on the specified atom in units of k_B T/AA.
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
VPUBLIC void killChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded charge maps.
#define APBS_TIMER_SETUP
APBS setup timer ID.
VPUBLIC Vclist * Vclist_ctor(Valist *alist, double max_radius, int npts[VAPBS_DIM], Vclist_DomainMode mode, double lower_corner[VAPBS_DIM], double upper_corner[VAPBS_DIM])
Construct the cell list object.
VPUBLIC void killForce(Vmem *mem, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Free memory from MG force calculation.
APOLparm_calcForce calcforce
char dielXpath[NOSH_MAXMOL][VMAX_ARGLEN]
VPUBLIC int loadKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the kappa maps given in NOsh into grid objects.
VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius)
Return the total solvent accessible surface area (SASA)
VPUBLIC int loadDielMaps(NOsh *nosh, Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL])
Load the dielectric maps given in NOsh into grid objects.
VPUBLIC int Vacc_wcaForceAtom(Vacc *thee, APOLparm *apolparm, Vclist *clist, Vatom *atom, double *force)
Return the WCA integral force.
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 int initMG(int icalc, NOsh *nosh, MGparm *mgparm, PBEparm *pbeparm, double realCenter[3], Vpbe *pbe[NOSH_MAXCALC], Valist *alist[NOSH_MAXMOL], Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL], Vgrid *kappaMap[NOSH_MAXMOL], Vgrid *chargeMap[NOSH_MAXMOL], Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC], Vgrid *potMap[NOSH_MAXMOL])
Initialize an MG calculation.
#define APBS_TIMER_SOLVER
APBS solver timer ID.
VPUBLIC void storeAtomEnergy(Vpmg *pmg, int icalc, double **atomEnergy, int *nenergy)
Store energy in arrays for future use.
VPUBLIC void killMG(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC])
Kill structures initialized during an MG calculation.
AtomData sub-class; stores atom data.
VPUBLIC int printApolForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data.
VPUBLIC void printFEPARM(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Print out FE-specific params loaded from input.
VPUBLIC void Vfetk_dtor(Vfetk **thee)
Object destructor.
VPUBLIC double Vpmg_qfAtomEnergy(Vpmg *thee, Vatom *atom)
Get the per-atom "fixed charge" contribution to the electrostatic energy.
VPUBLIC void killPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded potential maps.
VPUBLIC void Vfetk_setParameters(Vfetk *thee, PBEparm *pbeparm, FEMparm *feparm)
Set the parameter objects.
VPUBLIC void killFE(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vfetk *fetk[NOSH_MAXCALC], Gem *gm[NOSH_MAXMOL])
Kill structures initialized during an FE calculation.
Vdata_Format potfmt[NOSH_MAXMOL]
VPUBLIC Vrc_Codes Valist_readXML(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from an XML file.
VPUBLIC void Vpmg_dtor(Vpmg **thee)
Object destructor.
VPUBLIC int printElecEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data.
Vdata_Format kappafmt[NOSH_MAXMOL]
VPUBLIC int partFE(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Partition mesh (if applicable)
VPUBLIC void killKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded kappa maps.
APOLparm_calcEnergy calcenergy
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius)
Return the total solvent accessible volume (SAV)
enum eMGparm_CalcType MGparm_CalcType
Declare MGparm_CalcType type.
VPUBLIC void Vacc_atomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible area.
VPUBLIC int Vparam_readFlatFile(Vparam *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read a flat-file format parameter database.
VPUBLIC int Vpmg_fillArray(Vpmg *thee, double *vec, Vdata_Type type, double parm, Vhal_PBEType pbetype, PBEparm *pbeparm)
Fill the specified array with accessibility values.
VPUBLIC int Vstring_strcasecmp(const char *s1, const char *s2)
Case-insensitive string comparison (BSD standard)
VPUBLIC int energyFE(NOsh *nosh, int icalc, Vfetk *fetk[NOSH_MAXCALC], int *nenergy, double *totEnergy, double *qfEnergy, double *qmEnergy, double *dielEnergy)
Calculate electrostatic energies from FE solution.
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
VPUBLIC int printForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data (deprecated...see printElecForce)
VPUBLIC void printPBEPARM(PBEparm *pbeparm)
Print out generic PBE params loaded from input.
#define Vunit_Na
Avogadro's number.
VPUBLIC int writedataMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out observables from MG calculation to file.
VPUBLIC Vparam_AtomData * Vparam_getAtomData(Vparam *thee, char resName[VMAX_ARGLEN], char atomName[VMAX_ARGLEN])
Get atom data.
Parameter structure for FEM-specific variables from input files.
Structure to hold atomic forces.
Reads and assigns charge/radii parameters.
int printcalc[NOSH_MAXPRINT][NOSH_MAXPOP]
VPUBLIC double Vpmg_dielEnergy(Vpmg *thee, int extFlag)
Get the "polarization" contribution to the electrostatic energy.
VPUBLIC Vpmg * Vpmg_ctor(Vpmgp *pmgp, Vpbe *pbe, int focusFlag, Vpmg *pmgOLD, MGparm *mgparm, PBEparm_calcEnergy energyFlag)
Constructor for the Vpmg class (allocates new memory)
Contains public data members for Vpmgp class/module.
#define Vunit_kb
Boltzmann constant.
PBEparm_calcEnergy calcenergy
VPUBLIC int printEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data (deprecated...see printElecEnergy)
VPUBLIC double returnEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Access net local energy.
#define APBS_TIMER_ENERGY
APBS energy timer ID.
VPUBLIC int loadPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the potential maps given in NOsh into grid objects.
NOsh_MolFormat molfmt[NOSH_MAXMOL]
char chargepath[NOSH_MAXMOL][VMAX_ARGLEN]
Parameter structure for PBE variables from input files.
VPUBLIC int Vpmg_fillco(Vpmg *thee, Vsurf_Meth surfMeth, double splineWin, Vchrg_Meth chargeMeth, int useDielXMap, Vgrid *dielXMap, int useDielYMap, Vgrid *dielYMap, int useDielZMap, Vgrid *dielZMap, int useKappaMap, Vgrid *kappaMap, int usePotMap, Vgrid *potMap, int useChargeMap, Vgrid *chargeMap)
Fill the coefficient arrays prior to solving the equation.
VPUBLIC void killMolecules(NOsh *nosh, Valist *alist[NOSH_MAXMOL])
Destroy the loaded molecules.
PBEparm_calcForce calcforce
int printnarg[NOSH_MAXPRINT]
VPUBLIC int forceMG(Vmem *mem, NOsh *nosh, PBEparm *pbeparm, MGparm *mgparm, Vpmg *pmg, int *nforce, AtomForce **atomForce, Valist *alist[NOSH_MAXMOL])
Calculate forces from MG solution.
Surface object list of per-atom surface points.
VPUBLIC Vrc_Codes Valist_readPDB(Valist *thee, Vparam *param, Vio *sock)
Fill atom list with information from a PDB file.
char dielZpath[NOSH_MAXMOL][VMAX_ARGLEN]
FEMparm_EstType akeySOLVE
Parameter structure for MG-specific variables from input files.
Vdata_Format chargefmt[NOSH_MAXMOL]
VPUBLIC void Vpbe_dtor(Vpbe **thee)
Object destructor.
VPUBLIC void killDielMaps(NOsh *nosh, Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL])
Destroy the loaded dielectric.
VPUBLIC Vrc_Codes initFE(int icalc, NOsh *nosh, FEMparm *feparm, PBEparm *pbeparm, Vpbe *pbe[NOSH_MAXCALC], Valist *alist[NOSH_MAXMOL], Vfetk *fetk[NOSH_MAXCALC])
Initialize FE solver objects.
VPUBLIC int printElecForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data.
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
VPUBLIC void printMGPARM(MGparm *mgparm, double realCenter[3])
Print out MG-specific params loaded from input.
VPUBLIC int printApolEnergy(NOsh *nosh, int iprint)
Combine and pretty-print energy data.
Contains public data members for Vatom class/module.
char writematstem[VMAX_ARGLEN]
VPUBLIC int Vgrid_readDX(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read in data in OpenDX grid format.
Container class for list of atom objects.
VPUBLIC Vpmgp * Vpmgp_ctor(MGparm *mgparm)
Construct PMG parameter object and initialize to default values.
Class for parsing fixed format input files.
char meshpath[NOSH_MAXMOL][VMAX_ARGLEN]
VPUBLIC int Vparam_readXMLFile(Vparam *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read an XML format parameter database.
int elec2calc[NOSH_MAXCALC]
Contains public data members for Vfetk class/module.
VPUBLIC int Vpmg_ibForce(Vpmg *thee, double *force, int atomID, Vsurf_Meth srfm)
Calculate the osmotic pressure on the specified atom in units of k_B T/AA.
VPUBLIC int loadMolecules(NOsh *nosh, Vparam *param, Valist *alist[NOSH_MAXMOL])
Load the molecules given in NOsh into atom lists.
VPUBLIC void Vgrid_writeGZ(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out OpenDX data in GZIP format.
#define NOSH_MAXCALC
Maximum number of calculations in a run.
VPUBLIC int forceAPOL(Vacc *acc, Vmem *mem, APOLparm *apolparm, int *nforce, AtomForce **atomForce, Valist *alist, Vclist *clist)
Calculate non-polar forces.
Vdata_Type writetype[PBEPARM_MAXWRITE]
VPUBLIC int Vpmg_qfForce(Vpmg *thee, double *force, int atomID, Vchrg_Meth chgm)
Calculate the "charge-field" force on the specified atom in units of k_B T/AA.
VPUBLIC int solveFE(int icalc, PBEparm *pbeparm, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Solve-estimate-refine.
VPUBLIC Vrc_Codes Valist_readPQR(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from a PQR file.
VPUBLIC double Vfetk_energy(Vfetk *thee, int color, int nonlin)
Return the total electrostatic energy.
enum eVfetk_MeshLoad Vfetk_MeshLoad
Declare FEMparm_GuessType type.
VPUBLIC void startVio()
Wrapper to start MALOC Vio layer.
char parmpath[VMAX_ARGLEN]
VPUBLIC int Vpmg_solve(Vpmg *thee)
Solve the PBE using PMG.
VPUBLIC double Vpmg_energy(Vpmg *thee, int extFlag)
Get the total electrostatic energy.
VPUBLIC double Vpbe_getDeblen(Vpbe *thee)
Get Debye-Huckel screening length.
int printop[NOSH_MAXPRINT][NOSH_MAXPOP]
char molpath[NOSH_MAXMOL][VMAX_ARGLEN]
VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type)
Fill an array with the specified data.
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
char writestem[PBEPARM_MAXWRITE][VMAX_ARGLEN]
VPUBLIC void Valist_dtor(Valist **thee)
Destroys atom list object.
char apolname[NOSH_MAXCALC][VMAX_ARGLEN]
VPUBLIC void Vacc_atomdSAV(Vacc *thee, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible volume.
VPUBLIC int initAPOL(NOsh *nosh, Vmem *mem, Vparam *param, APOLparm *apolparm, int *nforce, AtomForce **atomForce, Valist *alist)
Upperlevel routine to the non-polar energy and force routines.
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 Valist * Valist_ctor()
Construct the atom list object.
#define APBS_TIMER_FORCE
APBS force timer ID.
Parameter structure for APOL-specific variables from input files.
VPUBLIC int writedataXML(NOsh *nosh, Vcom *com, const char *fname, double totEnergy[NOSH_MAXCALC], double qfEnergy[NOSH_MAXCALC], double qmEnergy[NOSH_MAXCALC], double dielEnergy[NOSH_MAXCALC], int nenergy[NOSH_MAXCALC], double *atomEnergy[NOSH_MAXCALC], int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Write out information to an XML file.
VPUBLIC VaccSurf * Vacc_atomSurf(Vacc *thee, Vatom *atom, VaccSurf *ref, double prad)
Set up an array of points corresponding to the SAS due to a particular atom.
VPUBLIC int energyMG(NOsh *nosh, int icalc, Vpmg *pmg, int *nenergy, double *totEnergy, double *qfEnergy, double *qmEnergy, double *dielEnergy)
Calculate electrostatic energies from MG solution.
VPUBLIC int preRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Pre-refine mesh before solve.
VPUBLIC int writedataFE(int rank, NOsh *nosh, PBEparm *pbeparm, Vfetk *fetk)
Write FEM data to files.
VPUBLIC Vgrid * Vgrid_ctor(int nx, int ny, int nz, double hx, double hy, double hzed, double xmin, double ymin, double zmin, double *data)
Construct Vgrid object with values obtained from Vpmg_readDX (for example)
VPUBLIC void killEnergy()
Kill arrays allocated for energies.
char dielYpath[NOSH_MAXMOL][VMAX_ARGLEN]
char kappapath[NOSH_MAXMOL][VMAX_ARGLEN]