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++) {
4795 (*atomForce)[0].sasaForce[j] = 0.0;
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;
struct sAPOLparm APOLparm
Declaration of the APOLparm class as the APOLparm structure.
enum eBEMparm_CalcType BEMparm_CalcType
Declare BEMparm_CalcType type.
struct sBEMparm BEMparm
Parameter structure for BEM-specific variables from input files.
struct sFEMparm FEMparm
Declaration of the FEMparm class as the FEMparm structure.
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.
struct sGEOFLOWparm GEOFLOWparm
Parameter structure for GEOFLOW-specific variables from input files.
struct sMGparm MGparm
Declaration of the MGparm class as the MGparm structure.
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.
struct sNOsh NOsh
Declaration of the NOsh class as the NOsh structure.
struct sPBAMparm PBAMparm
Parameter structure for PBAM-specific variables from input files.
#define CHR_MAXLEN
Number of things that can be written out in a single calculation.
struct sPBEparm PBEparm
Declaration of the PBEparm class as the PBEparm structure.
struct sPBSAMparm PBSAMparm
Parameter structure for PBSAM-specific variables from input files.
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.
struct sVaccSurf VaccSurf
Declaration of the VaccSurf class as the VaccSurf structure.
struct sVacc Vacc
Declaration of the Vacc class as the Vacc structure.
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.
struct sValist Valist
Declaration of the Valist class as the Valist structure.
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.
struct sVatom Vatom
Declaration of the Vatom class as the Vatom structure.
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
VPUBLIC void Vclist_dtor(Vclist **thee)
Destroy object.
struct sVclist Vclist
Declaration of the Vclist class as the Vclist structure.
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.
struct sVfetk Vfetk
Declaration of the Vfetk class as the Vfetk structure.
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.
struct sVgrid Vgrid
Declaration of the Vgrid class as the sVgrid structure.
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.
struct sVparam_AtomData Vparam_AtomData
Declaration of the Vparam_AtomData class as the sVparam_AtomData structure.
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.
struct sVpbe Vpbe
Declaration of the Vpbe class as the Vpbe structure.
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.
struct sVpmg Vpmg
Declaration of the Vpmg class as the Vpmg structure.
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...
struct sVpmgp Vpmgp
Declaration of the Vpmgp class as the sVpmgp structure.
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.
APOLparm_calcEnergy calcenergy
APOLparm_calcForce calcforce
FEMparm_EstType akeySOLVE
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]
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
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.