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++)
253 double sum,hx,hy,hzed,xmin,ymin,zmin;
257 Vnm_tprint( 1,
"Got paths for %d dielectric map sets\n",nosh->
ndiel);
262 for (i=0; i<nosh->
ndiel; i++) {
263 Vnm_tprint( 1,
"Reading x-shifted dielectric map data from %s:\n", nosh->
dielXpath[i]);
264 dielXMap[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
272 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
278 nx = dielXMap[i]->nx;
279 ny = dielXMap[i]->ny;
280 nz = dielXMap[i]->nz;
283 hx = dielXMap[i]->hx;
284 hy = dielXMap[i]->hy;
285 hzed = dielXMap[i]->hzed;
288 xmin = dielXMap[i]->xmin;
289 ymin = dielXMap[i]->ymin;
290 zmin = dielXMap[i]->zmin;
291 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
292 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
293 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
296 for (ii=0; ii<(nx*ny*nz); ii++)
297 sum += (dielXMap[i]->data[ii]);
298 sum = sum*hx*hy*hzed;
299 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
307 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
313 nx = dielXMap[i]->nx;
314 ny = dielXMap[i]->ny;
315 nz = dielXMap[i]->nz;
318 hx = dielXMap[i]->hx;
319 hy = dielXMap[i]->hy;
320 hzed = dielXMap[i]->hzed;
323 xmin = dielXMap[i]->xmin;
324 ymin = dielXMap[i]->ymin;
325 zmin = dielXMap[i]->zmin;
326 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
327 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
328 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
331 for (ii=0; ii<(nx*ny*nz); ii++)
332 sum += (dielXMap[i]->data[ii]);
333 sum = sum*hx*hy*hzed;
334 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
340 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
346 nx = dielXMap[i]->nx;
347 ny = dielXMap[i]->ny;
348 nz = dielXMap[i]->nz;
351 hx = dielXMap[i]->hx;
352 hy = dielXMap[i]->hy;
353 hzed = dielXMap[i]->hzed;
356 xmin = dielXMap[i]->xmin;
357 ymin = dielXMap[i]->ymin;
358 zmin = dielXMap[i]->zmin;
359 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
360 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
361 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
364 for (ii=0; ii<(nx*ny*nz); ii++)
365 sum += (dielXMap[i]->data[ii]);
366 sum = sum*hx*hy*hzed;
367 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
371 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
375 Vnm_tprint( 2,
"AVS input not supported yet!\n");
379 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
382 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
386 Vnm_tprint( 1,
"Reading y-shifted dielectric map data from \
388 dielYMap[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
396 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
402 nx = dielYMap[i]->nx;
403 ny = dielYMap[i]->ny;
404 nz = dielYMap[i]->nz;
407 hx = dielYMap[i]->hx;
408 hy = dielYMap[i]->hy;
409 hzed = dielYMap[i]->hzed;
412 xmin = dielYMap[i]->xmin;
413 ymin = dielYMap[i]->ymin;
414 zmin = dielYMap[i]->zmin;
415 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
416 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
417 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
420 for (ii=0; ii<(nx*ny*nz); ii++)
421 sum += (dielYMap[i]->data[ii]);
422 sum = sum*hx*hy*hzed;
423 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
430 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
436 nx = dielYMap[i]->nx;
437 ny = dielYMap[i]->ny;
438 nz = dielYMap[i]->nz;
441 hx = dielYMap[i]->hx;
442 hy = dielYMap[i]->hy;
443 hzed = dielYMap[i]->hzed;
446 xmin = dielYMap[i]->xmin;
447 ymin = dielYMap[i]->ymin;
448 zmin = dielYMap[i]->zmin;
449 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
450 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
451 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
454 for (ii=0; ii<(nx*ny*nz); ii++)
455 sum += (dielYMap[i]->data[ii]);
456 sum = sum*hx*hy*hzed;
457 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
462 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
468 nx = dielYMap[i]->nx;
469 ny = dielYMap[i]->ny;
470 nz = dielYMap[i]->nz;
473 hx = dielYMap[i]->hx;
474 hy = dielYMap[i]->hy;
475 hzed = dielYMap[i]->hzed;
478 xmin = dielYMap[i]->xmin;
479 ymin = dielYMap[i]->ymin;
480 zmin = dielYMap[i]->zmin;
481 Vnm_tprint(1,
" %d x %d x %d grid\n", nx, ny, nz);
482 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n", hx, hy, hzed);
483 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
486 for (ii=0; ii<(nx*ny*nz); ii++)
487 sum += (dielYMap[i]->data[ii]);
488 sum = sum*hx*hy*hzed;
489 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
493 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
497 Vnm_tprint( 2,
"AVS input not supported yet!\n");
501 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
504 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
509 Vnm_tprint( 1,
"Reading z-shifted dielectric map data from \
511 dielZMap[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
519 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
525 nx = dielZMap[i]->nx;
526 ny = dielZMap[i]->ny;
527 nz = dielZMap[i]->nz;
530 hx = dielZMap[i]->hx;
531 hy = dielZMap[i]->hy;
532 hzed = dielZMap[i]->hzed;
535 xmin = dielZMap[i]->xmin;
536 ymin = dielZMap[i]->ymin;
537 zmin = dielZMap[i]->zmin;
538 Vnm_tprint(1,
" %d x %d x %d grid\n",
540 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
542 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
545 for (ii=0; ii<(nx*ny*nz); ii++) sum += (dielZMap[i]->data[ii]);
546 sum = sum*hx*hy*hzed;
547 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
554 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
560 nx = dielZMap[i]->nx;
561 ny = dielZMap[i]->ny;
562 nz = dielZMap[i]->nz;
565 hx = dielZMap[i]->hx;
566 hy = dielZMap[i]->hy;
567 hzed = dielZMap[i]->hzed;
570 xmin = dielZMap[i]->xmin;
571 ymin = dielZMap[i]->ymin;
572 zmin = dielZMap[i]->zmin;
573 Vnm_tprint(1,
" %d x %d x %d grid\n",
575 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
577 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
580 for (ii=0; ii<(nx*ny*nz); ii++) sum += (dielZMap[i]->data[ii]);
581 sum = sum*hx*hy*hzed;
582 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
587 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
593 nx = dielZMap[i]->nx;
594 ny = dielZMap[i]->ny;
595 nz = dielZMap[i]->nz;
598 hx = dielZMap[i]->hx;
599 hy = dielZMap[i]->hy;
600 hzed = dielZMap[i]->hzed;
603 xmin = dielZMap[i]->xmin;
604 ymin = dielZMap[i]->ymin;
605 zmin = dielZMap[i]->zmin;
606 Vnm_tprint(1,
" %d x %d x %d grid\n",
608 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
610 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
613 for (ii=0; ii<(nx*ny*nz); ii++) sum += (dielZMap[i]->data[ii]);
614 sum = sum*hx*hy*hzed;
615 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
619 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
623 Vnm_tprint( 2,
"AVS input not supported yet!\n");
627 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
630 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
646 if (nosh->
ndiel > 0) {
648 Vnm_tprint( 1,
"Destroying %d dielectric map sets\n",
651 for (i=0; i<nosh->
ndiel; i++) {
674 Vnm_tprint( 1,
"Got paths for %d kappa maps\n", nosh->
nkappa);
677 for (i=0; i<nosh->
nkappa; i++) {
678 Vnm_tprint( 1,
"Reading kappa map data from %s:\n",
680 map[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
688 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
692 Vnm_tprint(1,
" %d x %d x %d grid\n",
693 map[i]->nx, map[i]->ny, map[i]->nz);
694 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
695 map[i]->hx, map[i]->hy, map[i]->hzed);
696 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
697 map[i]->xmin, map[i]->ymin, map[i]->zmin);
699 for (ii = 0, len = map[i]->nx * map[i]->ny * map[i]->nz;
703 sum += (map[i]->data[ii]);
705 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
706 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
713 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
717 Vnm_tprint(1,
" %d x %d x %d grid\n",
718 map[i]->nx, map[i]->ny, map[i]->nz);
719 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
720 map[i]->hx, map[i]->hy, map[i]->hzed);
721 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
722 map[i]->xmin, map[i]->ymin, map[i]->zmin);
724 for (ii = 0, len = map[i]->nx * map[i]->ny * map[i]->nz;
728 sum += (map[i]->data[ii]);
730 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
731 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
735 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
739 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
743 Vnm_tprint( 2,
"AVS input not supported yet!\n");
748 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
752 Vnm_tprint(1,
" %d x %d x %d grid\n",
753 map[i]->nx, map[i]->ny, map[i]->nz);
754 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
755 map[i]->hx, map[i]->hy, map[i]->hzed);
756 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
757 map[i]->xmin, map[i]->ymin, map[i]->zmin);
759 for (ii=0, len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
760 sum += (map[i]->data[ii]);
762 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
763 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
766 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
782 Vnm_tprint( 1,
"Destroying %d kappa maps\n", nosh->
nkappa);
803 Vnm_tprint( 1,
"Got paths for %d potential maps\n", nosh->
npot);
806 for (i=0; i<nosh->
npot; i++) {
807 Vnm_tprint( 1,
"Reading potential map data from %s:\n",
809 map[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
810 switch (nosh->
potfmt[i]) {
818 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
824 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
829 Vnm_tprint(1,
" %d x %d x %d grid\n",
830 map[i]->nx, map[i]->ny, map[i]->nz);
831 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
832 map[i]->hx, map[i]->hy, map[i]->hzed);
833 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
834 map[i]->xmin, map[i]->ymin, map[i]->zmin);
836 for (ii=0,len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
837 sum += (map[i]->data[ii]);
839 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
840 Vnm_tprint(1,
" Volume integral = %3.2e A^3\n", sum);
844 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
848 Vnm_tprint( 2,
"MCSF input not supported yet!\n");
852 Vnm_tprint( 2,
"AVS input not supported yet!\n");
855 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
871 if (nosh->
npot > 0) {
873 Vnm_tprint( 1,
"Destroying %d potential maps\n", nosh->
npot);
894 Vnm_tprint( 1,
"Got paths for %d charge maps\n", nosh->
ncharge);
897 for (i=0; i<nosh->
ncharge; i++) {
898 Vnm_tprint( 1,
"Reading charge map data from %s:\n",
900 map[i] =
Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
907 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
911 Vnm_tprint(1,
" %d x %d x %d grid\n",
912 map[i]->nx, map[i]->ny, map[i]->nz);
913 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
914 map[i]->hx, map[i]->hy, map[i]->hzed);
915 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
916 map[i]->xmin, map[i]->ymin, map[i]->zmin);
918 for (ii=0,len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
919 sum += (map[i]->data[ii]);
921 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
922 Vnm_tprint(1,
" Charge map integral = %3.2e e\n", sum);
928 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
932 Vnm_tprint(1,
" %d x %d x %d grid\n",
933 map[i]->nx, map[i]->ny, map[i]->nz);
934 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
935 map[i]->hx, map[i]->hy, map[i]->hzed);
936 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
937 map[i]->xmin, map[i]->ymin, map[i]->zmin);
939 for (ii=0,len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
940 sum += (map[i]->data[ii]);
942 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
943 Vnm_tprint(1,
" Charge map integral = %3.2e e\n", sum);
946 Vnm_tprint( 2,
"UHBD input not supported yet!\n");
949 Vnm_tprint( 2,
"AVS input not supported yet!\n");
952 Vnm_tprint(2,
"MCSF input not supported yet!\n");
956 Vnm_tprint( 2,
"Fatal error while reading from %s\n",
960 Vnm_tprint(1,
" %d x %d x %d grid\n",
961 map[i]->nx, map[i]->ny, map[i]->nz);
962 Vnm_tprint(1,
" (%g, %g, %g) A spacings\n",
963 map[i]->hx, map[i]->hy, map[i]->hzed);
964 Vnm_tprint(1,
" (%g, %g, %g) A lower corner\n",
965 map[i]->xmin, map[i]->ymin, map[i]->zmin);
967 for (ii=0,len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
968 sum += (map[i]->data[ii]);
970 sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
971 Vnm_tprint(1,
" Charge map integral = %3.2e e\n", sum);
974 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
992 Vnm_tprint( 1,
"Destroying %d charge maps\n", nosh->
ncharge);
1005 double ionstr = 0.0;
1007 for (i=0; i<pbeparm->
nion; i++)
1008 ionstr += 0.5*(VSQR(pbeparm->
ionq[i])*pbeparm->
ionc[i]);
1010 Vnm_tprint( 1,
" Molecule ID: %d\n", pbeparm->
molid);
1013 Vnm_tprint( 1,
" Nonlinear traditional PBE\n");
1016 Vnm_tprint( 1,
" Linearized traditional PBE\n");
1019 Vnm_tprint( 1,
" Nonlinear regularized PBE\n");
1020 Vnm_tprint( 2,
" ** Sorry, but Nathan broke the nonlinear regularized PBE implementation. **\n");
1021 Vnm_tprint( 2,
" ** Please let us know if you are interested in using it. **\n");
1025 Vnm_tprint( 1,
" Linearized regularized PBE\n");
1028 Vnm_tprint( 1,
" Nonlinear Size-Modified PBE\n");
1031 Vnm_tprint(2,
" Unknown PBE type (%d)!\n", pbeparm->
pbetype);
1035 Vnm_tprint( 1,
" Zero boundary conditions\n");
1037 Vnm_tprint( 1,
" Single Debye-Huckel sphere boundary \
1040 Vnm_tprint( 1,
" Multiple Debye-Huckel sphere boundary \
1043 Vnm_tprint( 1,
" Boundary conditions from focusing\n");
1045 Vnm_tprint( 1,
" Boundary conditions from potential map\n");
1047 Vnm_tprint( 1,
" Membrane potential boundary conditions.\n");
1049 Vnm_tprint( 1,
" %d ion species (%4.3f M ionic strength):\n",
1050 pbeparm->
nion, ionstr);
1051 for (i=0; i<pbeparm->
nion; i++) {
1052 Vnm_tprint( 1,
" %4.3f A-radius, %4.3f e-charge, \
1053%4.3f M concentration\n",
1058 Vnm_tprint( 1,
" Lattice spacing: %4.3f A (SMPBE) \n", pbeparm->
smvolume);
1059 Vnm_tprint( 1,
" Relative size parameter: %4.3f (SMPBE) \n", pbeparm->
smsize);
1062 Vnm_tprint( 1,
" Solute dielectric: %4.3f\n", pbeparm->
pdie);
1063 Vnm_tprint( 1,
" Solvent dielectric: %4.3f\n", pbeparm->
sdie);
1064 switch (pbeparm->
srfm) {
1066 Vnm_tprint( 1,
" Using \"molecular\" surface \
1067definition; no smoothing\n");
1068 Vnm_tprint( 1,
" Solvent probe radius: %4.3f A\n",
1072 Vnm_tprint( 1,
" Using \"molecular\" surface definition;\
1073harmonic average smoothing\n");
1074 Vnm_tprint( 1,
" Solvent probe radius: %4.3f A\n",
1078 Vnm_tprint( 1,
" Using spline-based surface definition;\
1079window = %4.3f\n", pbeparm->
swin);
1084 Vnm_tprint( 1,
" Temperature: %4.3f K\n", pbeparm->
temp);
1086energies will be calculated\n");
1088forces will be calculated \n");
1090solvent forces will be calculated\n");
1091 for (i=0; i<pbeparm->
numwrite; i++) {
1094 Vnm_tprint(1,
" Charge distribution to be written to ");
1097 Vnm_tprint(1,
" Potential to be written to ");
1100 Vnm_tprint(1,
" Molecular solvent accessibility \
1104 Vnm_tprint(1,
" Spline-based solvent accessibility \
1108 Vnm_tprint(1,
" van der Waals solvent accessibility \
1112 Vnm_tprint(1,
" Ion accessibility to be written to ");
1115 Vnm_tprint(1,
" Potential Laplacian to be written to ");
1118 Vnm_tprint(1,
" Energy density to be written to ");
1121 Vnm_tprint(1,
" Ion number density to be written to ");
1124 Vnm_tprint(1,
" Ion charge density to be written to ");
1127 Vnm_tprint(1,
" X-shifted dielectric map to be written \
1131 Vnm_tprint(1,
" Y-shifted dielectric map to be written \
1135 Vnm_tprint(1,
" Z-shifted dielectric map to be written \
1139 Vnm_tprint(1,
" Kappa map to be written to ");
1142 Vnm_tprint(1,
" Atom potentials to be written to ");
1145 Vnm_tprint(2,
" Invalid data type for writing!\n");
1150 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"dx");
1153 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"dxbin");
1156 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"dx.gz");
1159 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"grd");
1162 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"ucd");
1165 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"mcsf");
1168 Vnm_tprint(1,
"%s.%s\n", pbeparm->
writestem[i],
"txt");
1171 Vnm_tprint(2,
" Invalid format for writing!\n");
1181 switch (mgparm->
chgm) {
1183 Vnm_tprint(1,
" Using linear spline charge discretization.\n");
1186 Vnm_tprint(1,
" Using cubic spline charge discretization.\n");
1192 Vnm_tprint( 1,
" Partition overlap fraction = %g\n",
1194 Vnm_tprint( 1,
" Processor array = %d x %d x %d\n",
1197 Vnm_tprint( 1,
" Grid dimensions: %d x %d x %d\n",
1199 Vnm_tprint( 1,
" Grid spacings: %4.3f x %4.3f x %4.3f\n",
1201 Vnm_tprint( 1,
" Grid lengths: %4.3f x %4.3f x %4.3f\n",
1203 Vnm_tprint( 1,
" Grid center: (%4.3f, %4.3f, %4.3f)\n",
1204 realCenter[0], realCenter[1], realCenter[2]);
1205 Vnm_tprint( 1,
" Multigrid levels: %d\n", mgparm->
nlev);
1215 double realCenter[3],
1236 Vatom *atom = VNULL;
1237 Vgrid *theDielXMap = VNULL,
1238 *theDielYMap = VNULL,
1239 *theDielZMap = VNULL;
1240 Vgrid *theKappaMap = VNULL,
1242 *theChargeMap = VNULL;
1248 for (j=0; j<3; j++) realCenter[j] = mgparm->
center[j];
1252 myalist = alist[pbeparm->
molid-1];
1266 Vnm_tprint(0,
"Setting up PBE object...\n");
1268 sparm = pbeparm->
swin;
1270 sparm = pbeparm->
srad;
1272 if (pbeparm->
nion > 0) {
1273 iparm = pbeparm->
ionr[0];
1279 Vnm_tprint( 2,
"Can't focus first calculation!\n");
1291 pbeparm->
sdie, sparm, focusFlag, pbeparm->
sdens,
1296 Vnm_tprint(0,
"Setting up PDE object...\n");
1301 mgparm->
method = (mgparm->
useAqua == 1) ? VSOL_NewtonAqua : VSOL_Newton;
1307 mgparm->
method = (mgparm->
useAqua == 1) ? VSOL_CGMGAqua : VSOL_MG;
1311 Vnm_tprint(2,
"Sorry, LRPBE isn't supported with the MG solver!\n");
1314 Vnm_tprint(2,
"Sorry, NRPBE isn't supported with the MG solver!\n");
1318 Vnm_tprint(2,
" ** Sorry, due to numerical stability issues SMPBE is currently disabled. We apologize for the inconvenience.\n");
1319 Vnm_tprint(2,
" ** Please let us know if you would like to use it in the future.\n");
1335 Vnm_tprint(2,
"Error! Unknown PBE type (%d)!\n", pbeparm->
pbetype);
1338 Vnm_tprint(0,
"Setting PDE center to local center...\n");
1339 pmgp[icalc]->bcfl = pbeparm->
bcfl;
1340 pmgp[icalc]->xcent = realCenter[0];
1341 pmgp[icalc]->ycent = realCenter[1];
1342 pmgp[icalc]->zcent = realCenter[2];
1346 Vnm_tprint( 2,
"Can't focus first calculation!\n");
1351 pmg[icalc] =
Vpmg_ctor(pmgp[icalc], pbe[icalc], 1, pmg[icalc-1],
1357 if (icalc>0)
Vpmg_dtor(&(pmg[icalc-1]));
1358 pmg[icalc] =
Vpmg_ctor(pmgp[icalc], pbe[icalc], 0, VNULL, mgparm,
PCE_NO);
1366 theDielXMap = dielXMap[pbeparm->
dielMapID-1];
1368 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1375 theDielYMap = dielYMap[pbeparm->
dielMapID-1];
1377 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1384 theDielZMap = dielZMap[pbeparm->
dielMapID-1];
1386 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1393 theKappaMap = kappaMap[pbeparm->
kappaMapID-1];
1395 Vnm_print(2,
"Error! %d is not a valid kappa map ID!\n",
1402 thePotMap = potMap[pbeparm->
potMapID-1];
1404 Vnm_print(2,
"Error! %d is not a valid potential map ID!\n",
1413 Vnm_print(2,
"Error! %d is not a valid charge map ID!\n",
1420 Vnm_print(2,
"Warning: You specified 'bcfl map' in the input file, but no potential map was found.\n");
1421 Vnm_print(2,
" You must specify 'usemap pot' statement in the APBS input file!\n");
1422 Vnm_print(2,
"Bailing out ...\n");
1435 Vnm_print(2,
"initMG: problems setting up coefficients (fillco)!\n");
1441 Vnm_tprint(1,
" Debye length: %g A\n",
Vpbe_getDeblen(pbe[icalc]));
1448 bytesTotal = Vmem_bytesTotal();
1449 highWater = Vmem_highWaterTotal();
1452 Vnm_tprint( 1,
" Current memory usage: %4.3f MB total, \
1453%4.3f MB high water\n", (
double)(bytesTotal)/(1024.*1024.),
1454 (
double)(highWater)/(1024.*1024.));
1467 Vnm_tprint(1,
"Destroying multigrid structures.\n");
1480 for(i=0;i<nosh->
ncalc;i++){
1497 if (nosh != VNULL) {
1498 if (nosh->
bogus)
return 1;
1506 Vnm_print(2,
" Error during PDE solution!\n");
1510 Vnm_tprint( 1,
" Skipping solve for mg-dummy run; zeroing \
1515 for (i=0; i<nx*ny*nz; i++) pmg->
u[i] = 0.0;
1532 if (nosh->
bogus)
return 1;
1535 for (j=0; j<3; j++) {
1540 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part center = (%g, %g, %g)\n",
1546 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part lower corner = (%g, %g, %g)\n",
1547 __FILE__, __LINE__, partMin[0], partMin[1], partMin[2]);
1548 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part upper corner = (%g, %g, %g)\n",
1550 partMax[0], partMax[1], partMax[2]);
1553 for (j=0; j<3; j++) {
1554 partMin[j] = mgparm->
center[j] - 0.5*mgparm->
glen[j];
1555 partMax[j] = mgparm->
center[j] + 0.5*mgparm->
glen[j];
1596 if (nosh->
bogus == 0) {
1599 Vnm_tprint( 1,
" Total electrostatic energy = %1.12E kJ/mol\n",
1602 }
else *totEnergy = 0;
1610 Vnm_tprint( 1,
" Total electrostatic energy = %1.12E \
1612 Vnm_tprint( 1,
" Fixed charge energy = %g kJ/mol\n",
1614 Vnm_tprint( 1,
" Mobile charge energy = %g kJ/mol\n",
1616 Vnm_tprint( 1,
" Dielectric energy = %g kJ/mol\n",
1618 Vnm_tprint( 1,
" Per-atom energies:\n");
1625 Vnm_tprint( 1,
" Atom %d: %1.12E kJ/mol\n", i,
1629 }
else *nenergy = 0;
1655 Vnm_tprint( 1,
" Calculating forces...\n");
1662 for (j=0; j<3; j++) {
1663 (*atomForce)[0].
qfForce[j] = 0;
1664 (*atomForce)[0].ibForce[j] = 0;
1665 (*atomForce)[0].dbForce[j] = 0;
1668 if (nosh->
bogus == 0) {
1673 for (k=0; k<3; k++) {
1679 for (k=0; k<3; k++) {
1680 (*atomForce)[0].qfForce[k] += qfForce[k];
1681 (*atomForce)[0].ibForce[k] += ibForce[k];
1682 (*atomForce)[0].dbForce[k] += dbForce[k];
1686 Vnm_tprint( 1,
" Printing net forces for molecule %d (kJ/mol/A)\n",
1688 Vnm_tprint( 1,
" Legend:\n");
1689 Vnm_tprint( 1,
" qf -- fixed charge force\n");
1690 Vnm_tprint( 1,
" db -- dielectric boundary force\n");
1691 Vnm_tprint( 1,
" ib -- ionic boundary force\n");
1692 Vnm_tprint( 1,
" qf %4.3e %4.3e %4.3e\n",
1696 Vnm_tprint( 1,
" ib %4.3e %4.3e %4.3e\n",
1700 Vnm_tprint( 1,
" db %4.3e %4.3e %4.3e\n",
1707 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
1710 Vnm_tprint( 1,
" Printing per-atom forces for molecule %d (kJ/mol/A)\n",
1712 Vnm_tprint( 1,
" Legend:\n");
1713 Vnm_tprint( 1,
" tot n -- total force for atom n\n");
1714 Vnm_tprint( 1,
" qf n -- fixed charge force for atom n\n");
1715 Vnm_tprint( 1,
" db n -- dielectric boundary force for atom n\n");
1716 Vnm_tprint( 1,
" ib n -- ionic boundary force for atom n\n");
1719 if (nosh->
bogus == 0) {
1727 for (k=0; k<3; k++) {
1728 (*atomForce)[j].qfForce[k] = 0;
1729 (*atomForce)[j].ibForce[k] = 0;
1730 (*atomForce)[j].dbForce[k] = 0;
1734 Vnm_tprint( 1,
"mgF tot %d %4.3e %4.3e %4.3e\n", j,
1736 *((*atomForce)[j].qfForce[0]+(*atomForce)[j].ibForce[0]+
1737 (*atomForce)[j].dbForce[0]),
1739 *((*atomForce)[j].qfForce[1]+(*atomForce)[j].ibForce[1]+
1740 (*atomForce)[j].dbForce[1]),
1742 *((*atomForce)[j].qfForce[2]+(*atomForce)[j].ibForce[2]+
1743 (*atomForce)[j].dbForce[2]));
1744 Vnm_tprint( 1,
"mgF qf %d %4.3e %4.3e %4.3e\n", j,
1746 *(*atomForce)[j].qfForce[0],
1748 *(*atomForce)[j].qfForce[1],
1750 *(*atomForce)[j].qfForce[2]);
1751 Vnm_tprint( 1,
"mgF ib %d %4.3e %4.3e %4.3e\n", j,
1753 *(*atomForce)[j].ibForce[0],
1755 *(*atomForce)[j].ibForce[1],
1757 *(*atomForce)[j].ibForce[2]);
1758 Vnm_tprint( 1,
"mgF db %d %4.3e %4.3e %4.3e\n", j,
1760 *(*atomForce)[j].dbForce[0],
1762 *(*atomForce)[j].dbForce[1],
1764 *(*atomForce)[j].dbForce[2]);
1777 Vnm_tprint(1,
"No energy arrays to destroy.\n");
1788 Vnm_tprint(1,
"Destroying force arrays.\n");
1791 for (i=0; i<nosh->
ncalc; i++) {
1793 if (nforce[i] > 0) Vmem_free(mem, nforce[i],
sizeof(
AtomForce),
1794 (
void **)&(atomForce[i]));
1801 char writematstem[VMAX_ARGLEN];
1802 char outpath[VMAX_ARGLEN];
1806 if (nosh->
bogus)
return 1;
1809 strlenmax = VMAX_ARGLEN-14;
1811 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1813 Vnm_tprint(2,
" Not writing matrix!\n");
1816 sprintf(writematstem,
"%s-PE%d", pbeparm->
writematstem, rank);
1818 strlenmax = (int)(VMAX_ARGLEN)-1;
1820 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1822 Vnm_tprint(2,
" Not writing matrix!\n");
1833 strlenmax = VMAX_ARGLEN-5;
1835 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1837 Vnm_tprint(2,
" Not writing matrix!\n");
1840 sprintf(outpath,
"%s.%s", writematstem,
"mat");
1846 Vnm_tprint( 1,
" Writing Poisson operator matrix \
1847to %s...\n", outpath);
1851 Vnm_tprint( 1,
" Writing linearization of full \
1852Poisson-Boltzmann operator matrix to %s...\n", outpath);
1855 Vnm_tprint( 2,
" Bogus matrix specification\
1860 Vnm_tprint(0,
" Printing operator...\n");
1879 *atomEnergy = (
double *)Vmem_malloc(pmg->
vmem, *nenergy,
sizeof(
double));
1881 for (i=0; i<*nenergy; i++) {
1902 int ielec, icalc, i, j;
1903 char *timestring = VNULL;
1906 double conversion, ltenergy, gtenergy, scalar;
1908 if (nosh->
bogus)
return 1;
1914 file = fopen(fname,
"w");
1915 if (file == VNULL) {
1916 Vnm_print(2,
"writedataFlat: Problem opening virtual socket %s\n",
1924 timestring = ctime(&now);
1925 fprintf(file,
"%s\n", timestring);
1927 for (ielec=0; ielec<nosh->
nelec;ielec++) {
1937 fprintf(file,
"elec");
1939 fprintf(file,
" name %s\n", nosh->
elecname[ielec]);
1940 }
else fprintf(file,
"\n");
1942 switch (mgparm->
type) {
1944 fprintf(file,
" mg-dummy\n");
1947 fprintf(file,
" mg-manual\n");
1950 fprintf(file,
" mg-auto\n");
1953 fprintf(file,
" mg-para\n");
1959 fprintf(file,
" mol %d\n", pbeparm->
molid);
1960 fprintf(file,
" dime %d %d %d\n", mgparm->
dime[0], mgparm->
dime[1],\
1965 fprintf(file,
" npbe\n");
1968 fprintf(file,
" lpbe\n");
1974 if (pbeparm->
nion > 0) {
1975 for (i=0; i<pbeparm->
nion; i++) {
1976 fprintf(file,
" ion %4.3f %4.3f %4.3f\n",
1981 fprintf(file,
" pdie %4.3f\n", pbeparm->
pdie);
1982 fprintf(file,
" sdie %4.3f\n", pbeparm->
sdie);
1984 switch (pbeparm->
srfm) {
1986 fprintf(file,
" srfm mol\n");
1987 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1990 fprintf(file,
" srfm smol\n");
1991 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1994 fprintf(file,
" srfm spl2\n");
1995 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
2001 switch (pbeparm->
bcfl) {
2003 fprintf(file,
" bcfl zero\n");
2006 fprintf(file,
" bcfl sdh\n");
2009 fprintf(file,
" bcfl mdh\n");
2012 fprintf(file,
" bcfl focus\n");
2015 fprintf(file,
" bcfl map\n");
2018 fprintf(file,
" bcfl mem\n");
2024 fprintf(file,
" temp %4.3f\n", pbeparm->
temp);
2026 for (;icalc<=nosh->
elec2calc[ielec];icalc++){
2032 fprintf(file,
" calc\n");
2033 fprintf(file,
" id %i\n", (icalc+1));
2034 fprintf(file,
" grid %4.3f %4.3f %4.3f\n",
2036 fprintf(file,
" glen %4.3f %4.3f %4.3f\n",
2040 fprintf(file,
" totEnergy %1.12E kJ/mol\n",
2041 (totEnergy[icalc]*conversion));
2043 fprintf(file,
" totEnergy %1.12E kJ/mol\n",
2044 (totEnergy[icalc]*conversion));
2045 fprintf(file,
" qfEnergy %1.12E kJ/mol\n",
2046 (0.5*qfEnergy[icalc]*conversion));
2047 fprintf(file,
" qmEnergy %1.12E kJ/mol\n",
2048 (qmEnergy[icalc]*conversion));
2049 fprintf(file,
" dielEnergy %1.12E kJ/mol\n",
2050 (dielEnergy[icalc]*conversion));
2051 for (i=0; i<nenergy[icalc]; i++){
2052 fprintf(file,
" atom %i %1.12E kJ/mol\n", i,
2053 (0.5*atomEnergy[icalc][i]*conversion));
2059 fprintf(file,
" qfForce %1.12E %1.12E %1.12E kJ/mol/A\n",
2060 (atomForce[icalc][0].qfForce[0]*conversion),
2061 (atomForce[icalc][0].qfForce[1]*conversion),
2062 (atomForce[icalc][0].qfForce[2]*conversion));
2063 fprintf(file,
" ibForce %1.12E %1.12E %1.12E kJ/mol/A\n",
2064 (atomForce[icalc][0].ibForce[0]*conversion),
2065 (atomForce[icalc][0].ibForce[1]*conversion),
2066 (atomForce[icalc][0].ibForce[2]*conversion));
2067 fprintf(file,
" dbForce %1.12E %1.12E %1.12E kJ/mol/A\n",
2068 (atomForce[icalc][0].dbForce[0]*conversion),
2069 (atomForce[icalc][0].dbForce[1]*conversion),
2070 (atomForce[icalc][0].dbForce[2]*conversion));
2072 fprintf(file,
" end\n");
2075 fprintf(file,
"end\n");
2080for (i=0; i<nosh->
nprint; i++) {
2084 fprintf(file,
"print energy");
2085 fprintf(file,
" %d", nosh->
printcalc[i][0]+1);
2088 if (nosh->
printop[i][j-1] == 0) fprintf(file,
" +");
2089 else if (nosh->
printop[i][j-1] == 1) fprintf(file,
" -");
2090 fprintf(file,
" %d", nosh->
printcalc[i][j]+1);
2093 fprintf(file,
"\n");
2102 if (nosh->
printop[i][j-1] == 0) scalar = 1.0;
2103 else if (nosh->
printop[i][j-1] == 1) scalar = -1.0;
2108 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2111 fprintf(file,
" localEnergy %1.12E kJ/mol\n", \
2113 fprintf(file,
" globalEnergy %1.12E kJ/mol\nend\n", \
2135 int ielec, icalc, i, j;
2136 char *timestring = VNULL;
2140 double conversion, ltenergy, gtenergy, scalar;
2142 if (nosh->
bogus)
return 1;
2148 file = fopen(fname,
"w");
2149 if (file == VNULL) {
2150 Vnm_print(2,
"writedataXML: Problem opening virtual socket %s\n",
2155 fprintf(file,
"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
2156 fprintf(file,
"<APBS>\n");
2161 timestring = ctime(&now);
2162 for(c = timestring; *c !=
'\n'; c++);
2164 fprintf(file,
" <date>%s</date>\n", timestring);
2166 for (ielec=0; ielec<nosh->
nelec;ielec++){
2176 fprintf(file,
" <elec>\n");
2178 fprintf(file,
" <name>%s</name>\n", nosh->
elecname[ielec]);
2181 switch (mgparm->
type) {
2183 fprintf(file,
" <type>mg-dummy</type>\n");
2186 fprintf(file,
" <type>mg-manual</type>\n");
2189 fprintf(file,
" <type>mg-auto</type>\n");
2192 fprintf(file,
" <type>mg-para</type>\n");
2198 fprintf(file,
" <molid>%d</molid>\n", pbeparm->
molid);
2199 fprintf(file,
" <nx>%d</nx>\n", mgparm->
dime[0]);
2200 fprintf(file,
" <ny>%d</ny>\n", mgparm->
dime[1]);
2201 fprintf(file,
" <nz>%d</nz>\n", mgparm->
dime[2]);
2205 fprintf(file,
" <pbe>npbe</pbe>\n");
2208 fprintf(file,
" <pbe>lpbe</pbe>\n");
2214 if (pbeparm->
nion > 0) {
2215 for (i=0; i<pbeparm->
nion; i++) {
2216 fprintf(file,
" <ion>\n");
2217 fprintf(file,
" <radius>%4.3f A</radius>\n",
2219 fprintf(file,
" <charge>%4.3f A</charge>\n",
2221 fprintf(file,
" <concentration>%4.3f M</concentration>\n",
2223 fprintf(file,
" </ion>\n");
2228 fprintf(file,
" <pdie>%4.3f</pdie>\n", pbeparm->
pdie);
2229 fprintf(file,
" <sdie>%4.3f</sdie>\n", pbeparm->
sdie);
2231 switch (pbeparm->
srfm) {
2233 fprintf(file,
" <srfm>mol</srfm>\n");
2234 fprintf(file,
" <srad>%4.3f</srad>\n", pbeparm->
srad);
2237 fprintf(file,
" <srfm>smol</srfm>\n");
2238 fprintf(file,
" <srad>%4.3f</srad>\n", pbeparm->
srad);
2241 fprintf(file,
" <srfm>spl2</srfm>\n");
2247 switch (pbeparm->
bcfl) {
2249 fprintf(file,
" <bcfl>zero</bcfl>\n");
2252 fprintf(file,
" <bcfl>sdh</bcfl>\n");
2255 fprintf(file,
" <bcfl>mdh</bcfl>\n");
2258 fprintf(file,
" <bcfl>focus</bcfl>\n");
2261 fprintf(file,
" <bcfl>map</bcfl>\n");
2264 fprintf(file,
" <bcfl>mem</bcfl>\n");
2270 fprintf(file,
" <temp>%4.3f K</temp>\n", pbeparm->
temp);
2272 for (;icalc<=nosh->
elec2calc[ielec];icalc++){
2278 fprintf(file,
" <calc>\n");
2279 fprintf(file,
" <id>%i</id>\n", (icalc+1));
2280 fprintf(file,
" <hx>%4.3f A</hx>\n", mgparm->
grid[0]);
2281 fprintf(file,
" <hy>%4.3f A</hy>\n", mgparm->
grid[1]);
2282 fprintf(file,
" <hz>%4.3f A</hz>\n", mgparm->
grid[2]);
2283 fprintf(file,
" <xlen>%4.3f A</xlen>\n", mgparm->
glen[0]);
2284 fprintf(file,
" <ylen>%4.3f A</ylen>\n", mgparm->
glen[1]);
2285 fprintf(file,
" <zlen>%4.3f A</zlen>\n", mgparm->
glen[2]);
2288 fprintf(file,
" <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2289 (totEnergy[icalc]*conversion));
2291 fprintf(file,
" <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2292 (totEnergy[icalc]*conversion));
2293 fprintf(file,
" <qfEnergy>%1.12E kJ/mol</qfEnergy>\n",
2294 (0.5*qfEnergy[icalc]*conversion));
2295 fprintf(file,
" <qmEnergy>%1.12E kJ/mol</qmEnergy>\n",
2296 (qmEnergy[icalc]*conversion));
2297 fprintf(file,
" <dielEnergy>%1.12E kJ/mol</dielEnergy>\n",
2298 (dielEnergy[icalc]*conversion));
2299 for (i=0; i<nenergy[icalc]; i++){
2300 fprintf(file,
" <atom>\n");
2301 fprintf(file,
" <id>%i</id>\n", i+1);
2302 fprintf(file,
" <energy>%1.12E kJ/mol</energy>\n",
2303 (0.5*atomEnergy[icalc][i]*conversion));
2304 fprintf(file,
" </atom>\n");
2310 fprintf(file,
" <qfforce_x>%1.12E</qfforce_x>\n",
2311 atomForce[icalc][0].qfForce[0]*conversion);
2312 fprintf(file,
" <qfforce_y>%1.12E</qfforce_y>\n",
2313 atomForce[icalc][0].qfForce[1]*conversion);
2314 fprintf(file,
" <qfforce_z>%1.12E</qfforce_z>\n",
2315 atomForce[icalc][0].qfForce[2]*conversion);
2316 fprintf(file,
" <ibforce_x>%1.12E</ibforce_x>\n",
2317 atomForce[icalc][0].ibForce[0]*conversion);
2318 fprintf(file,
" <ibforce_y>%1.12E</ibforce_y>\n",
2319 atomForce[icalc][0].ibForce[1]*conversion);
2320 fprintf(file,
" <ibforce_z>%1.12E</ibforce_z>\n",
2321 atomForce[icalc][0].ibForce[2]*conversion);
2322 fprintf(file,
" <dbforce_x>%1.12E</dbforce_x>\n",
2323 atomForce[icalc][0].dbForce[0]*conversion);
2324 fprintf(file,
" <dbforce_y>%1.12E</dbforce_y>\n",
2325 atomForce[icalc][0].dbForce[1]*conversion);
2326 fprintf(file,
" <dbforce_z>%1.12E</dbforce_z>\n",
2327 atomForce[icalc][0].dbForce[2]*conversion);
2330 fprintf(file,
" </calc>\n");
2333 fprintf(file,
" </elec>\n");
2338for (i=0; i<nosh->
nprint; i++) {
2342 fprintf(file,
" <printEnergy>\n");
2343 fprintf(file,
" <equation>%d", nosh->
printcalc[i][0]+1);
2346 if (nosh->
printop[i][j-1] == 0) fprintf(file,
" +");
2347 else if (nosh->
printop[i][j-1] == 1) fprintf(file,
" -");
2348 fprintf(file,
" %d", nosh->
printcalc[i][j] +1);
2351 fprintf(file,
"</equation>\n");
2360 if (nosh->
printop[i][j-1] == 0) scalar = 1.0;
2361 else if (nosh->
printop[i][j-1] == 1) scalar = -1.0;
2366 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2367 fprintf(file,
" <localEnergy>%1.12E kJ/mol</localEnergy>\n", \
2369 fprintf(file,
" <globalEnergy>%1.12E kJ/mol</globalEnergy>\n", \
2372 fprintf(file,
" </printEnergy>\n");
2377fprintf(file,
"</APBS>\n");
2389 char writestem[VMAX_ARGLEN];
2390 char outpath[VMAX_ARGLEN];
2410 if (nosh->
bogus)
return 1;
2412 for (i=0; i<pbeparm->
numwrite; i++) {
2425 Vnm_tprint(1,
" Writing charge distribution to ");
2429 xmin = xcent - 0.5*(nx-1)*hx;
2430 ymin = ycent - 0.5*(ny-1)*hy;
2431 zmin = zcent - 0.5*(nz-1)*hzed;
2434 sprintf(title,
"CHARGE DISTRIBUTION (e)");
2439 Vnm_tprint(1,
" Writing potential to ");
2443 xmin = xcent - 0.5*(nx-1)*hx;
2444 ymin = ycent - 0.5*(ny-1)*hy;
2445 zmin = zcent - 0.5*(nz-1)*hzed;
2448 sprintf(title,
"POTENTIAL (kT/e)");
2453 Vnm_tprint(1,
" Writing molecular accessibility to ");
2457 xmin = xcent - 0.5*(nx-1)*hx;
2458 ymin = ycent - 0.5*(ny-1)*hy;
2459 zmin = zcent - 0.5*(nz-1)*hzed;
2463 "SOLVENT ACCESSIBILITY -- MOLECULAR (%4.3f PROBE)",
2469 Vnm_tprint(1,
" Writing spline-based accessibility to ");
2473 xmin = xcent - 0.5*(nx-1)*hx;
2474 ymin = ycent - 0.5*(ny-1)*hy;
2475 zmin = zcent - 0.5*(nz-1)*hzed;
2479 "SOLVENT ACCESSIBILITY -- SPLINE (%4.3f WINDOW)",
2485 Vnm_tprint(1,
" Writing van der Waals accessibility to ");
2489 xmin = xcent - 0.5*(nx-1)*hx;
2490 ymin = ycent - 0.5*(ny-1)*hy;
2491 zmin = zcent - 0.5*(nz-1)*hzed;
2494 sprintf(title,
"SOLVENT ACCESSIBILITY -- VAN DER WAALS");
2499 Vnm_tprint(1,
" Writing ion accessibility to ");
2503 xmin = xcent - 0.5*(nx-1)*hx;
2504 ymin = ycent - 0.5*(ny-1)*hy;
2505 zmin = zcent - 0.5*(nz-1)*hzed;
2509 "ION ACCESSIBILITY -- SPLINE (%4.3f RADIUS)",
2515 Vnm_tprint(1,
" Writing potential Laplacian to ");
2519 xmin = xcent - 0.5*(nx-1)*hx;
2520 ymin = ycent - 0.5*(ny-1)*hy;
2521 zmin = zcent - 0.5*(nz-1)*hzed;
2525 "POTENTIAL LAPLACIAN (kT/e/A^2)");
2530 Vnm_tprint(1,
" Writing energy density to ");
2534 xmin = xcent - 0.5*(nx-1)*hx;
2535 ymin = ycent - 0.5*(ny-1)*hy;
2536 zmin = zcent - 0.5*(nz-1)*hzed;
2539 sprintf(title,
"ENERGY DENSITY (kT/e/A)^2");
2544 Vnm_tprint(1,
" Writing number density to ");
2548 xmin = xcent - 0.5*(nx-1)*hx;
2549 ymin = ycent - 0.5*(ny-1)*hy;
2550 zmin = zcent - 0.5*(nz-1)*hzed;
2554 "ION NUMBER DENSITY (M)");
2559 Vnm_tprint(1,
" Writing charge density to ");
2563 xmin = xcent - 0.5*(nx-1)*hx;
2564 ymin = ycent - 0.5*(ny-1)*hy;
2565 zmin = zcent - 0.5*(nz-1)*hzed;
2569 "ION CHARGE DENSITY (e_c * M)");
2574 Vnm_tprint(1,
" Writing x-shifted dielectric map to ");
2578 xmin = xcent - 0.5*(nx-1)*hx;
2579 ymin = ycent - 0.5*(ny-1)*hy;
2580 zmin = zcent - 0.5*(nz-1)*hzed;
2584 "X-SHIFTED DIELECTRIC MAP");
2589 Vnm_tprint(1,
" Writing y-shifted dielectric map to ");
2593 xmin = xcent - 0.5*(nx-1)*hx;
2594 ymin = ycent - 0.5*(ny-1)*hy;
2595 zmin = zcent - 0.5*(nz-1)*hzed;
2599 "Y-SHIFTED DIELECTRIC MAP");
2604 Vnm_tprint(1,
" Writing z-shifted dielectric map to ");
2608 xmin = xcent - 0.5*(nx-1)*hx;
2609 ymin = ycent - 0.5*(ny-1)*hy;
2610 zmin = zcent - 0.5*(nz-1)*hzed;
2614 "Z-SHIFTED DIELECTRIC MAP");
2619 Vnm_tprint(1,
" Writing kappa map to ");
2623 xmin = xcent - 0.5*(nx-1)*hx;
2624 ymin = ycent - 0.5*(ny-1)*hy;
2625 zmin = zcent - 0.5*(nz-1)*hzed;
2634 Vnm_tprint(1,
" Writing atom potentials to ");
2638 xmin = xcent - 0.5*(nx-1)*hx;
2639 ymin = ycent - 0.5*(ny-1)*hy;
2640 zmin = zcent - 0.5*(nz-1)*hzed;
2648 Vnm_tprint(2,
"Invalid data type for writing!\n");
2654 sprintf(writestem,
"%s-PE%d", pbeparm->
writestem[i], rank);
2659 sprintf(writestem,
"%s", pbeparm->
writestem[i]);
2666 sprintf(outpath,
"%s.%s", writestem,
"dx");
2667 Vnm_tprint(1,
"%s\n", outpath);
2668 grid =
Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2676 sprintf(outpath,
"%s.%s", writestem,
"dxbin");
2677 Vnm_tprint(1,
"%s\n", outpath);
2678 grid =
Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2687 sprintf(outpath,
"%s.%s", writestem,
"ucd");
2688 Vnm_tprint(1,
"%s\n", outpath);
2689 Vnm_tprint(2,
"Sorry, AVS format isn't supported for \
2690uniform meshes yet!\n");
2694 sprintf(outpath,
"%s.%s", writestem,
"mcsf");
2695 Vnm_tprint(1,
"%s\n", outpath);
2696 Vnm_tprint(2,
"Sorry, MCSF format isn't supported for \
2697 uniform meshes yet!\n");
2701 sprintf(outpath,
"%s.%s", writestem,
"grd");
2702 Vnm_tprint(1,
"%s\n", outpath);
2703 grid =
Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2711 sprintf(outpath,
"%s.%s", writestem,
"dx.gz");
2712 Vnm_tprint(1,
"%s\n", outpath);
2713 grid =
Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2720 sprintf(outpath,
"%s.%s", writestem,
"txt");
2721 Vnm_tprint(1,
"%s\n", outpath);
2722 Vnm_print(0,
"routines: Opening virtual socket...\n");
2723 sock = Vio_ctor(
"FILE",
"ASC",VNULL,outpath,
"w");
2724 if (sock == VNULL) {
2725 Vnm_print(2,
"routines: Problem opening virtual socket %s\n",
2729 if (Vio_connect(sock, 0) < 0) {
2730 Vnm_print(2,
"routines: Problem connecting virtual socket %s\n",
2734 Vio_printf(sock,
"# Data from %s\n", PACKAGE_STRING);
2735 Vio_printf(sock,
"# \n");
2736 Vio_printf(sock,
"# %s\n", title);
2737 Vio_printf(sock,
"# \n");
2739 for(i=0;i<natoms;i++)
2740 Vio_printf(sock,
"%12.6e\n", pmg->
rwork[i]);
2743 Vnm_tprint(2,
"Bogus data format (%d)!\n",
2769 Vnm_tprint( 2,
" No energy available in Calculation %d\n", calcid+1);
2772 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++){
2775 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2776 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2797 Vnm_tprint( 2,
"Warning: The 'energy' print keyword is deprecated.\n" \
2798 " Use eilecEnergy for electrostatics energy calcs.\n\n");
2801 Vnm_tprint( 1,
"print energy %d ", nosh->
printcalc[iprint][0]+1);
2803 Vnm_tprint( 1,
"print energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2806 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2807 if (nosh->
printop[iprint][iarg-1] == 0)
2808 Vnm_tprint(1,
"+ ");
2809 else if (nosh->
printop[iprint][iarg-1] == 1)
2810 Vnm_tprint(1,
"- ");
2812 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2817 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2819 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2823 Vnm_tprint(1,
"end\n");
2829 Vnm_tprint( 2,
" Didn't calculate energy in Calculation \
2833 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2836 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2837 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2843 Vnm_tprint( 1,
" Local net energy (PE %d) = %1.12E kJ/mol\n",
2844 Vcom_rank(com), ltenergy);
2845 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
2846 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2847 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2866 Vnm_tprint( 1,
"\nprint energy %d ", nosh->
printcalc[iprint][0]+1);
2868 Vnm_tprint( 1,
"\nprint energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2871 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2872 if (nosh->
printop[iprint][iarg-1] == 0)
2873 Vnm_tprint(1,
"+ ");
2874 else if (nosh->
printop[iprint][iarg-1] == 1)
2875 Vnm_tprint(1,
"- ");
2877 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2882 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2884 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2888 Vnm_tprint(1,
"end\n");
2894 Vnm_tprint( 2,
" Didn't calculate energy in Calculation \
2898 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2901 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2902 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2908 Vnm_tprint( 1,
" Local net energy (PE %d) = %1.12E kJ/mol\n",
2909 Vcom_rank(com), ltenergy);
2910 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
2911 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2912 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2929 Vnm_tprint( 1,
"\nprint APOL energy %d ", nosh->
printcalc[iprint][0]+1);
2931 Vnm_tprint( 1,
"\nprint APOL energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2934 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2935 if (nosh->
printop[iprint][iarg-1] == 0)
2936 Vnm_tprint(1,
"+ ");
2937 else if (nosh->
printop[iprint][iarg-1] == 1)
2938 Vnm_tprint(1,
"- ");
2940 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2945 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2947 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2951 Vnm_tprint(1,
"end\n");
2959 Vnm_tprint( 2,
" Didn't calculate energy in Calculation #%d\n", calcid+1);
2962 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2967 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2968 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2970 gtenergy += (scalar * ((apolparm->
gamma*apolparm->
sasa) +
2976 Vnm_tprint( 1,
" Global net APOL energy = %1.12E kJ/mol\n", gtenergy);
3000 Vnm_tprint( 2,
"Warning: The 'force' print keyword is deprecated.\n" \
3001 " Use elecForce for electrostatics force calcs.\n\n");
3004 Vnm_tprint( 1,
"print force %d ", nosh->
printcalc[iprint][0]+1);
3006 Vnm_tprint( 1,
"print force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3009 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3010 if (nosh->
printop[iprint][iarg-1] == 0)
3011 Vnm_tprint(1,
"+ ");
3012 else if (nosh->
printop[iprint][iarg-1] == 1)
3013 Vnm_tprint(1,
"- ");
3015 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3020 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3022 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3026 Vnm_tprint(1,
"end\n");
3031 refnforce = nforce[calcid];
3033 if (refcalcforce ==
PCF_NO) {
3034 Vnm_tprint( 2,
" Didn't calculate force in calculation \
3038 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3041 Vnm_tprint(2,
" Inconsistent calcforce declarations in \
3046 if (nforce[calcid] != refnforce) {
3047 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \
3060 aforce = atomForce[calcid];
3066 for (ivc=0; ivc<3; ivc++) {
3075 for (ifr=0; ifr<refnforce; ifr++) {
3076 for (ivc=0; ivc<3; ivc++) {
3088 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3091 aforce = atomForce[calcid];
3093 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3094 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3099 for (ivc=0; ivc<3; ivc++) {
3108 for (ifr=0; ifr<refnforce; ifr++) {
3109 for (ivc=0; ivc<3; ivc++) {
3121 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
3122 for (ifr=0; ifr<refnforce; ifr++) {
3123 Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
3124 Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
3125 Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
3130 Vnm_tprint( 1,
" Local net fixed charge force = \
3131(%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
3132 lforce[0].qfForce[1], lforce[0].qfForce[2]);
3133 Vnm_tprint( 1,
" Local net ionic boundary force = \
3134(%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
3135 lforce[0].ibForce[1], lforce[0].ibForce[2]);
3136 Vnm_tprint( 1,
" Local net dielectric boundary force = \
3137(%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
3138 lforce[0].dbForce[1], lforce[0].dbForce[2]);
3140 for (ifr=0; ifr<refnforce; ifr++) {
3141 Vnm_tprint( 1,
" Local fixed charge force \
3142(atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
3143 lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
3144 Vnm_tprint( 1,
" Local ionic boundary force \
3145(atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
3146 lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
3147 Vnm_tprint( 1,
" Local dielectric boundary force \
3148(atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
3149 lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
3155 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A).\n");
3156 Vnm_tprint( 1,
" Legend:\n");
3157 Vnm_tprint( 1,
" tot -- Total force\n");
3158 Vnm_tprint( 1,
" qf -- Fixed charge force\n");
3159 Vnm_tprint( 1,
" db -- Dielectric boundary force\n");
3160 Vnm_tprint( 1,
" ib -- Ionic boundary force\n");
3162 for (ivc=0; ivc<3; ivc++) {
3168 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3169 totforce[1], totforce[2]);
3170 Vnm_tprint( 1,
" qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3171 gforce[0].qfForce[1], gforce[0].qfForce[2]);
3172 Vnm_tprint( 1,
" ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3173 gforce[0].ibForce[1], gforce[0].ibForce[2]);
3174 Vnm_tprint( 1,
" db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3175 gforce[0].dbForce[1], gforce[0].dbForce[2]);
3179 Vnm_tprint( 1,
" Printing per-atom forces (kJ/mol/A).\n");
3180 Vnm_tprint( 1,
" Legend:\n");
3181 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3182 Vnm_tprint( 1,
" qf n -- Fixed charge force for atom n\n");
3183 Vnm_tprint( 1,
" db n -- Dielectric boundary force for atom n\n");
3184 Vnm_tprint( 1,
" ib n -- Ionic boundary force for atom n\n");
3185 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3191 for (ifr=0; ifr<refnforce; ifr++) {
3192 Vnm_tprint( 1,
" qf %d %1.12E %1.12E %1.12E\n", ifr,
3193 gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3194 gforce[ifr].qfForce[2]);
3195 Vnm_tprint( 1,
" ib %d %1.12E %1.12E %1.12E\n", ifr,
3196 gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3197 gforce[ifr].ibForce[2]);
3198 Vnm_tprint( 1,
" db %d %1.12E %1.12E %1.12E\n", ifr,
3199 gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3200 gforce[ifr].dbForce[2]);
3201 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3202 (gforce[ifr].dbForce[0] \
3203 + gforce[ifr].ibForce[0] +
3204 gforce[ifr].qfForce[0]),
3205 (gforce[ifr].dbForce[1] \
3206 + gforce[ifr].ibForce[1] +
3207 gforce[ifr].qfForce[1]),
3208 (gforce[ifr].dbForce[2] \
3209 + gforce[ifr].ibForce[2] +
3210 gforce[ifr].qfForce[2]));
3211 for (ivc=0; ivc<3; ivc++) {
3212 totforce[ivc] += (gforce[ifr].
dbForce[ivc] \
3217 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3218 totforce[1], totforce[2]);
3221 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3222 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3247 Vnm_tprint( 1,
"print force %d ", nosh->
printcalc[iprint][0]+1);
3249 Vnm_tprint( 1,
"print force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3252 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3253 if (nosh->
printop[iprint][iarg-1] == 0)
3254 Vnm_tprint(1,
"+ ");
3255 else if (nosh->
printop[iprint][iarg-1] == 1)
3256 Vnm_tprint(1,
"- ");
3258 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3263 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3265 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3269 Vnm_tprint(1,
"end\n");
3274 refnforce = nforce[calcid];
3276 if (refcalcforce ==
PCF_NO) {
3277 Vnm_tprint( 2,
" Didn't calculate force in calculation \
3281 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3284 Vnm_tprint(2,
" Inconsistent calcforce declarations in \
3289 if (nforce[calcid] != refnforce) {
3290 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \
3303 aforce = atomForce[calcid];
3309 for (ivc=0; ivc<3; ivc++) {
3318 for (ifr=0; ifr<refnforce; ifr++) {
3319 for (ivc=0; ivc<3; ivc++) {
3331 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3334 aforce = atomForce[calcid];
3336 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3337 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3342 for (ivc=0; ivc<3; ivc++) {
3351 for (ifr=0; ifr<refnforce; ifr++) {
3352 for (ivc=0; ivc<3; ivc++) {
3364 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
3365 for (ifr=0; ifr<refnforce; ifr++) {
3366 Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
3367 Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
3368 Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
3373 Vnm_tprint( 1,
" Local net fixed charge force = \
3374(%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
3375 lforce[0].qfForce[1], lforce[0].qfForce[2]);
3376 Vnm_tprint( 1,
" Local net ionic boundary force = \
3377(%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
3378 lforce[0].ibForce[1], lforce[0].ibForce[2]);
3379 Vnm_tprint( 1,
" Local net dielectric boundary force = \
3380(%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
3381 lforce[0].dbForce[1], lforce[0].dbForce[2]);
3383 for (ifr=0; ifr<refnforce; ifr++) {
3384 Vnm_tprint( 1,
" Local fixed charge force \
3385(atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
3386 lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
3387 Vnm_tprint( 1,
" Local ionic boundary force \
3388(atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
3389 lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
3390 Vnm_tprint( 1,
" Local dielectric boundary force \
3391(atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
3392 lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
3398 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A).\n");
3399 Vnm_tprint( 1,
" Legend:\n");
3400 Vnm_tprint( 1,
" tot -- Total force\n");
3401 Vnm_tprint( 1,
" qf -- Fixed charge force\n");
3402 Vnm_tprint( 1,
" db -- Dielectric boundary force\n");
3403 Vnm_tprint( 1,
" ib -- Ionic boundary force\n");
3405 for (ivc=0; ivc<3; ivc++) {
3411 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3412 totforce[1], totforce[2]);
3413 Vnm_tprint( 1,
" qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3414 gforce[0].qfForce[1], gforce[0].qfForce[2]);
3415 Vnm_tprint( 1,
" ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3416 gforce[0].ibForce[1], gforce[0].ibForce[2]);
3417 Vnm_tprint( 1,
" db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3418 gforce[0].dbForce[1], gforce[0].dbForce[2]);
3422 Vnm_tprint( 1,
" Printing per-atom forces (kJ/mol/A).\n");
3423 Vnm_tprint( 1,
" Legend:\n");
3424 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3425 Vnm_tprint( 1,
" qf n -- Fixed charge force for atom n\n");
3426 Vnm_tprint( 1,
" db n -- Dielectric boundary force for atom n\n");
3427 Vnm_tprint( 1,
" ib n -- Ionic boundary force for atom n\n");
3428 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3434 for (ifr=0; ifr<refnforce; ifr++) {
3435 Vnm_tprint( 1,
" qf %d %1.12E %1.12E %1.12E\n", ifr,
3436 gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3437 gforce[ifr].qfForce[2]);
3438 Vnm_tprint( 1,
" ib %d %1.12E %1.12E %1.12E\n", ifr,
3439 gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3440 gforce[ifr].ibForce[2]);
3441 Vnm_tprint( 1,
" db %d %1.12E %1.12E %1.12E\n", ifr,
3442 gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3443 gforce[ifr].dbForce[2]);
3444 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3445 (gforce[ifr].dbForce[0] \
3446 + gforce[ifr].ibForce[0] +
3447 gforce[ifr].qfForce[0]),
3448 (gforce[ifr].dbForce[1] \
3449 + gforce[ifr].ibForce[1] +
3450 gforce[ifr].qfForce[1]),
3451 (gforce[ifr].dbForce[2] \
3452 + gforce[ifr].ibForce[2] +
3453 gforce[ifr].qfForce[2]));
3454 for (ivc=0; ivc<3; ivc++) {
3455 totforce[ivc] += (gforce[ifr].
dbForce[ivc] \
3460 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3461 totforce[1], totforce[2]);
3464 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3465 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3492 Vnm_tprint( 1,
"\nprint APOL force %d ", nosh->
printcalc[iprint][0]+1);
3494 Vnm_tprint( 1,
"\nprint APOL force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3497 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3498 if (nosh->
printop[iprint][iarg-1] == 0)
3499 Vnm_tprint(1,
"+ ");
3500 else if (nosh->
printop[iprint][iarg-1] == 1)
3501 Vnm_tprint(1,
"- ");
3503 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3508 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3510 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3514 Vnm_tprint(1,
"end\n");
3519 refnforce = nforce[calcid];
3521 if (refcalcforce ==
ACF_NO) {
3522 Vnm_tprint( 2,
" Didn't calculate force in calculation \
3526 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3529 Vnm_tprint(2,
" Inconsistent calcforce declarations in \
3534 if (nforce[calcid] != refnforce) {
3535 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \
3548 aforce = atomForce[calcid];
3554 for (ivc=0; ivc<3; ivc++) {
3560 for (ifr=0; ifr<refnforce; ifr++) {
3561 for (ivc=0; ivc<3; ivc++) {
3570 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3573 aforce = atomForce[calcid];
3575 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3576 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3581 for (ivc=0; ivc<3; ivc++) {
3587 for (ifr=0; ifr<refnforce; ifr++) {
3588 for (ivc=0; ivc<3; ivc++) {
3597 Vnm_tprint( 0,
"printForce: Performing global reduction (sum)\n");
3598 for (ifr=0; ifr<refnforce; ifr++) {
3599 Vcom_reduce(com, lforce[ifr].sasaForce, gforce[ifr].sasaForce, 3, 2, 0);
3600 Vcom_reduce(com, lforce[ifr].savForce, gforce[ifr].savForce, 3, 2, 0);
3601 Vcom_reduce(com, lforce[ifr].wcaForce, gforce[ifr].wcaForce, 3, 2, 0);
3605 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A)\n");
3606 Vnm_tprint( 1,
" Legend:\n");
3607 Vnm_tprint( 1,
" tot -- Total force\n");
3608 Vnm_tprint( 1,
" sasa -- SASA force\n");
3609 Vnm_tprint( 1,
" sav -- SAV force\n");
3610 Vnm_tprint( 1,
" wca -- WCA force\n\n");
3612 for (ivc=0; ivc<3; ivc++) {
3618 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3619 totforce[1], totforce[2]);
3620 Vnm_tprint( 1,
" sasa %1.12E %1.12E %1.12E\n", gforce[0].sasaForce[0],
3621 gforce[0].sasaForce[1], gforce[0].sasaForce[2]);
3622 Vnm_tprint( 1,
" sav %1.12E %1.12E %1.12E\n", gforce[0].savForce[0],
3623 gforce[0].savForce[1], gforce[0].savForce[2]);
3624 Vnm_tprint( 1,
" wca %1.12E %1.12E %1.12E\n", gforce[0].wcaForce[0],
3625 gforce[0].wcaForce[1], gforce[0].wcaForce[2]);
3629 Vnm_tprint( 1,
" Printing per atom forces (kJ/mol/A)\n");
3630 Vnm_tprint( 1,
" Legend:\n");
3631 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3632 Vnm_tprint( 1,
" sasa n -- SASA force for atom n\n");
3633 Vnm_tprint( 1,
" sav n -- SAV force for atom n\n");
3634 Vnm_tprint( 1,
" wca n -- WCA force for atom n\n");
3635 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3644 for (ifr=0; ifr<refnforce; ifr++) {
3645 Vnm_tprint( 1,
" sasa %d %1.12E %1.12E %1.12E\n", ifr,
3646 gforce[ifr].sasaForce[0], gforce[ifr].sasaForce[1],
3647 gforce[ifr].sasaForce[2]);
3648 Vnm_tprint( 1,
" sav %d %1.12E %1.12E %1.12E\n", ifr,
3649 gforce[ifr].savForce[0], gforce[ifr].savForce[1],
3650 gforce[ifr].savForce[2]);
3651 Vnm_tprint( 1,
" wca %d %1.12E %1.12E %1.12E\n", ifr,
3652 gforce[ifr].wcaForce[0], gforce[ifr].wcaForce[1],
3653 gforce[ifr].wcaForce[2]);
3654 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3655 (gforce[ifr].wcaForce[0] \
3656 + gforce[ifr].savForce[0] +
3657 gforce[ifr].sasaForce[0]),
3658 (gforce[ifr].wcaForce[1] \
3659 + gforce[ifr].savForce[1] +
3660 gforce[ifr].sasaForce[1]),
3661 (gforce[ifr].wcaForce[2] \
3662 + gforce[ifr].savForce[2] +
3663 gforce[ifr].sasaForce[2]));
3664 for (ivc=0; ivc<3; ivc++) {
3665 totforce[ivc] += (gforce[ifr].
wcaForce[ivc] \
3670 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3671 totforce[1], totforce[2]);
3674 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3675 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3692 Vnm_tprint(1,
"Destroying finite element structures.\n");
3695 for(i=0;i<nosh->
ncalc;i++){
3699 for (i=0; i<nosh->
nmesh; i++) {
3737 Vatom *atom = VNULL;
3739 Vnm_tstart(27,
"Setup timer");
3742 if (pbeparm->
useDielMap) Vnm_tprint(2,
"FEM ignoring dielectric map!\n");
3743 if (pbeparm->
useKappaMap) Vnm_tprint(2,
"FEM ignoring kappa map!\n");
3744 if (pbeparm->
useChargeMap) Vnm_tprint(2,
"FEM ignoring charge map!\n");
3747 Vnm_tprint(0,
"Re-centering mesh...\n");
3748 theMol = pbeparm->
molid-1;
3749 myalist = alist[theMol];
3750 for (j=0; j<3; j++) {
3751 if (theMol < nosh->nmol) {
3752 center[j] = (myalist)->center[j];
3754 Vnm_tprint(2,
"ERROR! Bogus molecule number (%d)!\n",
3782 Vnm_tprint(0,
"Setting up PBE object...\n");
3784 sparm = pbeparm->
swin;
3787 sparm = pbeparm->
srad;
3789 if (pbeparm->
nion > 0) {
3790 iparm = pbeparm->
ionr[0];
3796 pbeparm->
sdie, sparm, focusFlag, pbeparm->
sdens,
3801 Vnm_tprint(1,
" Debye length: %g A\n",
Vpbe_getDeblen(pbe[icalc]));
3804 Vnm_tprint(0,
"Setting up FEtk object...\n");
3811 Vnm_tprint(0,
"Setting up mesh...\n");
3814 imesh = feparm->
meshID-1;
3816 for (i=0; i<3; i++) {
3820 Vnm_print(0,
"Using mesh %d (%s) in calculation.\n", imesh+1,
3822 switch (nosh->
meshfmt[imesh]) {
3824 Vnm_tprint(2,
"DX finite element mesh input not supported yet!\n");
3827 Vnm_tprint(2,
"DXBIN finite element mesh input not supported yet!\n");
3830 Vnm_tprint( 2,
"UHBD finite element mesh input not supported!\n");
3833 Vnm_tprint( 2,
"AVS finite element mesh input not supported!\n");
3836 Vnm_tprint(1,
"Reading MCSF-format input finite element mesh from %s.\n",
3838 sock = Vio_ctor(
"FILE",
"ASC", VNULL, nosh->
meshpath[imesh],
"r");
3839 if (sock == VNULL) {
3840 Vnm_print(2,
"Problem opening virtual socket %s!\n",
3844 if (Vio_accept(sock, 0) < 0) {
3845 Vnm_print(2,
"Problem accepting virtual socket %s!\n",
3852 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
3858 for (i=0; i<3; i++) {
3859 center[i] = alist[theMol]->center[i];
3860 length[i] = feparm->
glen[i];
3865 vrc =
Vfetk_loadMesh(fetk[icalc], center, length, meshType, sock);
3867 Vnm_print(2,
"Error constructing finite element mesh!\n");
3871 Gem_shapeChk(fetk[icalc]->gm);
3875 for (j=0; j<2; j++) {
3879 AM_markRefine(fetk[icalc]->am, 0, -1, 0, 0.0);
3881 AM_refine(fetk[icalc]->am, 2, 0, feparm->
pkey);
3883 Gem_shapeChk(fetk[icalc]->gm);
3888 Vnm_tstop(27,
"Setup timer");
3891 bytesTotal = Vmem_bytesTotal();
3892 highWater = Vmem_highWaterTotal();
3895 Vnm_tprint( 1,
" Current memory usage: %4.3f MB total, \
3896%4.3f MB high water\n", (
double)(bytesTotal)/(1024.*1024.),
3897 (
double)(highWater)/(1024.*1024.));
3909 Vnm_tprint(1,
" Domain size: %g A x %g A x %g A\n",
3912 switch(feparm->
ekey) {
3914 Vnm_tprint(1,
" Per-simplex error tolerance: %g\n", feparm->
etol);
3917 Vnm_tprint(1,
" Global error tolerance: %g\n", feparm->
etol);
3920 Vnm_tprint(1,
" Fraction of simps to refine: %g\n", feparm->
etol);
3923 Vnm_tprint(2,
"Invalid ekey (%d)!\n", feparm->
ekey);
3929 Vnm_tprint(1,
" Uniform pre-solve refinement.\n");
3932 Vnm_tprint(1,
" Geometry-based pre-solve refinement.\n");
3935 Vnm_tprint(2,
"Invalid akeyPRE (%d)!\n", feparm->
akeyPRE);
3941 Vnm_tprint(1,
" Residual-based a posteriori refinement.\n");
3944 Vnm_tprint(1,
" Dual-based a posteriori refinement.\n");
3947 Vnm_tprint(1,
" Local-based a posteriori refinement.\n");
3950 Vnm_tprint(2,
"Invalid akeySOLVE (%d)!\n", feparm->
akeySOLVE);
3953 Vnm_tprint(1,
" Refinement of initial mesh to ~%d vertices\n",
3955 Vnm_tprint(1,
" Geometry-based refinment lower bound: %g A\n",
3957 Vnm_tprint(1,
" Maximum number of solve-estimate-refine cycles: %d\n",
3959 Vnm_tprint(1,
" Maximum number of vertices in mesh: %d\n",
3963 if (nosh->
bogus)
return;
3965 Vnm_tprint(1,
" HB linear solver: AM_hPcg\n");
3967 Vnm_tprint(1,
" Non-HB linear solver: ");
3968 switch (fetk[icalc]->lkey) {
3970 Vnm_print(1,
"SLU direct\n");
3973 Vnm_print(1,
"multigrid\n");
3976 Vnm_print(1,
"conjugate gradient\n");
3979 Vnm_print(1,
"BiCGStab\n");
3982 Vnm_print(1,
"???\n");
3987 Vnm_tprint(1,
" Linear solver tol.: %g\n", fetk[icalc]->ltol);
3988 Vnm_tprint(1,
" Linear solver max. iters.: %d\n", fetk[icalc]->lmax);
3989 Vnm_tprint(1,
" Linear solver preconditioner: ");
3990 switch (fetk[icalc]->lprec) {
3992 Vnm_print(1,
"identity\n");
3995 Vnm_print(1,
"diagonal\n");
3998 Vnm_print(1,
"multigrid\n");
4001 Vnm_print(1,
"???\n");
4004 Vnm_tprint(1,
" Nonlinear solver: ");
4005 switch (fetk[icalc]->nkey) {
4007 Vnm_print(1,
"newton\n");
4010 Vnm_print(1,
"incremental\n");
4013 Vnm_print(1,
"pseudo-arclength\n");
4016 Vnm_print(1,
"??? ");
4019 Vnm_tprint(1,
" Nonlinear solver tol.: %g\n", fetk[icalc]->ntol);
4020 Vnm_tprint(1,
" Nonlinear solver max. iters.: %d\n", fetk[icalc]->nmax);
4021 Vnm_tprint(1,
" Initial guess: ");
4022 switch (fetk[icalc]->gues) {
4024 Vnm_tprint(1,
"zero\n");
4027 Vnm_tprint(1,
"boundary function\n");
4030 Vnm_tprint(1,
"interpolated previous solution\n");
4033 Vnm_tprint(1,
"???\n");
4058 Vnm_tprint(1,
" Commencing uniform refinement to %d verts.\n",
4062 Vnm_tprint(1,
" Commencing geometry-based refinement to %d \
4066 Vnm_tprint(2,
"What? You can't do a posteriori error estimation \
4067before you solve!\n");
4082 Vnm_tprint(1,
" Initial mesh has %d vertices\n",
4083 Gem_numVV(fetk[icalc]->gm));
4089 nverts = Gem_numVV(fetk[icalc]->gm);
4091 Vnm_tprint(1,
" Hit vertex number limit.\n");
4094 marked = AM_markRefine(fetk[icalc]->am, feparm->
akeyPRE, -1,
4097 Vnm_tprint(1,
" Marked 0 simps; hit error/size tolerance.\n");
4100 Vnm_tprint(1,
" Have %d verts, marked %d. Refining...\n", nverts,
4102 AM_refine(fetk[icalc]->am, 0, 0, feparm->
pkey);
4105 nverts = Gem_numVV(fetk[icalc]->gm);
4106 Vnm_tprint(1,
" Done refining; have %d verts.\n", nverts);
4128 (pbeparm->
pbetype == PBE_NRPBE) ||
4151 Vnm_print(2,
"SORRY! DON'T USE HB!!!\n");
4155 AM_hlSolve(fetk[icalc]->am, meth, lkeyHB, fetk[icalc]->lmax,
4156 fetk[icalc]->ltol, fetk[icalc]->gues, fetk[icalc]->pjac);
4201 if (nosh->
bogus == 0) {
4203 (pbeparm->
pbetype==PBE_NRPBE) ||
4214 Vnm_tprint(1,
" Total electrostatic energy = %1.12E kJ/mol\n",
4221 Vnm_tprint(2,
"Error! Verbose energy evaluation not available for FEM yet!\n");
4222 Vnm_tprint(2,
"E-mail nathan.baker@pnl.gov if you want this.\n");
4247 nverts = Gem_numVV(fetk[icalc]->gm);
4248 if (nverts > feparm->
maxvert) {
4249 Vnm_tprint(1,
" Current number of vertices (%d) exceeds max (%d)!\n",
4253 Vnm_tprint(1,
" Mesh currently has %d vertices\n", nverts);
4257 Vnm_tprint(1,
" Commencing uniform refinement.\n");
4260 Vnm_tprint(1,
" Commencing geometry-based refinement.\n");
4263 Vnm_tprint(1,
" Commencing residual-based refinement.\n");
4266 Vnm_tprint(1,
" Commencing dual problem-based refinement.\n");
4269 Vnm_tprint(1,
" Commencing local-based refinement.\n.");
4272 Vnm_tprint(2,
" Error -- unknown refinement type (%d)!\n",
4279 marked = AM_markRefine(fetk[icalc]->am, feparm->
akeySOLVE, -1,
4282 Vnm_tprint(1,
" Marked 0 simps; hit error/size tolerance.\n");
4285 Vnm_tprint(1,
" Have %d verts, marked %d. Refining...\n", nverts,
4287 AM_refine(fetk[icalc]->am, 0, 0, feparm->
pkey);
4288 nverts = Gem_numVV(fetk[icalc]->gm);
4289 Vnm_tprint(1,
" Done refining; have %d verts.\n", nverts);
4291 Gem_shapeChk(fetk[icalc]->gm);
4306 char writestem[VMAX_ARGLEN];
4307 char outpath[VMAX_ARGLEN];
4313 if (nosh->
bogus)
return 1;
4318 for (i=0; i<pbeparm->
numwrite; i++) {
4326 Vnm_tprint(2,
" Sorry; can't write charge distribution for FEM!\n");
4332 Vnm_tprint(1,
" Writing potential to ");
4338 Vnm_tprint(1,
" Writing molecular accessibility to ");
4344 Vnm_tprint(1,
" Writing spline-based accessibility to ");
4350 Vnm_tprint(1,
" Writing van der Waals accessibility to ");
4356 Vnm_tprint(1,
" Writing ion accessibility to ");
4362 Vnm_tprint(2,
" Sorry; can't write charge distribution for FEM!\n");
4368 Vnm_tprint(2,
" Sorry; can't write energy density for FEM!\n");
4374 Vnm_tprint(1,
" Writing number density to ");
4380 Vnm_tprint(1,
" Writing charge density to ");
4386 Vnm_tprint(2,
" Sorry; can't write x-shifted dielectric map for FEM!\n");
4392 Vnm_tprint(2,
" Sorry; can't write y-shifted dielectric map for FEM!\n");
4398 Vnm_tprint(2,
" Sorry; can't write z-shifted dielectric map for FEM!\n");
4404 Vnm_tprint(1,
" Sorry; can't write kappa map for FEM!\n");
4410 Vnm_tprint(1,
" Sorry; can't write atom potentials for FEM!\n");
4416 Vnm_tprint(2,
"Invalid data type for writing!\n");
4421 if (!writeit)
return 0;
4425 sprintf(writestem,
"%s-PE%d", pbeparm->
writestem[i], rank);
4430 sprintf(writestem,
"%s", pbeparm->
writestem[i]);
4437 sprintf(outpath,
"%s.%s", writestem,
"dx");
4438 Vnm_tprint(1,
"%s\n", outpath);
4444 sprintf(outpath,
"%s.%s", writestem,
"dxbin");
4445 Vnm_tprint(1,
"%s\n", outpath);
4450 sprintf(outpath,
"%s.%s", writestem,
"ucd");
4451 Vnm_tprint(1,
"%s\n", outpath);
4456 Vnm_tprint(2,
"UHBD format not supported for FEM!\n");
4460 Vnm_tprint(2,
"MCSF format not supported yet!\n");
4464 Vnm_tprint(2,
"Bogus data format (%d)!\n",
4492 Vatom *atom = VNULL;
4535 for (i=0; i < natoms; i++) {
4541 if ((x+atomRadius) > xmax) xmax = x + atomRadius;
4542 if ((x-atomRadius) < xmin) xmin = x - atomRadius;
4543 if ((y+atomRadius) > ymax) ymax = y + atomRadius;
4544 if ((y-atomRadius) < ymin) ymin = y - atomRadius;
4545 if ((z+atomRadius) > zmax) zmax = z + atomRadius;
4546 if ((z-atomRadius) < zmin) zmin = z - atomRadius;
4547 disp[0] = (x - center[0]);
4548 disp[1] = (y - center[1]);
4549 disp[2] = (z - center[2]);
4550 dist = (disp[0]*disp[0]) + (disp[1]*disp[1]) + (disp[2]*disp[2]);
4551 dist = VSQRT(dist) + atomRadius;
4554 soluteXlen = xmax - xmin;
4555 soluteYlen = ymax - ymin;
4556 soluteZlen = zmax - zmin;
4559 Vnm_print(0,
"APOL: Setting up hash table and accessibility object...\n");
4560 nhash[0] = soluteXlen/0.5;
4561 nhash[1] = soluteYlen/0.5;
4562 nhash[2] = soluteZlen/0.5;
4563 for (i=0; i<3; i++) inhash[i] = (
int)(nhash[i]);
4566 if (inhash[i] < 3) inhash[i] = 3;
4567 if (inhash[i] > MAX_HASH_DIM) inhash[i] = MAX_HASH_DIM;
4571 srad = apolparm->
srad;
4572 sradPad = srad + (2*apolparm->
dpos);
4578 if (param == VNULL && (apolparm->
bconc != 0.0)) {
4579 Vnm_tprint(2,
"initAPOL: Got NULL Vparam object!\n");
4580 Vnm_tprint(2,
"initAPOL: You are performing an apolar calculation with the van der Waals integral term,\n");
4581 Vnm_tprint(2,
"initAPOL: this term requires van der Waals parameters which are not available from the \n");
4582 Vnm_tprint(2,
"initAPOL: PQR file. Therefore, you need to supply a parameter file with the parm keyword,\n");
4583 Vnm_tprint(2,
"initAPOL: for example,\n");
4584 Vnm_tprint(2,
"initAPOL: read parm flat amber94.dat end\n");
4585 Vnm_tprint(2,
"initAPOL: where the relevant parameter files can be found in apbs/tools/conversion/param/vparam.\n");
4589 if (apolparm->
bconc != 0.0){
4592 if (atomData == VNULL) {
4593 Vnm_tprint(2,
"initAPOL: Couldn't find parameters for WAT OW or WAT O!\n");
4594 Vnm_tprint(2,
"initAPOL: These parameters must be present in your file\n");
4595 Vnm_tprint(2,
"initAPOL: for apolar calculations.\n");
4605 rc =
forceAPOL(acc, mem, apolparm, nforce, atomForce, alist, clist);
4607 Vnm_print(2,
"Error in apolar force calculation!\n");
4619 if (VABS(apolparm->
gamma) > VSMALL) {
4623 for (i = 0; i < len; i++) {
4629 apolparm->
sasa = 0.0;
4631 for (i = 0; i < len; i++) {
4637 if (VABS(apolparm->
press) > VSMALL){
4640 apolparm->
sav = 0.0;
4644 if (VABS(apolparm->
bconc) > VSMALL) {
4646 for (i = 0; i < len; i++) {
4649 Vnm_print(2,
"Error in apolar energy calculation!\n");
4652 atomwcaEnergy[i] = energy;
4657 Vnm_print(2,
"Error in apolar energy calculation!\n");
4678 double atomwcaEnergy[],
4682 double energy = 0.0;
4686 Vnm_print(1,
"\nSolvent Accessible Surface Area (SASA) for each atom:\n");
4687 for (i = 0; i < numatoms; i++) {
4688 Vnm_print(1,
" SASA for atom %i: %1.12E\n", i, atomsasa[i]);
4691 Vnm_print(1,
"\nTotal solvent accessible surface area: %g A^2\n",sasa);
4698 Vnm_print(1,
"energyAPOL: Cannot calculate component energy, skipping.\n");
4701 energy = (apolparm->
gamma*sasa) + (apolparm->
press*sav)
4704 Vnm_print(1,
"\nSurface tension*area energies (gamma * SASA) for each atom:\n");
4705 for (i = 0; i < numatoms; i++) {
4706 Vnm_print(1,
" Surface tension*area energy for atom %i: %1.12E\n", i, apolparm->
gamma*atomsasa[i]);
4709 Vnm_print(1,
"\nTotal surface tension energy: %g kJ/mol\n", apolparm->
gamma*sasa);
4710 Vnm_print(1,
"\nTotal solvent accessible volume: %g A^3\n", sav);
4711 Vnm_print(1,
"\nTotal pressure*volume energy: %g kJ/mol\n", apolparm->
press*sav);
4712 Vnm_print(1,
"\nWCA dispersion Energies for each atom:\n");
4713 for (i = 0; i < numatoms; i++) {
4714 Vnm_print(1,
" WCA energy for atom %i: %1.12E\n", i, atomwcaEnergy[i]);
4717 Vnm_print(1,
"\nTotal WCA energy: %g kJ/mol\n", (apolparm->
wcaEnergy));
4718 Vnm_print(1,
"\nTotal non-polar energy = %1.12E kJ/mol\n", energy);
4722 Vnm_print(2,
"energyAPOL: Error in energyAPOL. Unknown option.\n");
4737 time_t ts, ts_main, ts_sub;
4755 Vatom *atom = VNULL;
4758 srad = apolparm->
srad;
4759 press = apolparm->
press;
4760 gamma = apolparm->
gamma;
4761 offset = apolparm->
dpos;
4762 bconc = apolparm->
bconc;
4767 Vnm_print(0,
"forceAPOL: Trying atom surf...\n");
4769 if (acc->
surf == VNULL) {
4771 for (i=0; i<natom; i++) {
4779 Vnm_print(0,
"forceAPOL: atom surf: Time elapsed: %f\n", ((
double)clock() - ts) / CLOCKS_PER_SEC);
4782 Vnm_print(0,
"forceAPOL: calcforce == ACF_TOTAL\n");
4786 if(*atomForce != VNULL){
4787 Vmem_free(mem,*nforce,
sizeof(
AtomForce), (
void **)atomForce);
4790 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4794 for (j=0; j<3; j++) {
4796 (*atomForce)[0].savForce[j] = 0.0;
4797 (*atomForce)[0].wcaForce[j] = 0.0;
4801 for (i=0; i<natom; i++) {
4810 if(VABS(gamma) > VSMALL) {
4813 if(VABS(press) > VSMALL) {
4816 if(VABS(bconc) > VSMALL) {
4821 (*atomForce)[0].sasaForce[j] += dSASA[j];
4822 (*atomForce)[0].savForce[j] += dSAV[j];
4823 (*atomForce)[0].wcaForce[j] += force[j];
4828 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A)\n");
4829 Vnm_tprint( 1,
" Legend:\n");
4830 Vnm_tprint( 1,
" sasa -- SASA force\n");
4831 Vnm_tprint( 1,
" sav -- SAV force\n");
4832 Vnm_tprint( 1,
" wca -- WCA force\n\n");
4834 Vnm_tprint( 1,
" sasa %4.3e %4.3e %4.3e\n",
4835 (*atomForce)[0].sasaForce[0],
4836 (*atomForce)[0].sasaForce[1],
4837 (*atomForce)[0].sasaForce[2]);
4838 Vnm_tprint( 1,
" sav %4.3e %4.3e %4.3e\n",
4839 (*atomForce)[0].savForce[0],
4840 (*atomForce)[0].savForce[1],
4841 (*atomForce)[0].savForce[2]);
4842 Vnm_tprint( 1,
" wca %4.3e %4.3e %4.3e\n",
4843 (*atomForce)[0].wcaForce[0],
4844 (*atomForce)[0].wcaForce[1],
4845 (*atomForce)[0].wcaForce[2]);
4847 Vnm_print(0,
"forceAPOL: calcforce == ACF_TOTAL: %f\n", ((
double)clock() - ts) / CLOCKS_PER_SEC);
4850 if(*atomForce == VNULL){
4851 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4854 Vmem_free(mem,*nforce,
sizeof(
AtomForce), (
void **)atomForce);
4855 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4860 Vnm_tprint( 1,
" Printing per atom forces (kJ/mol/A)\n");
4861 Vnm_tprint( 1,
" Legend:\n");
4862 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
4863 Vnm_tprint( 1,
" sasa n -- SASA force for atom n\n");
4864 Vnm_tprint( 1,
" sav n -- SAV force for atom n\n");
4865 Vnm_tprint( 1,
" wca n -- WCA force for atom n\n\n");
4867 Vnm_tprint( 1,
" gamma %f\n" \
4873 for (i=0; i<natom; i++) {
4883 for (j=0; j<3; j++) {
4884 (*atomForce)[i].sasaForce[j] = 0.0;
4885 (*atomForce)[i].savForce[j] = 0.0;
4886 (*atomForce)[i].wcaForce[j] = 0.0;
4889 if(VABS(gamma) > VSMALL)
Vacc_atomdSASA(acc, offset, srad, atom, dSASA);
4890 if(VABS(press) > VSMALL)
Vacc_atomdSAV(acc, srad, atom, dSAV);
4893 xF = -((gamma*dSASA[0]) + (press*dSAV[0]) + (bconc*force[0]));
4894 yF = -((gamma*dSASA[1]) + (press*dSAV[1]) + (bconc*force[1]));
4895 zF = -((gamma*dSASA[2]) + (press*dSAV[2]) + (bconc*force[2]));
4898 (*atomForce)[i].sasaForce[j] += dSASA[j];
4899 (*atomForce)[i].savForce[j] += dSAV[j];
4900 (*atomForce)[i].wcaForce[j] += force[j];
4904 Vnm_print( 1,
" tot %i %4.3e %4.3e %4.3e\n",
4909 Vnm_print( 1,
" sasa %i %4.3e %4.3e %4.3e\n",
4911 (*atomForce)[i].sasaForce[0],
4912 (*atomForce)[i].sasaForce[1],
4913 (*atomForce)[i].sasaForce[2]);
4914 Vnm_print( 1,
" sav %i %4.3e %4.3e %4.3e\n",
4916 (*atomForce)[i].savForce[0],
4917 (*atomForce)[i].savForce[1],
4918 (*atomForce)[i].savForce[2]);
4919 Vnm_print( 1,
" wca %i %4.3e %4.3e %4.3e\n",
4921 (*atomForce)[i].wcaForce[0],
4922 (*atomForce)[i].wcaForce[1],
4923 (*atomForce)[i].wcaForce[2]);
4933 Vnm_print(0,
"forceAPOL: Time elapsed: %f\n", ((
double)clock() - ts_main) / CLOCKS_PER_SEC);
4942VPUBLIC
int initBEM(
int icalc,
4964 Vnm_tprint(1,
"Destroying boundary element structures.\n");
4971void apbs2tabipb_(TABIPBparm* tabiparm,
4972 TABIPBvars* tabivars);
4985 if (nosh != VNULL) {
4986 if (nosh->
bogus)
return 1;
4991 TABIPBparm *tabiparm = (TABIPBparm*)calloc(1,
sizeof(TABIPBparm));
4993 sprintf(tabiparm->fpath,
"");
4994 strncpy(tabiparm->fname, nosh->
molpath[0],4);
4995 tabiparm->fname[4] =
'\0';
4996 tabiparm->density = pbeparm->
sdens;
4997 tabiparm->probe_radius = pbeparm->
srad;
4999 tabiparm->epsp = pbeparm->
pdie;
5000 tabiparm->epsw = pbeparm->
sdie;
5001 tabiparm->bulk_strength = 0.0;
5002 for (i=0; i<pbeparm->
nion; i++)
5003 tabiparm->bulk_strength += pbeparm->
ionc[i]
5004 *pbeparm->
ionq[i]*pbeparm->
ionq[i];
5005 tabiparm->temp = pbeparm->
temp;
5008 tabiparm->maxparnode = bemparm->
tree_n0;
5009 tabiparm->theta = bemparm->
mac;
5011 tabiparm->mesh_flag = bemparm->
mesh;
5015 tabiparm->output_datafile = bemparm->
outdata;
5017 TABIPBvars *tabivars = (TABIPBvars*)calloc(1,
sizeof(TABIPBvars));
5018 if ((tabivars->chrpos = (
double *) malloc(3 * tabiparm->number_of_lines *
sizeof(
double))) == NULL) {
5019 printf(
"Error in allocating t_chrpos!\n");
5021 if ((tabivars->atmchr = (
double *) malloc(tabiparm->number_of_lines *
sizeof(
double))) == NULL) {
5022 printf(
"Error in allocating t_atmchr!\n");
5024 if ((tabivars->atmrad = (
double *) malloc(tabiparm->number_of_lines *
sizeof(
double))) == NULL) {
5025 printf(
"Error in allocating t_atmrad!\n");
5030 for (i = 0; i < tabiparm->number_of_lines; i++){
5040 apbs2tabipb_(tabiparm, tabivars);
5043 free(tabivars->chrpos);
5044 free(tabivars->atmchr);
5045 free(tabivars->atmrad);
5046 free(tabivars->vert_ptl);
5047 free(tabivars->xvct);
5048 free_matrix(tabivars->vert);
5049 free_matrix(tabivars->snrm);
5050 free_matrix(tabivars->face);
5052 Vnm_tprint(1,
"\n\nReturning to APBS caller...\n\n");
5053 Vnm_tprint(1,
"Solvation energy and Coulombic energy in kJ/mol...\n\n");
5054 Vnm_tprint(1,
" Global net ELEC energy = %1.12E\n", tabivars->soleng);
5055 Vnm_tprint(1,
" Global net COULOMBIC energy = %1.12E\n\n", tabivars->couleng);
5065VPUBLIC
int setPartBEM(
NOsh *nosh,
5073 if (nosh->
bogus)
return 1;
5079VPUBLIC
int energyBEM(
NOsh *nosh,
5105VPUBLIC
int forceBEM(
5123 Vnm_tprint( 1,
" Calculating forces...\n");
5131VPUBLIC
void printBEMPARM(
BEMparm *bemparm) {
5136VPUBLIC
int writedataBEM(
int rank,
5145VPUBLIC
int writematBEM(
int rank,
NOsh *nosh,
PBEparm *pbeparm) {
5148 if (nosh->
bogus)
return 1;
5155#ifdef ENABLE_GEOFLOW
5167 if (nosh != VNULL) {
5168 if (nosh->
bogus)
return 1;
5173 struct GeometricFlowInput geoflowIn = getGeometricFlowParams();
5176 geoflowIn.m_boundaryCondition = (int) pbeparm->
bcfl ;
5177 geoflowIn.m_grid[0] = apolparm->
grid[0];
5178 geoflowIn.m_grid[1] = apolparm->
grid[1];
5179 geoflowIn.m_grid[2] = apolparm->
grid[2];
5180 geoflowIn.m_gamma = apolparm->
gamma;
5181 geoflowIn.m_pdie = pbeparm->
pdie ;
5182 geoflowIn.m_sdie = pbeparm->
sdie ;
5183 geoflowIn.m_press = apolparm->
press ;
5184 geoflowIn.m_tol = parm->
etol;
5185 geoflowIn.m_bconc = apolparm->
bconc ;
5186 geoflowIn.m_vdwdispersion = parm->vdw;
5187 geoflowIn.m_etolSolvation = .01 ;
5193 struct GeometricFlowOutput geoflowOut =
5194 runGeometricFlowWrapAPBS( geoflowIn, molecules[0] );
5196 Vnm_tprint( 1,
" Global net energy = %1.12E\n", geoflowOut.m_totalSolvation);
5197 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E\n", geoflowOut.m_elecSolvation);
5198 Vnm_tprint( 1,
" Global net APOL energy = %1.12E\n", geoflowOut.m_nonpolarSolvation);
5218 if (nosh != VNULL) {
5219 if (nosh->
bogus)
return 1;
5224 PBAMInput pbamIn = getPBAMParams();
5226 pbamIn.nmol_ = nosh->
nmol;
5229 pbamIn.temp_ = pbeparm->
temp;
5230 if (fabs(pbamIn.temp_-0.0) < 1e-3)
5232 printf(
"No temperature specified. Setting to 298.15K\n");
5233 pbamIn.temp_ = 298.15;
5237 pbamIn.idiel_ = pbeparm->
pdie;
5238 pbamIn.sdiel_ = pbeparm->
sdie;
5241 pbamIn.salt_ = parm->salt;
5244 strncpy(pbamIn.runType_, parm->runtype,
CHR_MAXLEN);
5245 strncpy(pbamIn.runName_, parm->runname,
CHR_MAXLEN);
5247 pbamIn.randOrient_ = parm->setrandorient;
5249 pbamIn.boxLen_ = parm->pbcboxlen;
5250 pbamIn.pbcType_ = parm->setpbcs;
5252 pbamIn.setunits_ = parm->setunits;
5253 if(parm->setunits == 1) strncpy(pbamIn.units_, parm->units,
CHR_MAXLEN);
5256 if (parm->setgridpt) pbamIn.gridPts_ = parm->gridpt;
5257 strncpy(pbamIn.map3D_, parm->map3dname,
CHR_MAXLEN);
5258 pbamIn.grid2Dct_ = parm->grid2Dct;
5259 for (i=0; i<pbamIn.grid2Dct_; i++)
5261 strncpy(pbamIn.grid2D_[i], parm->grid2Dname[i],
CHR_MAXLEN);
5262 strncpy(pbamIn.grid2Dax_[i], parm->grid2Dax[i],
CHR_MAXLEN);
5263 pbamIn.grid2Dloc_[i] = parm->grid2Dloc[i];
5265 strncpy(pbamIn.dxname_, parm->dxname,
CHR_MAXLEN);
5268 pbamIn.ntraj_ = parm->ntraj;
5269 strncpy(pbamIn.termCombine_, parm->termcombine,
CHR_MAXLEN);
5271 pbamIn.termct_ = parm->termct;
5272 pbamIn.contct_ = parm->confilct;
5274 if (strncmp(pbamIn.runType_,
"dynamics", 8)== 0)
5276 if (pbamIn.nmol_ > parm->diffct)
5278 Vnm_tprint(2,
"You need more diffusion information!\n");
5282 for (i=0; i<pbamIn.nmol_; i++)
5284 if (parm->xyzct[i] < parm->ntraj)
5286 Vnm_tprint(2,
"For molecule %d, you are missing trajectory!\n", i+1);
5289 for (j=0; j<pbamIn.ntraj_; j++)
5291 strncpy(pbamIn.xyzfil_[i][j], parm->xyzfil[i][j],
CHR_MAXLEN);
5296 for (i=0; i<pbamIn.nmol_; i++)
5298 strncpy(pbamIn.moveType_[i], parm->moveType[i],
CHR_MAXLEN);
5299 pbamIn.transDiff_[i] = parm->transDiff[i];
5300 pbamIn.rotDiff_[i] = parm->rotDiff[i];
5303 for (i=0; i<pbamIn.termct_; i++)
5305 strncpy(pbamIn.termnam_[i], parm->termnam[i],
CHR_MAXLEN);
5306 pbamIn.termnu_[i][0] = parm->termnu[i][0];
5307 pbamIn.termval_[i] = parm->termVal[i];
5310 for (i=0; i<pbamIn.contct_; i++)
5312 strncpy(pbamIn.confil_[i], parm->confil[i],
CHR_MAXLEN);
5318 printPBAMStruct( pbamIn );
5321 PBAMOutput pbamOut = runPBAMWrapAPBS( pbamIn, molecules, nosh->
nmol );
5323 Vnm_tprint(1,
"\n\nReturning to APBS caller...\n\n");
5325 if (!(strncmp(pbamIn.runType_,
"dynamics", 8) &&
5326 strncmp(pbamIn.runType_,
"energyforce", 11))) {
5328 if (!strncmp(pbamIn.units_,
"kcalmol", 7)) {
5330 Vnm_tprint(1,
"Interaction energy in kCal/mol...\n\n");
5332 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5333 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5334 i+1, pbamOut.energies_[i]);
5335 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5336 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5337 pbamOut.forces_[i][2]);
5338 if (pbamOut.energies_[i+1] == 0.)
break;
5341 }
else if (!strncmp(pbamIn.units_,
"jmol", 4)) {
5343 Vnm_tprint(1,
"Interaction energy in J/mol...\n\n");
5345 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5346 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5347 i+1, pbamOut.energies_[i]);
5348 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5349 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5350 pbamOut.forces_[i][2]);
5351 if (pbamOut.energies_[i+1] == 0.)
break;
5357 Vnm_tprint(1,
"Interaction energy in kT @ %6.2f K...\n\n", pbamIn.temp_);
5359 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5360 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5361 i+1, pbamOut.energies_[i]);
5362 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5363 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5364 pbamOut.forces_[i][2]);
5365 if (pbamOut.energies_[i+1] == 0.)
break;
5388 printf(
"solvePBSAM!!!\n");
5389 char fname_tp[VMAX_ARGLEN];
5390 if (nosh != VNULL) {
5391 if (nosh->
bogus)
return 1;
5396 PBSAMInput pbsamIn = getPBSAMParams();
5399 pbamIn.nmol_ = nosh->
nmol;
5402 pbamIn.temp_ = pbeparm->
temp;
5403 if (fabs(pbamIn.temp_-0.0) < 1e-3)
5405 printf(
"No temperature specified. Setting to 298.15K\n");
5406 pbamIn.temp_ = 298.15;
5410 pbamIn.idiel_ = pbeparm->
pdie;
5411 pbamIn.sdiel_ = pbeparm->
sdie;
5414 pbamIn.salt_ = parm->salt;
5417 strncpy(pbamIn.runType_, parm->runtype,
CHR_MAXLEN);
5418 strncpy(pbamIn.runName_, parm->runname,
CHR_MAXLEN);
5420 pbamIn.setunits_ = parm->setunits;
5421 if(parm->setunits == 1) strncpy(pbamIn.units_, parm->units,
CHR_MAXLEN);
5422 pbamIn.randOrient_ = parm->setrandorient;
5424 pbamIn.boxLen_ = parm->pbcboxlen;
5425 pbamIn.pbcType_ = parm->setpbcs;
5428 if (parm->setgridpt) pbamIn.gridPts_ = parm->gridpt;
5429 strncpy(pbamIn.map3D_, parm->map3dname,
CHR_MAXLEN);
5430 pbamIn.grid2Dct_ = parm->grid2Dct;
5431 for (i=0; i<pbamIn.grid2Dct_; i++)
5433 strncpy(pbamIn.grid2D_[i], parm->grid2Dname[i],
CHR_MAXLEN);
5434 strncpy(pbamIn.grid2Dax_[i], parm->grid2Dax[i],
CHR_MAXLEN);
5435 pbamIn.grid2Dloc_[i] = parm->grid2Dloc[i];
5437 strncpy(pbamIn.dxname_, parm->dxname,
CHR_MAXLEN);
5440 pbamIn.ntraj_ = parm->ntraj;
5441 strncpy(pbamIn.termCombine_, parm->termcombine,
CHR_MAXLEN);
5443 pbamIn.termct_ = parm->termct;
5444 pbamIn.contct_ = parm->confilct;
5446 if (strncmp(pbamIn.runType_,
"dynamics", 8)== 0)
5448 if (pbamIn.nmol_ > parm->diffct)
5450 Vnm_tprint(2,
"You need more diffusion information!\n");
5454 for (i=0; i<pbamIn.nmol_; i++)
5456 if (parm->xyzct[i] < parm->ntraj)
5458 Vnm_tprint(2,
"For molecule %d, you are missing trajectory!\n", i+1);
5461 for (j=0; j<pbamIn.ntraj_; j++)
5463 strncpy(pbamIn.xyzfil_[i][j], parm->xyzfil[i][j],
CHR_MAXLEN);
5468 for (i=0; i<pbamIn.nmol_; i++)
5470 strncpy(pbamIn.moveType_[i], parm->moveType[i],
CHR_MAXLEN);
5471 pbamIn.transDiff_[i] = parm->transDiff[i];
5472 pbamIn.rotDiff_[i] = parm->rotDiff[i];
5475 for (i=0; i<pbamIn.termct_; i++)
5477 strncpy(pbamIn.termnam_[i], parm->termnam[i],
CHR_MAXLEN);
5478 pbamIn.termnu_[i][0] = parm->termnu[i][0];
5479 pbamIn.termval_[i] = parm->termVal[i];
5482 for (i=0; i<pbamIn.contct_; i++)
5484 strncpy(pbamIn.confil_[i], parm->confil[i],
CHR_MAXLEN);
5490 pbsamIn.tolsp_ = samparm->tolsp;
5491 pbsamIn.imatct_ = samparm->imatct;
5492 pbsamIn.expct_ = samparm->expct;
5493 for (i=0; i<samparm->surfct; i++)
5495 strncpy(pbsamIn.surffil_[i], samparm->surffil[i],
CHR_MAXLEN);
5497 for (i=0; i<samparm->imatct; i++)
5499 strncpy(pbsamIn.imatfil_[i], samparm->imatfil[i],
CHR_MAXLEN);
5501 for (i=0; i<samparm->expct; i++)
5503 strncpy(pbsamIn.expfil_[i], samparm->expfil[i],
CHR_MAXLEN);
5507 if (samparm->setmsms == 1) {
5508 for (i=0; i<pbamIn.nmol_; i++) {
5510 for (j=0; j < VMAX_ARGLEN; j++)
5511 if (nosh->
molpath[i][j] ==
'\0')
break;
5515 char xyzr[VMAX_ARGLEN], surf[VMAX_ARGLEN], outname[VMAX_ARGLEN];
5516 for( k=0; k < j - 4; k++)
5518 xyzr[k] = nosh->
molpath[i][k];
5519 outname[k] = nosh->
molpath[i][k];
5520 surf[k] = nosh->
molpath[i][k];
5523 xyzr[k] =
'.'; surf[k] =
'.';
5524 xyzr[k+1] =
'x'; surf[k+1] =
'v';
5525 xyzr[k+2] =
'y'; surf[k+2] =
'e';
5526 xyzr[k+3] =
'z'; surf[k+3] =
'r';
5527 xyzr[k+4] =
'r'; surf[k+4] =
't';
5528 xyzr[k+5] =
'\0'; surf[k+5] =
'\0';;
5532 fp=fopen(xyzr,
"w");
5545 sprintf(fname_tp,
"msms.exe -if %s -prob %f -dens %f -of %s",
5546 xyzr, samparm->probe_radius,samparm->density, outname);
5548 sprintf(fname_tp,
"msms -if %s -prob %f -dens %f -of %s",
5549 xyzr, samparm->probe_radius,samparm->density, outname);
5552 printf(
"%s\n", fname_tp);
5554 printf(
"Running MSMS...\n");
5555 ierr = system(fname_tp);
5557 strncpy(pbsamIn.surffil_[i], surf,
CHR_MAXLEN);
5563 printPBSAMStruct( pbamIn, pbsamIn );
5566 PBAMOutput pbamOut = runPBSAMWrapAPBS(pbamIn, pbsamIn, molecules, nosh->
nmol);
5568 Vnm_tprint(1,
"\n\nReturning to APBS caller...\n\n");
5570 if (!(strncmp(pbamIn.runType_,
"dynamics", 8) &&
5571 strncmp(pbamIn.runType_,
"energyforce", 11))) {
5573 if (!strncmp(pbamIn.units_,
"kcalmol", 7)) {
5575 Vnm_tprint(1,
"Interaction energy in kCal/mol...\n\n");
5577 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5578 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5579 i+1, pbamOut.energies_[i]);
5580 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5581 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5582 pbamOut.forces_[i][2]);
5583 if (pbamOut.energies_[i+1] == 0.)
break;
5586 }
else if (!strncmp(pbamIn.units_,
"jmol", 4)) {
5588 Vnm_tprint(1,
"Interaction energy in J/mol...\n\n");
5590 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5591 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5592 i+1, pbamOut.energies_[i]);
5593 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5594 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5595 pbamOut.forces_[i][2]);
5596 if (pbamOut.energies_[i+1] == 0.)
break;
5602 Vnm_tprint(1,
"Interaction energy in kT @ %6.2f K...\n\n", pbamIn.temp_);
5604 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5605 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5606 i+1, pbamOut.energies_[i]);
5607 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5608 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5609 pbamOut.forces_[i][2]);
5610 if (pbamOut.energies_[i+1] == 0.)
break;
enum eBEMparm_CalcType BEMparm_CalcType
Declare BEMparm_CalcType type.
VPUBLIC void printPBEPARM(PBEparm *pbeparm)
Print out generic PBE params loaded from input.
VPUBLIC int loadPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the potential maps given in NOsh into grid objects.
VPUBLIC void killForce(Vmem *mem, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Free memory from MG force calculation.
VPUBLIC void storeAtomEnergy(Vpmg *pmg, int icalc, double **atomEnergy, int *nenergy)
Store energy in arrays for future use.
VPUBLIC int solveMG(NOsh *nosh, Vpmg *pmg, MGparm_CalcType type)
Solve the PBE with MG.
VPUBLIC int postRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Estimate error, mark mesh, and refine mesh after solve.
VPUBLIC int partFE(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Partition mesh (if applicable)
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 preRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Pre-refine mesh before solve.
VPUBLIC void startVio()
Wrapper to start MALOC Vio layer.
VPUBLIC void killKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded kappa maps.
VPUBLIC double returnEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Access net local energy.
VPUBLIC int printElecEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data.
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 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 int setPartMG(NOsh *nosh, MGparm *mgparm, Vpmg *pmg)
Set MG partitions for calculating observables and performing I/O.
VPUBLIC void killDielMaps(NOsh *nosh, Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL])
Destroy the loaded dielectric.
VPUBLIC void killMG(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC])
Kill structures initialized during an MG calculation.
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.
VPUBLIC int solveFE(int icalc, PBEparm *pbeparm, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Solve-estimate-refine.
VPUBLIC void killEnergy()
Kill arrays allocated for energies.
VPUBLIC int loadChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the charge maps given in NOsh into grid 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 void killMolecules(NOsh *nosh, Valist *alist[NOSH_MAXMOL])
Destroy the loaded molecules.
VPUBLIC void printFEPARM(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Print out FE-specific params loaded from input.
VPUBLIC int writedataMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out observables from MG calculation to file.
VPUBLIC int forceAPOL(Vacc *acc, Vmem *mem, APOLparm *apolparm, int *nforce, AtomForce **atomForce, Valist *alist, Vclist *clist)
Calculate non-polar forces.
VPUBLIC int printEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data (deprecated...see printElecEnergy)
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.
VPUBLIC int printApolForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data.
VPUBLIC int writematMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out operator matrix from MG calculation to file.
VPUBLIC int writedataFE(int rank, NOsh *nosh, PBEparm *pbeparm, Vfetk *fetk)
Write FEM data to files.
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 void printMGPARM(MGparm *mgparm, double realCenter[3])
Print out MG-specific params loaded from input.
VPUBLIC void killPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded potential maps.
VPUBLIC void killFE(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vfetk *fetk[NOSH_MAXCALC], Gem *gm[NOSH_MAXMOL])
Kill structures initialized during an FE calculation.
VPUBLIC int energyAPOL(APOLparm *apolparm, double sasa, double sav, double atomsasa[], double atomwcaEnergy[], int numatoms)
Calculate non-polar energies.
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 forceMG(Vmem *mem, NOsh *nosh, PBEparm *pbeparm, MGparm *mgparm, Vpmg *pmg, int *nforce, AtomForce **atomForce, Valist *alist[NOSH_MAXMOL])
Calculate forces from MG solution.
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 int printApolEnergy(NOsh *nosh, int iprint)
Combine and pretty-print energy data.
VPUBLIC int loadKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the kappa maps given in NOsh into grid objects.
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 Vparam * loadParameter(NOsh *nosh)
Loads and returns parameter object.
VPUBLIC void killChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded charge maps.
VPUBLIC int loadMolecules(NOsh *nosh, Vparam *param, Valist *alist[NOSH_MAXMOL])
Load the molecules given in NOsh into atom lists.
enum eMGparm_CalcType MGparm_CalcType
Declare MGparm_CalcType type.
#define NOSH_MAXCALC
Maximum number of calculations in a run.
#define NOSH_MAXMOL
Maximum number of molecules in a run.
#define CHR_MAXLEN
Number of things that can be written out in a single calculation.
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
int Vacc_wcaEnergyAtom(Vacc *thee, APOLparm *apolparm, Valist *alist, Vclist *clist, int iatom, double *value)
Calculate the WCA energy for an atom.
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius)
Return the total solvent accessible surface area (SASA)
VPUBLIC void Vacc_atomdSAV(Vacc *thee, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible volume.
VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom)
Return the atomic solvent accessible surface area (SASA)
VPUBLIC void Vacc_atomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible area.
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 Vacc_wcaForceAtom(Vacc *thee, APOLparm *apolparm, Vclist *clist, Vatom *atom, double *force)
Return the WCA integral force.
VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius)
Return the total solvent accessible volume (SAV)
VPUBLIC Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
VPUBLIC Vrc_Codes Valist_readPQR(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from a PQR file.
VPUBLIC Valist * Valist_ctor()
Construct the atom list object.
VPUBLIC void Valist_dtor(Valist **thee)
Destroys atom list object.
VPUBLIC Vrc_Codes Valist_readXML(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from an XML file.
VPUBLIC Vrc_Codes Valist_readPDB(Valist *thee, Vparam *param, Vio *sock)
Fill atom list with information from a PDB file.
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
VPUBLIC void Vclist_dtor(Vclist **thee)
Destroy object.
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.
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 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_dtor(Vfetk **thee)
Object destructor.
VPUBLIC Vfetk * Vfetk_ctor(Vpbe *pbe, Vhal_PBEType type)
Constructor for Vfetk object.
VPUBLIC double Vfetk_energy(Vfetk *thee, int color, int nonlin)
Return the total electrostatic energy.
VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type)
Fill an array with the specified data.
VPUBLIC void Vfetk_setAtomColors(Vfetk *thee)
Transfer color (partition ID) information frmo a partitioned mesh to the atoms.
VPUBLIC void Vfetk_setParameters(Vfetk *thee, PBEparm *pbeparm, FEMparm *feparm)
Set the parameter objects.
VPUBLIC void Vgrid_dtor(Vgrid **thee)
Object destructor.
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.
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 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 void Vgrid_writeDXBIN(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out the binary data in OpenDX grid format.
VPUBLIC int Vgrid_readGZ(Vgrid *thee, const char *fname)
Read in OpenDX data in GZIP format.
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.
VPUBLIC int Vgrid_readDXBIN(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read in binary data in OpenDX grid format.
#define APBS_TIMER_SETUP
APBS setup timer ID.
#define APBS_TIMER_SOLVER
APBS solver timer ID.
#define APBS_TIMER_ENERGY
APBS energy timer ID.
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
#define APBS_TIMER_FORCE
APBS force timer ID.
VPUBLIC Vparam * Vparam_ctor()
Construct the object.
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 Vparam_AtomData * Vparam_getAtomData(Vparam *thee, char resName[VMAX_ARGLEN], char atomName[VMAX_ARGLEN])
Get atom data.
VPUBLIC int Vparam_readXMLFile(Vparam *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read an XML format parameter database.
VPUBLIC void Vpbe_dtor(Vpbe **thee)
Object destructor.
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 double Vpbe_getDeblen(Vpbe *thee)
Get Debye-Huckel screening length.
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 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 void Vpmg_dtor(Vpmg **thee)
Object destructor.
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)
VPUBLIC double Vpmg_dielEnergy(Vpmg *thee, int extFlag)
Get the "polarization" contribution to the electrostatic energy.
VPUBLIC int Vpmg_solve(Vpmg *thee)
Solve the PBE using PMG.
VPUBLIC double Vpmg_qfAtomEnergy(Vpmg *thee, Vatom *atom)
Get the per-atom "fixed charge" contribution to the electrostatic energy.
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.
VPUBLIC double Vpmg_energy(Vpmg *thee, int extFlag)
Get the total electrostatic energy.
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 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 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 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 void Vpmgp_dtor(Vpmgp **thee)
Object destructor.
VPUBLIC Vpmgp * Vpmgp_ctor(MGparm *mgparm)
Construct PMG parameter object and initialize to default values.
VPUBLIC int Vstring_strcasecmp(const char *s1, const char *s2)
Case-insensitive string comparison (BSD standard)
#define Vunit_kb
Boltzmann constant.
#define Vunit_Na
Avogadro's number.
Header file for front end auxiliary routines.
Structure to hold atomic forces.
Reads and assigns charge/radii parameters.
Parameter structure for APOL-specific variables from input files.
APOLparm_calcEnergy calcenergy
APOLparm_calcForce calcforce
Parameter structure for BEM-specific variables from input files.
Parameter structure for FEM-specific variables from input files.
FEMparm_EstType akeySOLVE
Parameter structure for GEOFLOW-specific variables from input files.
Parameter structure for MG-specific variables from input files.
Class for parsing fixed format input files.
char apolname[NOSH_MAXCALC][VMAX_ARGLEN]
char chargepath[NOSH_MAXMOL][VMAX_ARGLEN]
char elecname[NOSH_MAXCALC][VMAX_ARGLEN]
char dielZpath[NOSH_MAXMOL][VMAX_ARGLEN]
Vdata_Format meshfmt[NOSH_MAXMOL]
char kappapath[NOSH_MAXMOL][VMAX_ARGLEN]
char dielXpath[NOSH_MAXMOL][VMAX_ARGLEN]
int printnarg[NOSH_MAXPRINT]
NOsh_calc * calc[NOSH_MAXCALC]
char meshpath[NOSH_MAXMOL][VMAX_ARGLEN]
char potpath[NOSH_MAXMOL][VMAX_ARGLEN]
Vdata_Format potfmt[NOSH_MAXMOL]
int printcalc[NOSH_MAXPRINT][NOSH_MAXPOP]
char molpath[NOSH_MAXMOL][VMAX_ARGLEN]
Vdata_Format kappafmt[NOSH_MAXMOL]
NOsh_MolFormat molfmt[NOSH_MAXMOL]
char parmpath[VMAX_ARGLEN]
int apol2calc[NOSH_MAXCALC]
Vdata_Format chargefmt[NOSH_MAXMOL]
NOsh_PrintType printwhat[NOSH_MAXPRINT]
int elec2calc[NOSH_MAXCALC]
Vdata_Format dielfmt[NOSH_MAXMOL]
char dielYpath[NOSH_MAXMOL][VMAX_ARGLEN]
int printop[NOSH_MAXPRINT][NOSH_MAXPOP]
Parameter structure for PBAM-specific variables from input files.
Parameter structure for PBE variables from input files.
PBEparm_calcEnergy calcenergy
char writematstem[VMAX_ARGLEN]
Vdata_Format writefmt[PBEPARM_MAXWRITE]
char writestem[PBEPARM_MAXWRITE][VMAX_ARGLEN]
Vdata_Type writetype[PBEPARM_MAXWRITE]
PBEparm_calcForce calcforce
Parameter structure for PBSAM-specific variables from input files.
Oracle for solvent- and ion-accessibility around a biomolecule.
Surface object list of per-atom surface points.
Container class for list of atom objects.
Contains public data members for Vatom class/module.
Contains public data members for Vfetk class/module.
Electrostatic potential oracle for Cartesian mesh data.
AtomData sub-class; stores atom data.
Contains public data members for Vpbe class/module.
Contains public data members for Vpmg class/module.
Contains public data members for Vpmgp class/module.
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.