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++)
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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",
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",
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",
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",
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",
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",
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",
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 \ 1067 definition; no smoothing\n");
1068 Vnm_tprint( 1,
" Solvent probe radius: %4.3f A\n",
1072 Vnm_tprint( 1,
" Using \"molecular\" surface definition;\ 1073 harmonic average smoothing\n");
1074 Vnm_tprint( 1,
" Solvent probe radius: %4.3f A\n",
1078 Vnm_tprint( 1,
" Using spline-based surface definition;\ 1079 window = %4.3f\n", pbeparm->
swin);
1084 Vnm_tprint( 1,
" Temperature: %4.3f K\n", pbeparm->
temp);
1086 energies will be calculated\n");
1088 forces will be calculated \n");
1090 solvent 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 \ 1101 to be written to ");
1104 Vnm_tprint(1,
" Spline-based solvent accessibility \ 1105 to be written to ");
1108 Vnm_tprint(1,
" van der Waals solvent accessibility \ 1109 to be written to ");
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");
1321 pbe[icalc]->smsize = pbeparm->
smsize;
1322 pbe[icalc]->smvolume = pbeparm->
smvolume;
1323 pbe[icalc]->ipkey = pmgp[icalc]->ipkey;
1327 Vnm_tprint(2,
"Error! Unknown PBE type (%d)!\n", pbeparm->
pbetype);
1330 Vnm_tprint(0,
"Setting PDE center to local center...\n");
1331 pmgp[icalc]->bcfl = pbeparm->
bcfl;
1332 pmgp[icalc]->xcent = realCenter[0];
1333 pmgp[icalc]->ycent = realCenter[1];
1334 pmgp[icalc]->zcent = realCenter[2];
1338 Vnm_tprint( 2,
"Can't focus first calculation!\n");
1343 pmg[icalc] =
Vpmg_ctor(pmgp[icalc], pbe[icalc], 1, pmg[icalc-1],
1349 if (icalc>0)
Vpmg_dtor(&(pmg[icalc-1]));
1350 pmg[icalc] =
Vpmg_ctor(pmgp[icalc], pbe[icalc], 0, VNULL, mgparm,
PCE_NO);
1358 theDielXMap = dielXMap[pbeparm->
dielMapID-1];
1360 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1367 theDielYMap = dielYMap[pbeparm->
dielMapID-1];
1369 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1376 theDielZMap = dielZMap[pbeparm->
dielMapID-1];
1378 Vnm_print(2,
"Error! %d is not a valid dielectric map ID!\n",
1385 theKappaMap = kappaMap[pbeparm->
kappaMapID-1];
1387 Vnm_print(2,
"Error! %d is not a valid kappa map ID!\n",
1394 thePotMap = potMap[pbeparm->
potMapID-1];
1396 Vnm_print(2,
"Error! %d is not a valid potential map ID!\n",
1405 Vnm_print(2,
"Error! %d is not a valid charge map ID!\n",
1412 Vnm_print(2,
"Warning: You specified 'bcfl map' in the input file, but no potential map was found.\n");
1413 Vnm_print(2,
" You must specify 'usemap pot' statement in the APBS input file!\n");
1414 Vnm_print(2,
"Bailing out ...\n");
1427 Vnm_print(2,
"initMG: problems setting up coefficients (fillco)!\n");
1433 Vnm_tprint(1,
" Debye length: %g A\n",
Vpbe_getDeblen(pbe[icalc]));
1440 bytesTotal = Vmem_bytesTotal();
1441 highWater = Vmem_highWaterTotal();
1444 Vnm_tprint( 1,
" Current memory usage: %4.3f MB total, \ 1445 %4.3f MB high water\n", (
double)(bytesTotal)/(1024.*1024.),
1446 (
double)(highWater)/(1024.*1024.));
1459 Vnm_tprint(1,
"Destroying multigrid structures.\n");
1472 for(i=0;i<nosh->
ncalc;i++){
1489 if (nosh != VNULL) {
1490 if (nosh->
bogus)
return 1;
1498 Vnm_print(2,
" Error during PDE solution!\n");
1502 Vnm_tprint( 1,
" Skipping solve for mg-dummy run; zeroing \ 1507 for (i=0; i<
nx*
ny*
nz; i++) pmg->
u[i] = 0.0;
1524 if (nosh->
bogus)
return 1;
1527 for (j=0; j<3; j++) {
1532 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part center = (%g, %g, %g)\n",
1538 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part lower corner = (%g, %g, %g)\n",
1539 __FILE__, __LINE__, partMin[0], partMin[1], partMin[2]);
1540 Vnm_tprint(1,
"setPartMG (%s, %d): Disj part upper corner = (%g, %g, %g)\n",
1542 partMax[0], partMax[1], partMax[2]);
1545 for (j=0; j<3; j++) {
1546 partMin[j] = mgparm->
center[j] - 0.5*mgparm->
glen[j];
1547 partMax[j] = mgparm->
center[j] + 0.5*mgparm->
glen[j];
1588 if (nosh->
bogus == 0) {
1591 Vnm_tprint( 1,
" Total electrostatic energy = %1.12E kJ/mol\n",
1594 }
else *totEnergy = 0;
1602 Vnm_tprint( 1,
" Total electrostatic energy = %1.12E \ 1604 Vnm_tprint( 1,
" Fixed charge energy = %g kJ/mol\n",
1606 Vnm_tprint( 1,
" Mobile charge energy = %g kJ/mol\n",
1608 Vnm_tprint( 1,
" Dielectric energy = %g kJ/mol\n",
1610 Vnm_tprint( 1,
" Per-atom energies:\n");
1617 Vnm_tprint( 1,
" Atom %d: %1.12E kJ/mol\n", i,
1621 }
else *nenergy = 0;
1647 Vnm_tprint( 1,
" Calculating forces...\n");
1654 for (j=0; j<3; j++) {
1655 (*atomForce)[0].
qfForce[j] = 0;
1656 (*atomForce)[0].ibForce[j] = 0;
1657 (*atomForce)[0].dbForce[j] = 0;
1660 if (nosh->
bogus == 0) {
1665 for (k=0; k<3; k++) {
1671 for (k=0; k<3; k++) {
1672 (*atomForce)[0].qfForce[k] += qfForce[k];
1673 (*atomForce)[0].ibForce[k] += ibForce[k];
1674 (*atomForce)[0].dbForce[k] += dbForce[k];
1678 Vnm_tprint( 1,
" Printing net forces for molecule %d (kJ/mol/A)\n",
1680 Vnm_tprint( 1,
" Legend:\n");
1681 Vnm_tprint( 1,
" qf -- fixed charge force\n");
1682 Vnm_tprint( 1,
" db -- dielectric boundary force\n");
1683 Vnm_tprint( 1,
" ib -- ionic boundary force\n");
1684 Vnm_tprint( 1,
" qf %4.3e %4.3e %4.3e\n",
1688 Vnm_tprint( 1,
" ib %4.3e %4.3e %4.3e\n",
1692 Vnm_tprint( 1,
" db %4.3e %4.3e %4.3e\n",
1699 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
1702 Vnm_tprint( 1,
" Printing per-atom forces for molecule %d (kJ/mol/A)\n",
1704 Vnm_tprint( 1,
" Legend:\n");
1705 Vnm_tprint( 1,
" tot n -- total force for atom n\n");
1706 Vnm_tprint( 1,
" qf n -- fixed charge force for atom n\n");
1707 Vnm_tprint( 1,
" db n -- dielectric boundary force for atom n\n");
1708 Vnm_tprint( 1,
" ib n -- ionic boundary force for atom n\n");
1711 if (nosh->
bogus == 0) {
1719 for (k=0; k<3; k++) {
1720 (*atomForce)[j].qfForce[k] = 0;
1721 (*atomForce)[j].ibForce[k] = 0;
1722 (*atomForce)[j].dbForce[k] = 0;
1726 Vnm_tprint( 1,
"mgF tot %d %4.3e %4.3e %4.3e\n", j,
1728 *((*atomForce)[j].qfForce[0]+(*atomForce)[j].ibForce[0]+
1729 (*atomForce)[j].dbForce[0]),
1731 *((*atomForce)[j].qfForce[1]+(*atomForce)[j].ibForce[1]+
1732 (*atomForce)[j].dbForce[1]),
1734 *((*atomForce)[j].qfForce[2]+(*atomForce)[j].ibForce[2]+
1735 (*atomForce)[j].dbForce[2]));
1736 Vnm_tprint( 1,
"mgF qf %d %4.3e %4.3e %4.3e\n", j,
1738 *(*atomForce)[j].qfForce[0],
1740 *(*atomForce)[j].qfForce[1],
1742 *(*atomForce)[j].qfForce[2]);
1743 Vnm_tprint( 1,
"mgF ib %d %4.3e %4.3e %4.3e\n", j,
1745 *(*atomForce)[j].ibForce[0],
1747 *(*atomForce)[j].ibForce[1],
1749 *(*atomForce)[j].ibForce[2]);
1750 Vnm_tprint( 1,
"mgF db %d %4.3e %4.3e %4.3e\n", j,
1752 *(*atomForce)[j].dbForce[0],
1754 *(*atomForce)[j].dbForce[1],
1756 *(*atomForce)[j].dbForce[2]);
1769 Vnm_tprint(1,
"No energy arrays to destroy.\n");
1780 Vnm_tprint(1,
"Destroying force arrays.\n");
1783 for (i=0; i<nosh->
ncalc; i++) {
1785 if (nforce[i] > 0) Vmem_free(mem, nforce[i],
sizeof(
AtomForce),
1786 (
void **)&(atomForce[i]));
1793 char writematstem[VMAX_ARGLEN];
1794 char outpath[VMAX_ARGLEN];
1798 if (nosh->
bogus)
return 1;
1801 strlenmax = VMAX_ARGLEN-14;
1803 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1805 Vnm_tprint(2,
" Not writing matrix!\n");
1808 sprintf(writematstem,
"%s-PE%d", pbeparm->
writematstem, rank);
1810 strlenmax = (int)(VMAX_ARGLEN)-1;
1812 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1814 Vnm_tprint(2,
" Not writing matrix!\n");
1825 strlenmax = VMAX_ARGLEN-5;
1827 Vnm_tprint(2,
" Matrix name (%s) too long (%d char max)!\n",
1829 Vnm_tprint(2,
" Not writing matrix!\n");
1832 sprintf(outpath,
"%s.%s", writematstem,
"mat");
1838 Vnm_tprint( 1,
" Writing Poisson operator matrix \ 1839 to %s...\n", outpath);
1843 Vnm_tprint( 1,
" Writing linearization of full \ 1844 Poisson-Boltzmann operator matrix to %s...\n", outpath);
1847 Vnm_tprint( 2,
" Bogus matrix specification\ 1852 Vnm_tprint(0,
" Printing operator...\n");
1871 *atomEnergy = (
double *)Vmem_malloc(pmg->
vmem, *nenergy,
sizeof(
double));
1873 for (i=0; i<*nenergy; i++) {
1894 int ielec, icalc, i, j;
1895 char *timestring = VNULL;
1898 double conversion, ltenergy, gtenergy, scalar;
1900 if (nosh->
bogus)
return 1;
1906 file = fopen(fname,
"w");
1907 if (file == VNULL) {
1908 Vnm_print(2,
"writedataFlat: Problem opening virtual socket %s\n",
1916 timestring = ctime(&now);
1917 fprintf(file,
"%s\n", timestring);
1919 for (ielec=0; ielec<nosh->
nelec;ielec++) {
1929 fprintf(file,
"elec");
1931 fprintf(file,
" name %s\n", nosh->
elecname[ielec]);
1932 }
else fprintf(file,
"\n");
1934 switch (mgparm->
type) {
1936 fprintf(file,
" mg-dummy\n");
1939 fprintf(file,
" mg-manual\n");
1942 fprintf(file,
" mg-auto\n");
1945 fprintf(file,
" mg-para\n");
1951 fprintf(file,
" mol %d\n", pbeparm->
molid);
1952 fprintf(file,
" dime %d %d %d\n", mgparm->
dime[0], mgparm->
dime[1],\
1957 fprintf(file,
" npbe\n");
1960 fprintf(file,
" lpbe\n");
1966 if (pbeparm->
nion > 0) {
1967 for (i=0; i<pbeparm->
nion; i++) {
1968 fprintf(file,
" ion %4.3f %4.3f %4.3f\n",
1973 fprintf(file,
" pdie %4.3f\n", pbeparm->
pdie);
1974 fprintf(file,
" sdie %4.3f\n", pbeparm->
sdie);
1976 switch (pbeparm->
srfm) {
1978 fprintf(file,
" srfm mol\n");
1979 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1982 fprintf(file,
" srfm smol\n");
1983 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1986 fprintf(file,
" srfm spl2\n");
1987 fprintf(file,
" srad %4.3f\n", pbeparm->
srad);
1993 switch (pbeparm->
bcfl) {
1995 fprintf(file,
" bcfl zero\n");
1998 fprintf(file,
" bcfl sdh\n");
2001 fprintf(file,
" bcfl mdh\n");
2004 fprintf(file,
" bcfl focus\n");
2007 fprintf(file,
" bcfl map\n");
2010 fprintf(file,
" bcfl mem\n");
2016 fprintf(file,
" temp %4.3f\n", pbeparm->
temp);
2018 for (;icalc<=nosh->
elec2calc[ielec];icalc++){
2024 fprintf(file,
" calc\n");
2025 fprintf(file,
" id %i\n", (icalc+1));
2026 fprintf(file,
" grid %4.3f %4.3f %4.3f\n",
2028 fprintf(file,
" glen %4.3f %4.3f %4.3f\n",
2032 fprintf(file,
" totEnergy %1.12E kJ/mol\n",
2033 (totEnergy[icalc]*conversion));
2035 fprintf(file,
" totEnergy %1.12E kJ/mol\n",
2036 (totEnergy[icalc]*conversion));
2037 fprintf(file,
" qfEnergy %1.12E kJ/mol\n",
2038 (0.5*qfEnergy[icalc]*conversion));
2039 fprintf(file,
" qmEnergy %1.12E kJ/mol\n",
2040 (qmEnergy[icalc]*conversion));
2041 fprintf(file,
" dielEnergy %1.12E kJ/mol\n",
2042 (dielEnergy[icalc]*conversion));
2043 for (i=0; i<nenergy[icalc]; i++){
2044 fprintf(file,
" atom %i %1.12E kJ/mol\n", i,
2045 (0.5*atomEnergy[icalc][i]*conversion));
2051 fprintf(file,
" qfForce %1.12E %1.12E %1.12E kJ/mol/A\n",
2052 (atomForce[icalc][0].qfForce[0]*conversion),
2053 (atomForce[icalc][0].qfForce[1]*conversion),
2054 (atomForce[icalc][0].qfForce[2]*conversion));
2055 fprintf(file,
" ibForce %1.12E %1.12E %1.12E kJ/mol/A\n",
2056 (atomForce[icalc][0].ibForce[0]*conversion),
2057 (atomForce[icalc][0].ibForce[1]*conversion),
2058 (atomForce[icalc][0].ibForce[2]*conversion));
2059 fprintf(file,
" dbForce %1.12E %1.12E %1.12E kJ/mol/A\n",
2060 (atomForce[icalc][0].dbForce[0]*conversion),
2061 (atomForce[icalc][0].dbForce[1]*conversion),
2062 (atomForce[icalc][0].dbForce[2]*conversion));
2064 fprintf(file,
" end\n");
2067 fprintf(file,
"end\n");
2072 for (i=0; i<nosh->
nprint; i++) {
2076 fprintf(file,
"print energy");
2077 fprintf(file,
" %d", nosh->
printcalc[i][0]+1);
2080 if (nosh->
printop[i][j-1] == 0) fprintf(file,
" +");
2081 else if (nosh->
printop[i][j-1] == 1) fprintf(file,
" -");
2082 fprintf(file,
" %d", nosh->
printcalc[i][j]+1);
2085 fprintf(file,
"\n");
2094 if (nosh->
printop[i][j-1] == 0) scalar = 1.0;
2095 else if (nosh->
printop[i][j-1] == 1) scalar = -1.0;
2100 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2103 fprintf(file,
" localEnergy %1.12E kJ/mol\n", \
2105 fprintf(file,
" globalEnergy %1.12E kJ/mol\nend\n", \
2127 int ielec, icalc, i, j;
2128 char *timestring = VNULL;
2132 double conversion, ltenergy, gtenergy, scalar;
2134 if (nosh->
bogus)
return 1;
2140 file = fopen(fname,
"w");
2141 if (file == VNULL) {
2142 Vnm_print(2,
"writedataXML: Problem opening virtual socket %s\n",
2147 fprintf(file,
"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
2148 fprintf(file,
"<APBS>\n");
2153 timestring = ctime(&now);
2154 for(c = timestring; *c !=
'\n'; c++);
2156 fprintf(file,
" <date>%s</date>\n", timestring);
2158 for (ielec=0; ielec<nosh->
nelec;ielec++){
2168 fprintf(file,
" <elec>\n");
2170 fprintf(file,
" <name>%s</name>\n", nosh->
elecname[ielec]);
2173 switch (mgparm->
type) {
2175 fprintf(file,
" <type>mg-dummy</type>\n");
2178 fprintf(file,
" <type>mg-manual</type>\n");
2181 fprintf(file,
" <type>mg-auto</type>\n");
2184 fprintf(file,
" <type>mg-para</type>\n");
2190 fprintf(file,
" <molid>%d</molid>\n", pbeparm->
molid);
2191 fprintf(file,
" <nx>%d</nx>\n", mgparm->
dime[0]);
2192 fprintf(file,
" <ny>%d</ny>\n", mgparm->
dime[1]);
2193 fprintf(file,
" <nz>%d</nz>\n", mgparm->
dime[2]);
2197 fprintf(file,
" <pbe>npbe</pbe>\n");
2200 fprintf(file,
" <pbe>lpbe</pbe>\n");
2206 if (pbeparm->
nion > 0) {
2207 for (i=0; i<pbeparm->
nion; i++) {
2208 fprintf(file,
" <ion>\n");
2209 fprintf(file,
" <radius>%4.3f A</radius>\n",
2211 fprintf(file,
" <charge>%4.3f A</charge>\n",
2213 fprintf(file,
" <concentration>%4.3f M</concentration>\n",
2215 fprintf(file,
" </ion>\n");
2220 fprintf(file,
" <pdie>%4.3f</pdie>\n", pbeparm->
pdie);
2221 fprintf(file,
" <sdie>%4.3f</sdie>\n", pbeparm->
sdie);
2223 switch (pbeparm->
srfm) {
2225 fprintf(file,
" <srfm>mol</srfm>\n");
2226 fprintf(file,
" <srad>%4.3f</srad>\n", pbeparm->
srad);
2229 fprintf(file,
" <srfm>smol</srfm>\n");
2230 fprintf(file,
" <srad>%4.3f</srad>\n", pbeparm->
srad);
2233 fprintf(file,
" <srfm>spl2</srfm>\n");
2239 switch (pbeparm->
bcfl) {
2241 fprintf(file,
" <bcfl>zero</bcfl>\n");
2244 fprintf(file,
" <bcfl>sdh</bcfl>\n");
2247 fprintf(file,
" <bcfl>mdh</bcfl>\n");
2250 fprintf(file,
" <bcfl>focus</bcfl>\n");
2253 fprintf(file,
" <bcfl>map</bcfl>\n");
2256 fprintf(file,
" <bcfl>mem</bcfl>\n");
2262 fprintf(file,
" <temp>%4.3f K</temp>\n", pbeparm->
temp);
2264 for (;icalc<=nosh->
elec2calc[ielec];icalc++){
2270 fprintf(file,
" <calc>\n");
2271 fprintf(file,
" <id>%i</id>\n", (icalc+1));
2272 fprintf(file,
" <hx>%4.3f A</hx>\n", mgparm->
grid[0]);
2273 fprintf(file,
" <hy>%4.3f A</hy>\n", mgparm->
grid[1]);
2274 fprintf(file,
" <hz>%4.3f A</hz>\n", mgparm->
grid[2]);
2275 fprintf(file,
" <xlen>%4.3f A</xlen>\n", mgparm->
glen[0]);
2276 fprintf(file,
" <ylen>%4.3f A</ylen>\n", mgparm->
glen[1]);
2277 fprintf(file,
" <zlen>%4.3f A</zlen>\n", mgparm->
glen[2]);
2280 fprintf(file,
" <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2281 (totEnergy[icalc]*conversion));
2283 fprintf(file,
" <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2284 (totEnergy[icalc]*conversion));
2285 fprintf(file,
" <qfEnergy>%1.12E kJ/mol</qfEnergy>\n",
2286 (0.5*qfEnergy[icalc]*conversion));
2287 fprintf(file,
" <qmEnergy>%1.12E kJ/mol</qmEnergy>\n",
2288 (qmEnergy[icalc]*conversion));
2289 fprintf(file,
" <dielEnergy>%1.12E kJ/mol</dielEnergy>\n",
2290 (dielEnergy[icalc]*conversion));
2291 for (i=0; i<nenergy[icalc]; i++){
2292 fprintf(file,
" <atom>\n");
2293 fprintf(file,
" <id>%i</id>\n", i+1);
2294 fprintf(file,
" <energy>%1.12E kJ/mol</energy>\n",
2295 (0.5*atomEnergy[icalc][i]*conversion));
2296 fprintf(file,
" </atom>\n");
2302 fprintf(file,
" <qfforce_x>%1.12E</qfforce_x>\n",
2303 atomForce[icalc][0].qfForce[0]*conversion);
2304 fprintf(file,
" <qfforce_y>%1.12E</qfforce_y>\n",
2305 atomForce[icalc][0].qfForce[1]*conversion);
2306 fprintf(file,
" <qfforce_z>%1.12E</qfforce_z>\n",
2307 atomForce[icalc][0].qfForce[2]*conversion);
2308 fprintf(file,
" <ibforce_x>%1.12E</ibforce_x>\n",
2309 atomForce[icalc][0].ibForce[0]*conversion);
2310 fprintf(file,
" <ibforce_y>%1.12E</ibforce_y>\n",
2311 atomForce[icalc][0].ibForce[1]*conversion);
2312 fprintf(file,
" <ibforce_z>%1.12E</ibforce_z>\n",
2313 atomForce[icalc][0].ibForce[2]*conversion);
2314 fprintf(file,
" <dbforce_x>%1.12E</dbforce_x>\n",
2315 atomForce[icalc][0].dbForce[0]*conversion);
2316 fprintf(file,
" <dbforce_y>%1.12E</dbforce_y>\n",
2317 atomForce[icalc][0].dbForce[1]*conversion);
2318 fprintf(file,
" <dbforce_z>%1.12E</dbforce_z>\n",
2319 atomForce[icalc][0].dbForce[2]*conversion);
2322 fprintf(file,
" </calc>\n");
2325 fprintf(file,
" </elec>\n");
2330 for (i=0; i<nosh->
nprint; i++) {
2334 fprintf(file,
" <printEnergy>\n");
2335 fprintf(file,
" <equation>%d", nosh->
printcalc[i][0]+1);
2338 if (nosh->
printop[i][j-1] == 0) fprintf(file,
" +");
2339 else if (nosh->
printop[i][j-1] == 1) fprintf(file,
" -");
2340 fprintf(file,
" %d", nosh->
printcalc[i][j] +1);
2343 fprintf(file,
"</equation>\n");
2352 if (nosh->
printop[i][j-1] == 0) scalar = 1.0;
2353 else if (nosh->
printop[i][j-1] == 1) scalar = -1.0;
2358 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2359 fprintf(file,
" <localEnergy>%1.12E kJ/mol</localEnergy>\n", \
2361 fprintf(file,
" <globalEnergy>%1.12E kJ/mol</globalEnergy>\n", \
2364 fprintf(file,
" </printEnergy>\n");
2369 fprintf(file,
"</APBS>\n");
2381 char writestem[VMAX_ARGLEN];
2382 char outpath[VMAX_ARGLEN];
2402 if (nosh->
bogus)
return 1;
2404 for (i=0; i<pbeparm->
numwrite; i++) {
2417 Vnm_tprint(1,
" Writing charge distribution to ");
2426 sprintf(title,
"CHARGE DISTRIBUTION (e)");
2431 Vnm_tprint(1,
" Writing potential to ");
2440 sprintf(title,
"POTENTIAL (kT/e)");
2445 Vnm_tprint(1,
" Writing molecular accessibility to ");
2455 "SOLVENT ACCESSIBILITY -- MOLECULAR (%4.3f PROBE)",
2461 Vnm_tprint(1,
" Writing spline-based accessibility to ");
2471 "SOLVENT ACCESSIBILITY -- SPLINE (%4.3f WINDOW)",
2477 Vnm_tprint(1,
" Writing van der Waals accessibility to ");
2486 sprintf(title,
"SOLVENT ACCESSIBILITY -- VAN DER WAALS");
2491 Vnm_tprint(1,
" Writing ion accessibility to ");
2501 "ION ACCESSIBILITY -- SPLINE (%4.3f RADIUS)",
2507 Vnm_tprint(1,
" Writing potential Laplacian to ");
2517 "POTENTIAL LAPLACIAN (kT/e/A^2)");
2522 Vnm_tprint(1,
" Writing energy density to ");
2531 sprintf(title,
"ENERGY DENSITY (kT/e/A)^2");
2536 Vnm_tprint(1,
" Writing number density to ");
2546 "ION NUMBER DENSITY (M)");
2551 Vnm_tprint(1,
" Writing charge density to ");
2561 "ION CHARGE DENSITY (e_c * M)");
2566 Vnm_tprint(1,
" Writing x-shifted dielectric map to ");
2576 "X-SHIFTED DIELECTRIC MAP");
2581 Vnm_tprint(1,
" Writing y-shifted dielectric map to ");
2591 "Y-SHIFTED DIELECTRIC MAP");
2596 Vnm_tprint(1,
" Writing z-shifted dielectric map to ");
2606 "Z-SHIFTED DIELECTRIC MAP");
2611 Vnm_tprint(1,
" Writing kappa map to ");
2626 Vnm_tprint(1,
" Writing atom potentials to ");
2640 Vnm_tprint(2,
"Invalid data type for writing!\n");
2646 sprintf(writestem,
"%s-PE%d", pbeparm->
writestem[i], rank);
2651 sprintf(writestem,
"%s", pbeparm->
writestem[i]);
2658 sprintf(outpath,
"%s.%s", writestem,
"dx");
2659 Vnm_tprint(1,
"%s\n", outpath);
2668 sprintf(outpath,
"%s.%s", writestem,
"dxbin");
2669 Vnm_tprint(1,
"%s\n", outpath);
2679 sprintf(outpath,
"%s.%s", writestem,
"ucd");
2680 Vnm_tprint(1,
"%s\n", outpath);
2681 Vnm_tprint(2,
"Sorry, AVS format isn't supported for \ 2682 uniform meshes yet!\n");
2686 sprintf(outpath,
"%s.%s", writestem,
"mcsf");
2687 Vnm_tprint(1,
"%s\n", outpath);
2688 Vnm_tprint(2,
"Sorry, MCSF format isn't supported for \ 2689 uniform meshes yet!\n");
2693 sprintf(outpath,
"%s.%s", writestem,
"grd");
2694 Vnm_tprint(1,
"%s\n", outpath);
2703 sprintf(outpath,
"%s.%s", writestem,
"dx.gz");
2704 Vnm_tprint(1,
"%s\n", outpath);
2712 sprintf(outpath,
"%s.%s", writestem,
"txt");
2713 Vnm_tprint(1,
"%s\n", outpath);
2714 Vnm_print(0,
"routines: Opening virtual socket...\n");
2715 sock = Vio_ctor(
"FILE",
"ASC",VNULL,outpath,
"w");
2716 if (sock == VNULL) {
2717 Vnm_print(2,
"routines: Problem opening virtual socket %s\n",
2721 if (Vio_connect(sock, 0) < 0) {
2722 Vnm_print(2,
"routines: Problem connecting virtual socket %s\n",
2726 Vio_printf(sock,
"# Data from %s\n", PACKAGE_STRING);
2727 Vio_printf(sock,
"# \n");
2728 Vio_printf(sock,
"# %s\n", title);
2729 Vio_printf(sock,
"# \n");
2731 for(i=0;i<natoms;i++)
2732 Vio_printf(sock,
"%12.6e\n", pmg->
rwork[i]);
2735 Vnm_tprint(2,
"Bogus data format (%d)!\n",
2761 Vnm_tprint( 2,
" No energy available in Calculation %d\n", calcid+1);
2764 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++){
2767 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2768 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2789 Vnm_tprint( 2,
"Warning: The 'energy' print keyword is deprecated.\n" \
2790 " Use eilecEnergy for electrostatics energy calcs.\n\n");
2793 Vnm_tprint( 1,
"print energy %d ", nosh->
printcalc[iprint][0]+1);
2795 Vnm_tprint( 1,
"print energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2798 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2799 if (nosh->
printop[iprint][iarg-1] == 0)
2800 Vnm_tprint(1,
"+ ");
2801 else if (nosh->
printop[iprint][iarg-1] == 1)
2802 Vnm_tprint(1,
"- ");
2804 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2809 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2811 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2815 Vnm_tprint(1,
"end\n");
2821 Vnm_tprint( 2,
" Didn't calculate energy in Calculation \ 2825 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2828 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2829 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2835 Vnm_tprint( 1,
" Local net energy (PE %d) = %1.12E kJ/mol\n",
2836 Vcom_rank(com), ltenergy);
2837 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
2838 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2839 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2858 Vnm_tprint( 1,
"\nprint energy %d ", nosh->
printcalc[iprint][0]+1);
2860 Vnm_tprint( 1,
"\nprint energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2863 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2864 if (nosh->
printop[iprint][iarg-1] == 0)
2865 Vnm_tprint(1,
"+ ");
2866 else if (nosh->
printop[iprint][iarg-1] == 1)
2867 Vnm_tprint(1,
"- ");
2869 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2874 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2876 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2880 Vnm_tprint(1,
"end\n");
2886 Vnm_tprint( 2,
" Didn't calculate energy in Calculation \ 2890 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2893 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2894 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2900 Vnm_tprint( 1,
" Local net energy (PE %d) = %1.12E kJ/mol\n",
2901 Vcom_rank(com), ltenergy);
2902 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
2903 Vcom_reduce(com, <energy, >energy, 1, 2, 0);
2904 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2921 Vnm_tprint( 1,
"\nprint APOL energy %d ", nosh->
printcalc[iprint][0]+1);
2923 Vnm_tprint( 1,
"\nprint APOL energy %d (%s) ", nosh->
printcalc[iprint][0]+1,
2926 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2927 if (nosh->
printop[iprint][iarg-1] == 0)
2928 Vnm_tprint(1,
"+ ");
2929 else if (nosh->
printop[iprint][iarg-1] == 1)
2930 Vnm_tprint(1,
"- ");
2932 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
2937 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
2939 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
2943 Vnm_tprint(1,
"end\n");
2951 Vnm_tprint( 2,
" Didn't calculate energy in Calculation #%d\n", calcid+1);
2954 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
2959 if (nosh->
printop[iprint][iarg-1] == 0) scalar = 1.0;
2960 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
2962 gtenergy += (scalar * ((apolparm->
gamma*apolparm->
sasa) +
2968 Vnm_tprint( 1,
" Global net APOL energy = %1.12E kJ/mol\n", gtenergy);
2992 Vnm_tprint( 2,
"Warning: The 'force' print keyword is deprecated.\n" \
2993 " Use elecForce for electrostatics force calcs.\n\n");
2996 Vnm_tprint( 1,
"print force %d ", nosh->
printcalc[iprint][0]+1);
2998 Vnm_tprint( 1,
"print force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3001 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3002 if (nosh->
printop[iprint][iarg-1] == 0)
3003 Vnm_tprint(1,
"+ ");
3004 else if (nosh->
printop[iprint][iarg-1] == 1)
3005 Vnm_tprint(1,
"- ");
3007 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3012 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3014 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3018 Vnm_tprint(1,
"end\n");
3023 refnforce = nforce[calcid];
3025 if (refcalcforce ==
PCF_NO) {
3026 Vnm_tprint( 2,
" Didn't calculate force in calculation \ 3030 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3033 Vnm_tprint(2,
" Inconsistent calcforce declarations in \ 3038 if (nforce[calcid] != refnforce) {
3039 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \ 3052 aforce = atomForce[calcid];
3058 for (ivc=0; ivc<3; ivc++) {
3067 for (ifr=0; ifr<refnforce; ifr++) {
3068 for (ivc=0; ivc<3; ivc++) {
3080 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3083 aforce = atomForce[calcid];
3085 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3086 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3091 for (ivc=0; ivc<3; ivc++) {
3100 for (ifr=0; ifr<refnforce; ifr++) {
3101 for (ivc=0; ivc<3; ivc++) {
3113 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
3114 for (ifr=0; ifr<refnforce; ifr++) {
3115 Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
3116 Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
3117 Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
3122 Vnm_tprint( 1,
" Local net fixed charge force = \ 3123 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
3124 lforce[0].qfForce[1], lforce[0].qfForce[2]);
3125 Vnm_tprint( 1,
" Local net ionic boundary force = \ 3126 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
3127 lforce[0].ibForce[1], lforce[0].ibForce[2]);
3128 Vnm_tprint( 1,
" Local net dielectric boundary force = \ 3129 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
3130 lforce[0].dbForce[1], lforce[0].dbForce[2]);
3132 for (ifr=0; ifr<refnforce; ifr++) {
3133 Vnm_tprint( 1,
" Local fixed charge force \ 3134 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
3135 lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
3136 Vnm_tprint( 1,
" Local ionic boundary force \ 3137 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
3138 lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
3139 Vnm_tprint( 1,
" Local dielectric boundary force \ 3140 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
3141 lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
3147 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A).\n");
3148 Vnm_tprint( 1,
" Legend:\n");
3149 Vnm_tprint( 1,
" tot -- Total force\n");
3150 Vnm_tprint( 1,
" qf -- Fixed charge force\n");
3151 Vnm_tprint( 1,
" db -- Dielectric boundary force\n");
3152 Vnm_tprint( 1,
" ib -- Ionic boundary force\n");
3154 for (ivc=0; ivc<3; ivc++) {
3160 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3161 totforce[1], totforce[2]);
3162 Vnm_tprint( 1,
" qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3163 gforce[0].qfForce[1], gforce[0].qfForce[2]);
3164 Vnm_tprint( 1,
" ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3165 gforce[0].ibForce[1], gforce[0].ibForce[2]);
3166 Vnm_tprint( 1,
" db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3167 gforce[0].dbForce[1], gforce[0].dbForce[2]);
3171 Vnm_tprint( 1,
" Printing per-atom forces (kJ/mol/A).\n");
3172 Vnm_tprint( 1,
" Legend:\n");
3173 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3174 Vnm_tprint( 1,
" qf n -- Fixed charge force for atom n\n");
3175 Vnm_tprint( 1,
" db n -- Dielectric boundary force for atom n\n");
3176 Vnm_tprint( 1,
" ib n -- Ionic boundary force for atom n\n");
3177 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3183 for (ifr=0; ifr<refnforce; ifr++) {
3184 Vnm_tprint( 1,
" qf %d %1.12E %1.12E %1.12E\n", ifr,
3185 gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3186 gforce[ifr].qfForce[2]);
3187 Vnm_tprint( 1,
" ib %d %1.12E %1.12E %1.12E\n", ifr,
3188 gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3189 gforce[ifr].ibForce[2]);
3190 Vnm_tprint( 1,
" db %d %1.12E %1.12E %1.12E\n", ifr,
3191 gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3192 gforce[ifr].dbForce[2]);
3193 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3194 (gforce[ifr].dbForce[0] \
3195 + gforce[ifr].ibForce[0] +
3196 gforce[ifr].qfForce[0]),
3197 (gforce[ifr].dbForce[1] \
3198 + gforce[ifr].ibForce[1] +
3199 gforce[ifr].qfForce[1]),
3200 (gforce[ifr].dbForce[2] \
3201 + gforce[ifr].ibForce[2] +
3202 gforce[ifr].qfForce[2]));
3203 for (ivc=0; ivc<3; ivc++) {
3204 totforce[ivc] += (gforce[ifr].
dbForce[ivc] \
3209 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3210 totforce[1], totforce[2]);
3213 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3214 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3239 Vnm_tprint( 1,
"print force %d ", nosh->
printcalc[iprint][0]+1);
3241 Vnm_tprint( 1,
"print force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3244 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3245 if (nosh->
printop[iprint][iarg-1] == 0)
3246 Vnm_tprint(1,
"+ ");
3247 else if (nosh->
printop[iprint][iarg-1] == 1)
3248 Vnm_tprint(1,
"- ");
3250 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3255 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3257 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3261 Vnm_tprint(1,
"end\n");
3266 refnforce = nforce[calcid];
3268 if (refcalcforce ==
PCF_NO) {
3269 Vnm_tprint( 2,
" Didn't calculate force in calculation \ 3273 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3276 Vnm_tprint(2,
" Inconsistent calcforce declarations in \ 3281 if (nforce[calcid] != refnforce) {
3282 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \ 3295 aforce = atomForce[calcid];
3301 for (ivc=0; ivc<3; ivc++) {
3310 for (ifr=0; ifr<refnforce; ifr++) {
3311 for (ivc=0; ivc<3; ivc++) {
3323 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3326 aforce = atomForce[calcid];
3328 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3329 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3334 for (ivc=0; ivc<3; ivc++) {
3343 for (ifr=0; ifr<refnforce; ifr++) {
3344 for (ivc=0; ivc<3; ivc++) {
3356 Vnm_tprint( 0,
"printEnergy: Performing global reduction (sum)\n");
3357 for (ifr=0; ifr<refnforce; ifr++) {
3358 Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
3359 Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
3360 Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
3365 Vnm_tprint( 1,
" Local net fixed charge force = \ 3366 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
3367 lforce[0].qfForce[1], lforce[0].qfForce[2]);
3368 Vnm_tprint( 1,
" Local net ionic boundary force = \ 3369 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
3370 lforce[0].ibForce[1], lforce[0].ibForce[2]);
3371 Vnm_tprint( 1,
" Local net dielectric boundary force = \ 3372 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
3373 lforce[0].dbForce[1], lforce[0].dbForce[2]);
3375 for (ifr=0; ifr<refnforce; ifr++) {
3376 Vnm_tprint( 1,
" Local fixed charge force \ 3377 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
3378 lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
3379 Vnm_tprint( 1,
" Local ionic boundary force \ 3380 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
3381 lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
3382 Vnm_tprint( 1,
" Local dielectric boundary force \ 3383 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
3384 lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
3390 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A).\n");
3391 Vnm_tprint( 1,
" Legend:\n");
3392 Vnm_tprint( 1,
" tot -- Total force\n");
3393 Vnm_tprint( 1,
" qf -- Fixed charge force\n");
3394 Vnm_tprint( 1,
" db -- Dielectric boundary force\n");
3395 Vnm_tprint( 1,
" ib -- Ionic boundary force\n");
3397 for (ivc=0; ivc<3; ivc++) {
3403 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3404 totforce[1], totforce[2]);
3405 Vnm_tprint( 1,
" qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3406 gforce[0].qfForce[1], gforce[0].qfForce[2]);
3407 Vnm_tprint( 1,
" ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3408 gforce[0].ibForce[1], gforce[0].ibForce[2]);
3409 Vnm_tprint( 1,
" db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3410 gforce[0].dbForce[1], gforce[0].dbForce[2]);
3414 Vnm_tprint( 1,
" Printing per-atom forces (kJ/mol/A).\n");
3415 Vnm_tprint( 1,
" Legend:\n");
3416 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3417 Vnm_tprint( 1,
" qf n -- Fixed charge force for atom n\n");
3418 Vnm_tprint( 1,
" db n -- Dielectric boundary force for atom n\n");
3419 Vnm_tprint( 1,
" ib n -- Ionic boundary force for atom n\n");
3420 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3426 for (ifr=0; ifr<refnforce; ifr++) {
3427 Vnm_tprint( 1,
" qf %d %1.12E %1.12E %1.12E\n", ifr,
3428 gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3429 gforce[ifr].qfForce[2]);
3430 Vnm_tprint( 1,
" ib %d %1.12E %1.12E %1.12E\n", ifr,
3431 gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3432 gforce[ifr].ibForce[2]);
3433 Vnm_tprint( 1,
" db %d %1.12E %1.12E %1.12E\n", ifr,
3434 gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3435 gforce[ifr].dbForce[2]);
3436 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3437 (gforce[ifr].dbForce[0] \
3438 + gforce[ifr].ibForce[0] +
3439 gforce[ifr].qfForce[0]),
3440 (gforce[ifr].dbForce[1] \
3441 + gforce[ifr].ibForce[1] +
3442 gforce[ifr].qfForce[1]),
3443 (gforce[ifr].dbForce[2] \
3444 + gforce[ifr].ibForce[2] +
3445 gforce[ifr].qfForce[2]));
3446 for (ivc=0; ivc<3; ivc++) {
3447 totforce[ivc] += (gforce[ifr].
dbForce[ivc] \
3452 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3453 totforce[1], totforce[2]);
3456 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3457 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3484 Vnm_tprint( 1,
"\nprint APOL force %d ", nosh->
printcalc[iprint][0]+1);
3486 Vnm_tprint( 1,
"\nprint APOL force %d (%s) ", nosh->
printcalc[iprint][0]+1,
3489 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3490 if (nosh->
printop[iprint][iarg-1] == 0)
3491 Vnm_tprint(1,
"+ ");
3492 else if (nosh->
printop[iprint][iarg-1] == 1)
3493 Vnm_tprint(1,
"- ");
3495 Vnm_tprint( 2,
"Undefined PRINT operation!\n");
3500 Vnm_tprint( 1,
"%d ", nosh->
printcalc[iprint][iarg]+1);
3502 Vnm_tprint( 1,
"%d (%s) ", nosh->
printcalc[iprint][iarg]+1,
3506 Vnm_tprint(1,
"end\n");
3511 refnforce = nforce[calcid];
3513 if (refcalcforce ==
ACF_NO) {
3514 Vnm_tprint( 2,
" Didn't calculate force in calculation \ 3518 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3521 Vnm_tprint(2,
" Inconsistent calcforce declarations in \ 3526 if (nforce[calcid] != refnforce) {
3527 Vnm_tprint(2,
" Inconsistent number of forces evaluated in \ 3540 aforce = atomForce[calcid];
3546 for (ivc=0; ivc<3; ivc++) {
3552 for (ifr=0; ifr<refnforce; ifr++) {
3553 for (ivc=0; ivc<3; ivc++) {
3562 for (iarg=1; iarg<nosh->
printnarg[iprint]; iarg++) {
3565 aforce = atomForce[calcid];
3567 if (nosh->
printop[iprint][iarg-1] == 0) scalar = +1.0;
3568 else if (nosh->
printop[iprint][iarg-1] == 1) scalar = -1.0;
3573 for (ivc=0; ivc<3; ivc++) {
3579 for (ifr=0; ifr<refnforce; ifr++) {
3580 for (ivc=0; ivc<3; ivc++) {
3589 Vnm_tprint( 0,
"printForce: Performing global reduction (sum)\n");
3590 for (ifr=0; ifr<refnforce; ifr++) {
3591 Vcom_reduce(com, lforce[ifr].sasaForce, gforce[ifr].sasaForce, 3, 2, 0);
3592 Vcom_reduce(com, lforce[ifr].savForce, gforce[ifr].savForce, 3, 2, 0);
3593 Vcom_reduce(com, lforce[ifr].wcaForce, gforce[ifr].wcaForce, 3, 2, 0);
3597 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A)\n");
3598 Vnm_tprint( 1,
" Legend:\n");
3599 Vnm_tprint( 1,
" tot -- Total force\n");
3600 Vnm_tprint( 1,
" sasa -- SASA force\n");
3601 Vnm_tprint( 1,
" sav -- SAV force\n");
3602 Vnm_tprint( 1,
" wca -- WCA force\n\n");
3604 for (ivc=0; ivc<3; ivc++) {
3610 Vnm_tprint( 1,
" tot %1.12E %1.12E %1.12E\n", totforce[0],
3611 totforce[1], totforce[2]);
3612 Vnm_tprint( 1,
" sasa %1.12E %1.12E %1.12E\n", gforce[0].sasaForce[0],
3613 gforce[0].sasaForce[1], gforce[0].sasaForce[2]);
3614 Vnm_tprint( 1,
" sav %1.12E %1.12E %1.12E\n", gforce[0].savForce[0],
3615 gforce[0].savForce[1], gforce[0].savForce[2]);
3616 Vnm_tprint( 1,
" wca %1.12E %1.12E %1.12E\n", gforce[0].wcaForce[0],
3617 gforce[0].wcaForce[1], gforce[0].wcaForce[2]);
3621 Vnm_tprint( 1,
" Printing per atom forces (kJ/mol/A)\n");
3622 Vnm_tprint( 1,
" Legend:\n");
3623 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
3624 Vnm_tprint( 1,
" sasa n -- SASA force for atom n\n");
3625 Vnm_tprint( 1,
" sav n -- SAV force for atom n\n");
3626 Vnm_tprint( 1,
" wca n -- WCA force for atom n\n");
3627 Vnm_tprint( 1,
" tot all -- Total force for system\n");
3636 for (ifr=0; ifr<refnforce; ifr++) {
3637 Vnm_tprint( 1,
" sasa %d %1.12E %1.12E %1.12E\n", ifr,
3638 gforce[ifr].sasaForce[0], gforce[ifr].sasaForce[1],
3639 gforce[ifr].sasaForce[2]);
3640 Vnm_tprint( 1,
" sav %d %1.12E %1.12E %1.12E\n", ifr,
3641 gforce[ifr].savForce[0], gforce[ifr].savForce[1],
3642 gforce[ifr].savForce[2]);
3643 Vnm_tprint( 1,
" wca %d %1.12E %1.12E %1.12E\n", ifr,
3644 gforce[ifr].wcaForce[0], gforce[ifr].wcaForce[1],
3645 gforce[ifr].wcaForce[2]);
3646 Vnm_tprint( 1,
" tot %d %1.12E %1.12E %1.12E\n", ifr,
3647 (gforce[ifr].wcaForce[0] \
3648 + gforce[ifr].savForce[0] +
3649 gforce[ifr].sasaForce[0]),
3650 (gforce[ifr].wcaForce[1] \
3651 + gforce[ifr].savForce[1] +
3652 gforce[ifr].sasaForce[1]),
3653 (gforce[ifr].wcaForce[2] \
3654 + gforce[ifr].savForce[2] +
3655 gforce[ifr].sasaForce[2]));
3656 for (ivc=0; ivc<3; ivc++) {
3657 totforce[ivc] += (gforce[ifr].
wcaForce[ivc] \
3662 Vnm_tprint( 1,
" tot all %1.12E %1.12E %1.12E\n", totforce[0],
3663 totforce[1], totforce[2]);
3666 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&lforce));
3667 Vmem_free(VNULL, refnforce,
sizeof(
AtomForce), (
void **)(&gforce));
3684 Vnm_tprint(1,
"Destroying finite element structures.\n");
3687 for(i=0;i<nosh->
ncalc;i++){
3691 for (i=0; i<nosh->
nmesh; i++) {
3729 Vatom *atom = VNULL;
3731 Vnm_tstart(27,
"Setup timer");
3734 if (pbeparm->
useDielMap) Vnm_tprint(2,
"FEM ignoring dielectric map!\n");
3735 if (pbeparm->
useKappaMap) Vnm_tprint(2,
"FEM ignoring kappa map!\n");
3736 if (pbeparm->
useChargeMap) Vnm_tprint(2,
"FEM ignoring charge map!\n");
3739 Vnm_tprint(0,
"Re-centering mesh...\n");
3740 theMol = pbeparm->
molid-1;
3741 myalist = alist[theMol];
3742 for (j=0; j<3; j++) {
3743 if (theMol < nosh->nmol) {
3744 center[j] = (myalist)->center[j];
3746 Vnm_tprint(2,
"ERROR! Bogus molecule number (%d)!\n",
3774 Vnm_tprint(0,
"Setting up PBE object...\n");
3776 sparm = pbeparm->
swin;
3779 sparm = pbeparm->
srad;
3781 if (pbeparm->
nion > 0) {
3782 iparm = pbeparm->
ionr[0];
3788 pbeparm->
sdie, sparm, focusFlag, pbeparm->
sdens,
3793 Vnm_tprint(1,
" Debye length: %g A\n",
Vpbe_getDeblen(pbe[icalc]));
3796 Vnm_tprint(0,
"Setting up FEtk object...\n");
3803 Vnm_tprint(0,
"Setting up mesh...\n");
3806 imesh = feparm->
meshID-1;
3808 for (i=0; i<3; i++) {
3812 Vnm_print(0,
"Using mesh %d (%s) in calculation.\n", imesh+1,
3814 switch (nosh->
meshfmt[imesh]) {
3816 Vnm_tprint(2,
"DX finite element mesh input not supported yet!\n");
3819 Vnm_tprint(2,
"DXBIN finite element mesh input not supported yet!\n");
3822 Vnm_tprint( 2,
"UHBD finite element mesh input not supported!\n");
3825 Vnm_tprint( 2,
"AVS finite element mesh input not supported!\n");
3828 Vnm_tprint(1,
"Reading MCSF-format input finite element mesh from %s.\n",
3830 sock = Vio_ctor(
"FILE",
"ASC", VNULL, nosh->
meshpath[imesh],
"r");
3831 if (sock == VNULL) {
3832 Vnm_print(2,
"Problem opening virtual socket %s!\n",
3836 if (Vio_accept(sock, 0) < 0) {
3837 Vnm_print(2,
"Problem accepting virtual socket %s!\n",
3844 Vnm_tprint( 2,
"Invalid data format (%d)!\n",
3850 for (i=0; i<3; i++) {
3851 center[i] = alist[theMol]->center[i];
3852 length[i] = feparm->
glen[i];
3857 vrc =
Vfetk_loadMesh(fetk[icalc], center, length, meshType, sock);
3859 Vnm_print(2,
"Error constructing finite element mesh!\n");
3863 Gem_shapeChk(fetk[icalc]->gm);
3867 for (j=0; j<2; j++) {
3871 AM_markRefine(fetk[icalc]->am, 0, -1, 0, 0.0);
3873 AM_refine(fetk[icalc]->am, 2, 0, feparm->
pkey);
3875 Gem_shapeChk(fetk[icalc]->gm);
3880 Vnm_tstop(27,
"Setup timer");
3883 bytesTotal = Vmem_bytesTotal();
3884 highWater = Vmem_highWaterTotal();
3887 Vnm_tprint( 1,
" Current memory usage: %4.3f MB total, \ 3888 %4.3f MB high water\n", (
double)(bytesTotal)/(1024.*1024.),
3889 (
double)(highWater)/(1024.*1024.));
3901 Vnm_tprint(1,
" Domain size: %g A x %g A x %g A\n",
3904 switch(feparm->
ekey) {
3906 Vnm_tprint(1,
" Per-simplex error tolerance: %g\n", feparm->
etol);
3909 Vnm_tprint(1,
" Global error tolerance: %g\n", feparm->
etol);
3912 Vnm_tprint(1,
" Fraction of simps to refine: %g\n", feparm->
etol);
3915 Vnm_tprint(2,
"Invalid ekey (%d)!\n", feparm->
ekey);
3921 Vnm_tprint(1,
" Uniform pre-solve refinement.\n");
3924 Vnm_tprint(1,
" Geometry-based pre-solve refinement.\n");
3927 Vnm_tprint(2,
"Invalid akeyPRE (%d)!\n", feparm->
akeyPRE);
3933 Vnm_tprint(1,
" Residual-based a posteriori refinement.\n");
3936 Vnm_tprint(1,
" Dual-based a posteriori refinement.\n");
3939 Vnm_tprint(1,
" Local-based a posteriori refinement.\n");
3942 Vnm_tprint(2,
"Invalid akeySOLVE (%d)!\n", feparm->
akeySOLVE);
3945 Vnm_tprint(1,
" Refinement of initial mesh to ~%d vertices\n",
3947 Vnm_tprint(1,
" Geometry-based refinment lower bound: %g A\n",
3949 Vnm_tprint(1,
" Maximum number of solve-estimate-refine cycles: %d\n",
3951 Vnm_tprint(1,
" Maximum number of vertices in mesh: %d\n",
3955 if (nosh->
bogus)
return;
3957 Vnm_tprint(1,
" HB linear solver: AM_hPcg\n");
3959 Vnm_tprint(1,
" Non-HB linear solver: ");
3960 switch (fetk[icalc]->lkey) {
3962 Vnm_print(1,
"SLU direct\n");
3965 Vnm_print(1,
"multigrid\n");
3968 Vnm_print(1,
"conjugate gradient\n");
3971 Vnm_print(1,
"BiCGStab\n");
3974 Vnm_print(1,
"???\n");
3979 Vnm_tprint(1,
" Linear solver tol.: %g\n", fetk[icalc]->ltol);
3980 Vnm_tprint(1,
" Linear solver max. iters.: %d\n", fetk[icalc]->lmax);
3981 Vnm_tprint(1,
" Linear solver preconditioner: ");
3982 switch (fetk[icalc]->lprec) {
3984 Vnm_print(1,
"identity\n");
3987 Vnm_print(1,
"diagonal\n");
3990 Vnm_print(1,
"multigrid\n");
3993 Vnm_print(1,
"???\n");
3996 Vnm_tprint(1,
" Nonlinear solver: ");
3997 switch (fetk[icalc]->nkey) {
3999 Vnm_print(1,
"newton\n");
4002 Vnm_print(1,
"incremental\n");
4005 Vnm_print(1,
"pseudo-arclength\n");
4008 Vnm_print(1,
"??? ");
4011 Vnm_tprint(1,
" Nonlinear solver tol.: %g\n", fetk[icalc]->ntol);
4012 Vnm_tprint(1,
" Nonlinear solver max. iters.: %d\n", fetk[icalc]->nmax);
4013 Vnm_tprint(1,
" Initial guess: ");
4014 switch (fetk[icalc]->gues) {
4016 Vnm_tprint(1,
"zero\n");
4019 Vnm_tprint(1,
"boundary function\n");
4022 Vnm_tprint(1,
"interpolated previous solution\n");
4025 Vnm_tprint(1,
"???\n");
4050 Vnm_tprint(1,
" Commencing uniform refinement to %d verts.\n",
4054 Vnm_tprint(1,
" Commencing geometry-based refinement to %d \ 4058 Vnm_tprint(2,
"What? You can't do a posteriori error estimation \ 4059 before you solve!\n");
4074 Vnm_tprint(1,
" Initial mesh has %d vertices\n",
4075 Gem_numVV(fetk[icalc]->gm));
4081 nverts = Gem_numVV(fetk[icalc]->gm);
4083 Vnm_tprint(1,
" Hit vertex number limit.\n");
4086 marked = AM_markRefine(fetk[icalc]->am, feparm->
akeyPRE, -1,
4089 Vnm_tprint(1,
" Marked 0 simps; hit error/size tolerance.\n");
4092 Vnm_tprint(1,
" Have %d verts, marked %d. Refining...\n", nverts,
4094 AM_refine(fetk[icalc]->am, 0, 0, feparm->
pkey);
4097 nverts = Gem_numVV(fetk[icalc]->gm);
4098 Vnm_tprint(1,
" Done refining; have %d verts.\n", nverts);
4120 (pbeparm->
pbetype == PBE_NRPBE) ||
4143 Vnm_print(2,
"SORRY! DON'T USE HB!!!\n");
4147 AM_hlSolve(fetk[icalc]->am,
meth, lkeyHB, fetk[icalc]->lmax,
4148 fetk[icalc]->ltol, fetk[icalc]->gues, fetk[icalc]->pjac);
4193 if (nosh->
bogus == 0) {
4195 (pbeparm->
pbetype==PBE_NRPBE) ||
4206 Vnm_tprint(1,
" Total electrostatic energy = %1.12E kJ/mol\n",
4213 Vnm_tprint(2,
"Error! Verbose energy evaluation not available for FEM yet!\n");
4214 Vnm_tprint(2,
"E-mail nathan.baker@pnl.gov if you want this.\n");
4239 nverts = Gem_numVV(fetk[icalc]->gm);
4240 if (nverts > feparm->
maxvert) {
4241 Vnm_tprint(1,
" Current number of vertices (%d) exceeds max (%d)!\n",
4245 Vnm_tprint(1,
" Mesh currently has %d vertices\n", nverts);
4249 Vnm_tprint(1,
" Commencing uniform refinement.\n");
4252 Vnm_tprint(1,
" Commencing geometry-based refinement.\n");
4255 Vnm_tprint(1,
" Commencing residual-based refinement.\n");
4258 Vnm_tprint(1,
" Commencing dual problem-based refinement.\n");
4261 Vnm_tprint(1,
" Commencing local-based refinement.\n.");
4264 Vnm_tprint(2,
" Error -- unknown refinement type (%d)!\n",
4271 marked = AM_markRefine(fetk[icalc]->am, feparm->
akeySOLVE, -1,
4274 Vnm_tprint(1,
" Marked 0 simps; hit error/size tolerance.\n");
4277 Vnm_tprint(1,
" Have %d verts, marked %d. Refining...\n", nverts,
4279 AM_refine(fetk[icalc]->am, 0, 0, feparm->
pkey);
4280 nverts = Gem_numVV(fetk[icalc]->gm);
4281 Vnm_tprint(1,
" Done refining; have %d verts.\n", nverts);
4283 Gem_shapeChk(fetk[icalc]->gm);
4298 char writestem[VMAX_ARGLEN];
4299 char outpath[VMAX_ARGLEN];
4305 if (nosh->
bogus)
return 1;
4310 for (i=0; i<pbeparm->
numwrite; i++) {
4318 Vnm_tprint(2,
" Sorry; can't write charge distribution for FEM!\n");
4324 Vnm_tprint(1,
" Writing potential to ");
4330 Vnm_tprint(1,
" Writing molecular accessibility to ");
4336 Vnm_tprint(1,
" Writing spline-based accessibility to ");
4342 Vnm_tprint(1,
" Writing van der Waals accessibility to ");
4348 Vnm_tprint(1,
" Writing ion accessibility to ");
4354 Vnm_tprint(2,
" Sorry; can't write charge distribution for FEM!\n");
4360 Vnm_tprint(2,
" Sorry; can't write energy density for FEM!\n");
4366 Vnm_tprint(1,
" Writing number density to ");
4372 Vnm_tprint(1,
" Writing charge density to ");
4378 Vnm_tprint(2,
" Sorry; can't write x-shifted dielectric map for FEM!\n");
4384 Vnm_tprint(2,
" Sorry; can't write y-shifted dielectric map for FEM!\n");
4390 Vnm_tprint(2,
" Sorry; can't write z-shifted dielectric map for FEM!\n");
4396 Vnm_tprint(1,
" Sorry; can't write kappa map for FEM!\n");
4402 Vnm_tprint(1,
" Sorry; can't write atom potentials for FEM!\n");
4408 Vnm_tprint(2,
"Invalid data type for writing!\n");
4413 if (!writeit)
return 0;
4417 sprintf(writestem,
"%s-PE%d", pbeparm->
writestem[i], rank);
4422 sprintf(writestem,
"%s", pbeparm->
writestem[i]);
4429 sprintf(outpath,
"%s.%s", writestem,
"dx");
4430 Vnm_tprint(1,
"%s\n", outpath);
4436 sprintf(outpath,
"%s.%s", writestem,
"dxbin");
4437 Vnm_tprint(1,
"%s\n", outpath);
4442 sprintf(outpath,
"%s.%s", writestem,
"ucd");
4443 Vnm_tprint(1,
"%s\n", outpath);
4448 Vnm_tprint(2,
"UHBD format not supported for FEM!\n");
4452 Vnm_tprint(2,
"MCSF format not supported yet!\n");
4456 Vnm_tprint(2,
"Bogus data format (%d)!\n",
4484 Vatom *atom = VNULL;
4527 for (i=0; i < natoms; i++) {
4533 if ((x+atomRadius) >
xmax)
xmax = x + atomRadius;
4534 if ((x-atomRadius) <
xmin)
xmin = x - atomRadius;
4535 if ((y+atomRadius) >
ymax)
ymax = y + atomRadius;
4536 if ((y-atomRadius) <
ymin)
ymin = y - atomRadius;
4537 if ((z+atomRadius) >
zmax)
zmax = z + atomRadius;
4538 if ((z-atomRadius) <
zmin)
zmin = z - atomRadius;
4539 disp[0] = (x - center[0]);
4540 disp[1] = (y - center[1]);
4541 disp[2] = (z - center[2]);
4542 dist = (disp[0]*disp[0]) + (disp[1]*disp[1]) + (disp[2]*disp[2]);
4543 dist = VSQRT(dist) + atomRadius;
4551 Vnm_print(0,
"APOL: Setting up hash table and accessibility object...\n");
4552 nhash[0] = soluteXlen/0.5;
4553 nhash[1] = soluteYlen/0.5;
4554 nhash[2] = soluteZlen/0.5;
4555 for (i=0; i<3; i++) inhash[i] = (
int)(nhash[i]);
4558 if (inhash[i] < 3) inhash[i] = 3;
4559 if (inhash[i] > MAX_HASH_DIM) inhash[i] = MAX_HASH_DIM;
4563 srad = apolparm->
srad;
4564 sradPad = srad + (2*apolparm->
dpos);
4570 if (param == VNULL && (apolparm->
bconc != 0.0)) {
4571 Vnm_tprint(2,
"initAPOL: Got NULL Vparam object!\n");
4572 Vnm_tprint(2,
"initAPOL: You are performing an apolar calculation with the van der Waals integral term,\n");
4573 Vnm_tprint(2,
"initAPOL: this term requires van der Waals parameters which are not available from the \n");
4574 Vnm_tprint(2,
"initAPOL: PQR file. Therefore, you need to supply a parameter file with the parm keyword,\n");
4575 Vnm_tprint(2,
"initAPOL: for example,\n");
4576 Vnm_tprint(2,
"initAPOL: read parm flat amber94.dat end\n");
4577 Vnm_tprint(2,
"initAPOL: where the relevant parameter files can be found in apbs/tools/conversion/param/vparam.\n");
4581 if (apolparm->
bconc != 0.0){
4584 if (atomData == VNULL) {
4585 Vnm_tprint(2,
"initAPOL: Couldn't find parameters for WAT OW or WAT O!\n");
4586 Vnm_tprint(2,
"initAPOL: These parameters must be present in your file\n");
4587 Vnm_tprint(2,
"initAPOL: for apolar calculations.\n");
4597 rc =
forceAPOL(acc, mem, apolparm, nforce, atomForce, alist, clist);
4599 Vnm_print(2,
"Error in apolar force calculation!\n");
4611 if (VABS(apolparm->
gamma) > VSMALL) {
4615 for (i = 0; i < len; i++) {
4621 apolparm->
sasa = 0.0;
4623 for (i = 0; i < len; i++) {
4629 if (VABS(apolparm->
press) > VSMALL){
4632 apolparm->
sav = 0.0;
4636 if (VABS(apolparm->
bconc) > VSMALL) {
4638 for (i = 0; i < len; i++) {
4641 Vnm_print(2,
"Error in apolar energy calculation!\n");
4644 atomwcaEnergy[i] = energy;
4649 Vnm_print(2,
"Error in apolar energy calculation!\n");
4670 double atomwcaEnergy[],
4674 double energy = 0.0;
4678 Vnm_print(1,
"\nSolvent Accessible Surface Area (SASA) for each atom:\n");
4679 for (i = 0; i < numatoms; i++) {
4680 Vnm_print(1,
" SASA for atom %i: %1.12E\n", i, atomsasa[i]);
4683 Vnm_print(1,
"\nTotal solvent accessible surface area: %g A^2\n",sasa);
4690 Vnm_print(1,
"energyAPOL: Cannot calculate component energy, skipping.\n");
4693 energy = (apolparm->
gamma*sasa) + (apolparm->
press*sav)
4696 Vnm_print(1,
"\nSurface tension*area energies (gamma * SASA) for each atom:\n");
4697 for (i = 0; i < numatoms; i++) {
4698 Vnm_print(1,
" Surface tension*area energy for atom %i: %1.12E\n", i, apolparm->
gamma*atomsasa[i]);
4701 Vnm_print(1,
"\nTotal surface tension energy: %g kJ/mol\n", apolparm->
gamma*sasa);
4702 Vnm_print(1,
"\nTotal solvent accessible volume: %g A^3\n", sav);
4703 Vnm_print(1,
"\nTotal pressure*volume energy: %g kJ/mol\n", apolparm->
press*sav);
4704 Vnm_print(1,
"\nWCA dispersion Energies for each atom:\n");
4705 for (i = 0; i < numatoms; i++) {
4706 Vnm_print(1,
" WCA energy for atom %i: %1.12E\n", i, atomwcaEnergy[i]);
4709 Vnm_print(1,
"\nTotal WCA energy: %g kJ/mol\n", (apolparm->
wcaEnergy));
4710 Vnm_print(1,
"\nTotal non-polar energy = %1.12E kJ/mol\n", energy);
4714 Vnm_print(2,
"energyAPOL: Error in energyAPOL. Unknown option.\n");
4729 time_t ts, ts_main, ts_sub;
4747 Vatom *atom = VNULL;
4750 srad = apolparm->
srad;
4751 press = apolparm->
press;
4752 gamma = apolparm->
gamma;
4753 offset = apolparm->
dpos;
4754 bconc = apolparm->
bconc;
4759 Vnm_print(0,
"forceAPOL: Trying atom surf...\n");
4761 if (acc->
surf == VNULL) {
4763 for (i=0; i<natom; i++) {
4771 Vnm_print(0,
"forceAPOL: atom surf: Time elapsed: %f\n", ((
double)clock() - ts) / CLOCKS_PER_SEC);
4774 Vnm_print(0,
"forceAPOL: calcforce == ACF_TOTAL\n");
4778 if(*atomForce != VNULL){
4779 Vmem_free(mem,*nforce,
sizeof(
AtomForce), (
void **)atomForce);
4782 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4786 for (j=0; j<3; j++) {
4788 (*atomForce)[0].savForce[j] = 0.0;
4789 (*atomForce)[0].wcaForce[j] = 0.0;
4793 for (i=0; i<natom; i++) {
4802 if(VABS(gamma) > VSMALL) {
4805 if(VABS(press) > VSMALL) {
4808 if(VABS(bconc) > VSMALL) {
4813 (*atomForce)[0].sasaForce[j] += dSASA[j];
4814 (*atomForce)[0].savForce[j] += dSAV[j];
4815 (*atomForce)[0].wcaForce[j] += force[j];
4820 Vnm_tprint( 1,
" Printing net forces (kJ/mol/A)\n");
4821 Vnm_tprint( 1,
" Legend:\n");
4822 Vnm_tprint( 1,
" sasa -- SASA force\n");
4823 Vnm_tprint( 1,
" sav -- SAV force\n");
4824 Vnm_tprint( 1,
" wca -- WCA force\n\n");
4826 Vnm_tprint( 1,
" sasa %4.3e %4.3e %4.3e\n",
4827 (*atomForce)[0].sasaForce[0],
4828 (*atomForce)[0].sasaForce[1],
4829 (*atomForce)[0].sasaForce[2]);
4830 Vnm_tprint( 1,
" sav %4.3e %4.3e %4.3e\n",
4831 (*atomForce)[0].savForce[0],
4832 (*atomForce)[0].savForce[1],
4833 (*atomForce)[0].savForce[2]);
4834 Vnm_tprint( 1,
" wca %4.3e %4.3e %4.3e\n",
4835 (*atomForce)[0].wcaForce[0],
4836 (*atomForce)[0].wcaForce[1],
4837 (*atomForce)[0].wcaForce[2]);
4839 Vnm_print(0,
"forceAPOL: calcforce == ACF_TOTAL: %f\n", ((
double)clock() - ts) / CLOCKS_PER_SEC);
4842 if(*atomForce == VNULL){
4843 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4846 Vmem_free(mem,*nforce,
sizeof(
AtomForce), (
void **)atomForce);
4847 *atomForce = (
AtomForce *)Vmem_malloc(mem, *nforce,
4852 Vnm_tprint( 1,
" Printing per atom forces (kJ/mol/A)\n");
4853 Vnm_tprint( 1,
" Legend:\n");
4854 Vnm_tprint( 1,
" tot n -- Total force for atom n\n");
4855 Vnm_tprint( 1,
" sasa n -- SASA force for atom n\n");
4856 Vnm_tprint( 1,
" sav n -- SAV force for atom n\n");
4857 Vnm_tprint( 1,
" wca n -- WCA force for atom n\n\n");
4859 Vnm_tprint( 1,
" gamma %f\n" \
4865 for (i=0; i<natom; i++) {
4875 for (j=0; j<3; j++) {
4876 (*atomForce)[i].sasaForce[j] = 0.0;
4877 (*atomForce)[i].savForce[j] = 0.0;
4878 (*atomForce)[i].wcaForce[j] = 0.0;
4881 if(VABS(gamma) > VSMALL)
Vacc_atomdSASA(acc, offset, srad, atom, dSASA);
4882 if(VABS(press) > VSMALL)
Vacc_atomdSAV(acc, srad, atom, dSAV);
4885 xF = -((gamma*dSASA[0]) + (press*dSAV[0]) + (bconc*force[0]));
4886 yF = -((gamma*dSASA[1]) + (press*dSAV[1]) + (bconc*force[1]));
4887 zF = -((gamma*dSASA[2]) + (press*dSAV[2]) + (bconc*force[2]));
4890 (*atomForce)[i].sasaForce[j] += dSASA[j];
4891 (*atomForce)[i].savForce[j] += dSAV[j];
4892 (*atomForce)[i].wcaForce[j] += force[j];
4896 Vnm_print( 1,
" tot %i %4.3e %4.3e %4.3e\n",
4901 Vnm_print( 1,
" sasa %i %4.3e %4.3e %4.3e\n",
4903 (*atomForce)[i].sasaForce[0],
4904 (*atomForce)[i].sasaForce[1],
4905 (*atomForce)[i].sasaForce[2]);
4906 Vnm_print( 1,
" sav %i %4.3e %4.3e %4.3e\n",
4908 (*atomForce)[i].savForce[0],
4909 (*atomForce)[i].savForce[1],
4910 (*atomForce)[i].savForce[2]);
4911 Vnm_print( 1,
" wca %i %4.3e %4.3e %4.3e\n",
4913 (*atomForce)[i].wcaForce[0],
4914 (*atomForce)[i].wcaForce[1],
4915 (*atomForce)[i].wcaForce[2]);
4925 Vnm_print(0,
"forceAPOL: Time elapsed: %f\n", ((
double)clock() - ts_main) / CLOCKS_PER_SEC);
4934 VPUBLIC
int initBEM(
int icalc,
4956 Vnm_tprint(1,
"Destroying boundary element structures.\n");
4963 void apbs2tabipb_(TABIPBparm* tabiparm,
4964 TABIPBvars* tabivars);
4977 if (nosh != VNULL) {
4978 if (nosh->
bogus)
return 1;
4983 TABIPBparm *tabiparm = (TABIPBparm*)calloc(1,
sizeof(TABIPBparm));
4985 sprintf(tabiparm->fpath,
"");
4986 strncpy(tabiparm->fname, nosh->
molpath[0],4);
4987 tabiparm->fname[4] =
'\0';
4988 tabiparm->density = pbeparm->
sdens;
4989 tabiparm->probe_radius = pbeparm->
srad;
4991 tabiparm->epsp = pbeparm->
pdie;
4992 tabiparm->epsw = pbeparm->
sdie;
4993 tabiparm->bulk_strength = 0.0;
4994 for (i=0; i<pbeparm->
nion; i++)
4995 tabiparm->bulk_strength += pbeparm->
ionc[i]
4996 *pbeparm->
ionq[i]*pbeparm->
ionq[i];
4997 tabiparm->temp = pbeparm->
temp;
5000 tabiparm->maxparnode = bemparm->
tree_n0;
5001 tabiparm->theta = bemparm->
mac;
5003 tabiparm->mesh_flag = bemparm->
mesh;
5007 tabiparm->output_datafile = bemparm->
outdata;
5009 TABIPBvars *tabivars = (TABIPBvars*)calloc(1,
sizeof(TABIPBvars));
5010 if ((tabivars->chrpos = (
double *) malloc(3 * tabiparm->number_of_lines *
sizeof(
double))) == NULL) {
5011 printf(
"Error in allocating t_chrpos!\n");
5013 if ((tabivars->atmchr = (
double *) malloc(tabiparm->number_of_lines *
sizeof(
double))) == NULL) {
5014 printf(
"Error in allocating t_atmchr!\n");
5016 if ((tabivars->atmrad = (
double *) malloc(tabiparm->number_of_lines *
sizeof(
double))) == NULL) {
5017 printf(
"Error in allocating t_atmrad!\n");
5022 for (i = 0; i < tabiparm->number_of_lines; i++){
5032 apbs2tabipb_(tabiparm, tabivars);
5035 free(tabivars->chrpos);
5036 free(tabivars->atmchr);
5037 free(tabivars->atmrad);
5038 free(tabivars->vert_ptl);
5039 free(tabivars->xvct);
5040 free_matrix(tabivars->vert);
5041 free_matrix(tabivars->snrm);
5042 free_matrix(tabivars->face);
5044 Vnm_tprint(1,
"\n\nReturning to APBS caller...\n\n");
5045 Vnm_tprint(1,
"Solvation energy and Coulombic energy in kCal/mol...\n\n");
5046 Vnm_tprint(1,
" Global net ELEC energy = %1.12E\n", tabivars->soleng);
5047 Vnm_tprint(1,
" Global net COULOMBIC energy = %1.12E\n\n", tabivars->couleng);
5057 VPUBLIC
int setPartBEM(
NOsh *nosh,
5065 if (nosh->
bogus)
return 1;
5071 VPUBLIC
int energyBEM(
NOsh *nosh,
5097 VPUBLIC
int forceBEM(
5115 Vnm_tprint( 1,
" Calculating forces...\n");
5123 VPUBLIC
void printBEMPARM(
BEMparm *bemparm) {
5128 VPUBLIC
int writedataBEM(
int rank,
5137 VPUBLIC
int writematBEM(
int rank,
NOsh *nosh,
PBEparm *pbeparm) {
5140 if (nosh->
bogus)
return 1;
5147 #ifdef ENABLE_GEOFLOW 5159 if (nosh != VNULL) {
5160 if (nosh->
bogus)
return 1;
5165 struct GeometricFlowInput geoflowIn = getGeometricFlowParams();
5168 geoflowIn.m_boundaryCondition = (int) pbeparm->
bcfl ;
5169 geoflowIn.m_grid[0] = apolparm->
grid[0];
5170 geoflowIn.m_grid[1] = apolparm->
grid[1];
5171 geoflowIn.m_grid[2] = apolparm->
grid[2];
5172 geoflowIn.m_gamma = apolparm->
gamma;
5173 geoflowIn.m_pdie = pbeparm->
pdie ;
5174 geoflowIn.m_sdie = pbeparm->
sdie ;
5175 geoflowIn.m_press = apolparm->
press ;
5176 geoflowIn.m_tol = parm->
etol;
5177 geoflowIn.m_bconc = apolparm->
bconc ;
5178 geoflowIn.m_vdwdispersion = parm->vdw;
5179 geoflowIn.m_etolSolvation = .01 ;
5185 struct GeometricFlowOutput geoflowOut =
5186 runGeometricFlowWrapAPBS( geoflowIn, molecules[0] );
5188 Vnm_tprint( 1,
" Global net energy = %1.12E\n", geoflowOut.m_totalSolvation);
5189 Vnm_tprint( 1,
" Global net ELEC energy = %1.12E\n", geoflowOut.m_elecSolvation);
5190 Vnm_tprint( 1,
" Global net APOL energy = %1.12E\n", geoflowOut.m_nonpolarSolvation);
5210 printf(
"solvePBAM!!!\n");
5211 if (nosh != VNULL) {
5212 if (nosh->
bogus)
return 1;
5217 PBAMInput pbamIn = getPBAMParams();
5219 pbamIn.nmol_ = nosh->
nmol;
5222 pbamIn.temp_ = pbeparm->
temp;
5223 if (fabs(pbamIn.temp_-0.0) < 1e-3)
5225 printf(
"No temperature specified. Setting to 298.15K\n");
5226 pbamIn.temp_ = 298.15;
5230 pbamIn.idiel_ = pbeparm->
pdie;
5231 pbamIn.sdiel_ = pbeparm->
sdie;
5234 pbamIn.salt_ = parm->salt;
5237 strncpy(pbamIn.runType_, parm->runtype,
CHR_MAXLEN);
5238 strncpy(pbamIn.runName_, parm->runname,
CHR_MAXLEN);
5240 pbamIn.randOrient_ = parm->setrandorient;
5242 pbamIn.boxLen_ = parm->pbcboxlen;
5243 pbamIn.pbcType_ = parm->setpbcs;
5245 pbamIn.setunits_ = parm->setunits;
5246 if(parm->setunits == 1) strncpy(pbamIn.units_, parm->units,
CHR_MAXLEN);
5249 if (parm->setgridpt) pbamIn.gridPts_ = parm->gridpt;
5250 strncpy(pbamIn.map3D_, parm->map3dname,
CHR_MAXLEN);
5251 pbamIn.grid2Dct_ = parm->grid2Dct;
5252 for (i=0; i<pbamIn.grid2Dct_; i++)
5254 strncpy(pbamIn.grid2D_[i], parm->grid2Dname[i],
CHR_MAXLEN);
5255 strncpy(pbamIn.grid2Dax_[i], parm->grid2Dax[i],
CHR_MAXLEN);
5256 pbamIn.grid2Dloc_[i] = parm->grid2Dloc[i];
5258 strncpy(pbamIn.dxname_, parm->dxname,
CHR_MAXLEN);
5261 pbamIn.ntraj_ = parm->ntraj;
5262 strncpy(pbamIn.termCombine_, parm->termcombine,
CHR_MAXLEN);
5264 pbamIn.termct_ = parm->termct;
5265 pbamIn.contct_ = parm->confilct;
5267 if (strncmp(pbamIn.runType_,
"dynamics", 8)== 0)
5269 if (pbamIn.nmol_ > parm->diffct)
5271 Vnm_tprint(2,
"You need more diffusion information!\n");
5275 for (i=0; i<pbamIn.nmol_; i++)
5277 if (parm->xyzct[i] < parm->ntraj)
5279 Vnm_tprint(2,
"For molecule %d, you are missing trajectory!\n", i+1);
5282 for (j=0; j<pbamIn.ntraj_; j++)
5284 strncpy(pbamIn.xyzfil_[i][j], parm->xyzfil[i][j],
CHR_MAXLEN);
5289 for (i=0; i<pbamIn.nmol_; i++)
5291 strncpy(pbamIn.moveType_[i], parm->moveType[i],
CHR_MAXLEN);
5292 pbamIn.transDiff_[i] = parm->transDiff[i];
5293 pbamIn.rotDiff_[i] = parm->rotDiff[i];
5296 for (i=0; i<pbamIn.termct_; i++)
5298 strncpy(pbamIn.termnam_[i], parm->termnam[i],
CHR_MAXLEN);
5299 pbamIn.termnu_[i][0] = parm->termnu[i][0];
5300 pbamIn.termval_[i] = parm->termVal[i];
5303 for (i=0; i<pbamIn.contct_; i++)
5305 strncpy(pbamIn.confil_[i], parm->confil[i],
CHR_MAXLEN);
5311 printPBAMStruct( pbamIn );
5314 PBAMOutput pbamOut = runPBAMWrapAPBS( pbamIn, molecules, nosh->
nmol );
5316 Vnm_tprint(1,
"\n\nReturning to APBS caller...\n\n");
5318 if (!(strncmp(pbamIn.runType_,
"dynamics", 8) &&
5319 strncmp(pbamIn.runType_,
"energyforce", 11))) {
5321 if (!strncmp(pbamIn.units_,
"kcalmol", 7)) {
5323 Vnm_tprint(1,
"Interaction energy in kCal/mol...\n\n");
5325 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5326 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5327 i+1, pbamOut.energies_[i]);
5328 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5329 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5330 pbamOut.forces_[i][2]);
5331 if (pbamOut.energies_[i+1] == 0.)
break;
5334 }
else if (!strncmp(pbamIn.units_,
"jmol", 4)) {
5336 Vnm_tprint(1,
"Interaction energy in J/mol...\n\n");
5338 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5339 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5340 i+1, pbamOut.energies_[i]);
5341 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5342 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5343 pbamOut.forces_[i][2]);
5344 if (pbamOut.energies_[i+1] == 0.)
break;
5350 Vnm_tprint(1,
"Interaction energy in kT @ %6.2f K...\n\n", pbamIn.temp_);
5352 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5353 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5354 i+1, pbamOut.energies_[i]);
5355 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5356 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5357 pbamOut.forces_[i][2]);
5358 if (pbamOut.energies_[i+1] == 0.)
break;
5381 printf(
"solvePBSAM!!!\n");
5382 char fname_tp[VMAX_ARGLEN];
5383 if (nosh != VNULL) {
5384 if (nosh->
bogus)
return 1;
5389 PBSAMInput pbsamIn = getPBSAMParams();
5392 pbamIn.nmol_ = nosh->
nmol;
5395 pbamIn.temp_ = pbeparm->
temp;
5396 if (fabs(pbamIn.temp_-0.0) < 1e-3)
5398 printf(
"No temperature specified. Setting to 298.15K\n");
5399 pbamIn.temp_ = 298.15;
5403 pbamIn.idiel_ = pbeparm->
pdie;
5404 pbamIn.sdiel_ = pbeparm->
sdie;
5407 pbamIn.salt_ = parm->salt;
5410 strncpy(pbamIn.runType_, parm->runtype,
CHR_MAXLEN);
5411 strncpy(pbamIn.runName_, parm->runname,
CHR_MAXLEN);
5413 pbamIn.setunits_ = parm->setunits;
5414 if(parm->setunits == 1) strncpy(pbamIn.units_, parm->units,
CHR_MAXLEN);
5415 pbamIn.randOrient_ = parm->setrandorient;
5417 pbamIn.boxLen_ = parm->pbcboxlen;
5418 pbamIn.pbcType_ = parm->setpbcs;
5421 if (parm->setgridpt) pbamIn.gridPts_ = parm->gridpt;
5422 strncpy(pbamIn.map3D_, parm->map3dname,
CHR_MAXLEN);
5423 pbamIn.grid2Dct_ = parm->grid2Dct;
5424 for (i=0; i<pbamIn.grid2Dct_; i++)
5426 strncpy(pbamIn.grid2D_[i], parm->grid2Dname[i],
CHR_MAXLEN);
5427 strncpy(pbamIn.grid2Dax_[i], parm->grid2Dax[i],
CHR_MAXLEN);
5428 pbamIn.grid2Dloc_[i] = parm->grid2Dloc[i];
5430 strncpy(pbamIn.dxname_, parm->dxname,
CHR_MAXLEN);
5433 pbamIn.ntraj_ = parm->ntraj;
5434 strncpy(pbamIn.termCombine_, parm->termcombine,
CHR_MAXLEN);
5436 pbamIn.termct_ = parm->termct;
5437 pbamIn.contct_ = parm->confilct;
5439 if (strncmp(pbamIn.runType_,
"dynamics", 8)== 0)
5441 if (pbamIn.nmol_ > parm->diffct)
5443 Vnm_tprint(2,
"You need more diffusion information!\n");
5447 for (i=0; i<pbamIn.nmol_; i++)
5449 if (parm->xyzct[i] < parm->ntraj)
5451 Vnm_tprint(2,
"For molecule %d, you are missing trajectory!\n", i+1);
5454 for (j=0; j<pbamIn.ntraj_; j++)
5456 strncpy(pbamIn.xyzfil_[i][j], parm->xyzfil[i][j],
CHR_MAXLEN);
5461 for (i=0; i<pbamIn.nmol_; i++)
5463 strncpy(pbamIn.moveType_[i], parm->moveType[i],
CHR_MAXLEN);
5464 pbamIn.transDiff_[i] = parm->transDiff[i];
5465 pbamIn.rotDiff_[i] = parm->rotDiff[i];
5468 for (i=0; i<pbamIn.termct_; i++)
5470 strncpy(pbamIn.termnam_[i], parm->termnam[i],
CHR_MAXLEN);
5471 pbamIn.termnu_[i][0] = parm->termnu[i][0];
5472 pbamIn.termval_[i] = parm->termVal[i];
5475 for (i=0; i<pbamIn.contct_; i++)
5477 strncpy(pbamIn.confil_[i], parm->confil[i],
CHR_MAXLEN);
5483 pbsamIn.tolsp_ = samparm->tolsp;
5484 pbsamIn.imatct_ = samparm->imatct;
5485 pbsamIn.expct_ = samparm->expct;
5486 for (i=0; i<samparm->surfct; i++)
5488 strncpy(pbsamIn.surffil_[i], samparm->surffil[i],
CHR_MAXLEN);
5490 for (i=0; i<samparm->imatct; i++)
5492 strncpy(pbsamIn.imatfil_[i], samparm->imatfil[i],
CHR_MAXLEN);
5494 for (i=0; i<samparm->expct; i++)
5496 strncpy(pbsamIn.expfil_[i], samparm->expfil[i],
CHR_MAXLEN);
5500 if (samparm->setmsms == 1) {
5501 for (i=0; i<pbamIn.nmol_; i++) {
5503 for (j=0; j < VMAX_ARGLEN; j++)
5504 if (nosh->
molpath[i][j] ==
'\0')
break;
5507 char xyzr[j+2], surf[j+2], outname[j-4];
5508 for( k=0; k < j - 4; k++)
5510 xyzr[k] = nosh->
molpath[i][k];
5511 outname[k] = nosh->
molpath[i][k];
5512 surf[k] = nosh->
molpath[i][k];
5515 xyzr[k] =
'.'; surf[k] =
'.';
5516 xyzr[k+1] =
'x'; surf[k+1] =
'v';
5517 xyzr[k+2] =
'y'; surf[k+2] =
'e';
5518 xyzr[k+3] =
'z'; surf[k+3] =
'r';
5519 xyzr[k+4] =
'r'; surf[k+4] =
't';
5520 xyzr[k+5] =
'\0'; surf[k+5] =
'\0';;
5524 fp=fopen(xyzr,
"w");
5537 sprintf(fname_tp,
"msms.exe -if %s -prob %f -dens %f -of %s",
5538 xyzr, samparm->probe_radius,samparm->density, outname);
5540 sprintf(fname_tp,
"msms -if %s -prob %f -dens %f -of %s",
5541 xyzr, samparm->probe_radius,samparm->density, outname);
5544 printf(
"%s\n", fname_tp);
5546 printf(
"Running MSMS...\n");
5547 ierr = system(fname_tp);
5549 strncpy(pbsamIn.surffil_[i], surf,
CHR_MAXLEN);
5555 printPBSAMStruct( pbamIn, pbsamIn );
5558 PBAMOutput pbamOut = runPBSAMWrapAPBS(pbamIn, pbsamIn, molecules, nosh->
nmol);
5560 Vnm_tprint(1,
"\n\nReturning to APBS caller...\n\n");
5562 if (!(strncmp(pbamIn.runType_,
"dynamics", 8) &&
5563 strncmp(pbamIn.runType_,
"energyforce", 11))) {
5565 if (!strncmp(pbamIn.units_,
"kcalmol", 7)) {
5567 Vnm_tprint(1,
"Interaction energy in kCal/mol...\n\n");
5569 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5570 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5571 i+1, pbamOut.energies_[i]);
5572 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5573 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5574 pbamOut.forces_[i][2]);
5575 if (pbamOut.energies_[i+1] == 0.)
break;
5578 }
else if (!strncmp(pbamIn.units_,
"jmol", 4)) {
5580 Vnm_tprint(1,
"Interaction energy in J/mol...\n\n");
5582 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5583 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5584 i+1, pbamOut.energies_[i]);
5585 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5586 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5587 pbamOut.forces_[i][2]);
5588 if (pbamOut.energies_[i+1] == 0.)
break;
5594 Vnm_tprint(1,
"Interaction energy in kT @ %6.2f K...\n\n", pbamIn.temp_);
5596 for (
int i = 0; i < PBAMPARM_MAXMOL; i++) {
5597 Vnm_tprint(1,
" Molecule %d: Global net ELEC energy = %1.12E\n",
5598 i+1, pbamOut.energies_[i]);
5599 Vnm_tprint(1,
" Force = (%1.6E, %1.6E, %1.6E)\n\n",
5600 pbamOut.forces_[i][0], pbamOut.forces_[i][1],
5601 pbamOut.forces_[i][2]);
5602 if (pbamOut.energies_[i+1] == 0.)
break;
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 Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
VPUBLIC double Vpmg_qmEnergy(Vpmg *thee, int extFlag)
Get the "mobile charge" contribution to the electrostatic energy.
VPUBLIC double Vpmg_qfEnergy(Vpmg *thee, int extFlag)
Get the "fixed charge" contribution to the electrostatic energy.
Vdata_Format meshfmt[NOSH_MAXMOL]
VPUBLIC void Vgrid_writeDX(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out the data in OpenDX grid format.
Parameter structure for GEOFLOW-specific variables from input files.
int apol2calc[NOSH_MAXCALC]
int Vacc_wcaEnergyAtom(Vacc *thee, APOLparm *apolparm, Valist *alist, Vclist *clist, int iatom, double *value)
Calculate the WCA energy for an atom.
VPUBLIC Vparam * loadParameter(NOsh *nosh)
Loads and returns parameter object.
VPUBLIC void Vfetk_setAtomColors(Vfetk *thee)
Transfer color (partition ID) information frmo a partitioned mesh to the atoms.
VPUBLIC void Vgrid_writeUHBD(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out the data in UHBD grid format.
VPUBLIC Vfetk * Vfetk_ctor(Vpbe *pbe, Vhal_PBEType type)
Constructor for Vfetk object.
#define CHR_MAXLEN
Number of things that can be written out in a single calculation.
VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom)
Return the atomic solvent accessible surface area (SASA)
#define NOSH_MAXMOL
Maximum number of molecules in a run.
Contains public data members for Vpbe class/module.
Oracle for solvent- and ion-accessibility around a biomolecule.
Vdata_Format dielfmt[NOSH_MAXMOL]
VPUBLIC Vparam * Vparam_ctor()
Construct the object.
VPUBLIC int loadChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the charge maps given in NOsh into grid objects.
VPUBLIC void Vpmgp_dtor(Vpmgp **thee)
Object destructor.
VPUBLIC void Vgrid_dtor(Vgrid **thee)
Object destructor.
enum eBEMparm_CalcType BEMparm_CalcType
Declare BEMparm_CalcType type.
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
char potpath[NOSH_MAXMOL][VMAX_ARGLEN]
Parameter structure for PBAM-specific variables from input files.
VPUBLIC int energyAPOL(APOLparm *apolparm, double sasa, double sav, double atomsasa[], double atomwcaEnergy[], int numatoms)
Calculate non-polar energies.
NOsh_PrintType printwhat[NOSH_MAXPRINT]
Electrostatic potential oracle for Cartesian mesh data.
VPUBLIC int Vgrid_readGZ(Vgrid *thee, const char *fname)
Read in OpenDX data in GZIP format.
VPUBLIC int writedataFlat(NOsh *nosh, Vcom *com, const char *fname, double totEnergy[NOSH_MAXCALC], double qfEnergy[NOSH_MAXCALC], double qmEnergy[NOSH_MAXCALC], double dielEnergy[NOSH_MAXCALC], int nenergy[NOSH_MAXCALC], double *atomEnergy[NOSH_MAXCALC], int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Write out information to a flat file.
Vdata_Format writefmt[PBEPARM_MAXWRITE]
char elecname[NOSH_MAXCALC][VMAX_ARGLEN]
VPUBLIC int postRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Estimate error, mark mesh, and refine mesh after solve.
NOsh_calc * calc[NOSH_MAXCALC]
Contains public data members for Vpmg class/module.
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
VPUBLIC int solveMG(NOsh *nosh, Vpmg *pmg, MGparm_CalcType type)
Solve the PBE with MG.
VPUBLIC void Vclist_dtor(Vclist **thee)
Destroy object.
VPUBLIC void Vpmg_setPart(Vpmg *thee, double lowerCorner[3], double upperCorner[3], int bflags[6])
Set partition information which restricts the calculation of observables to a (rectangular) subset of...
VPUBLIC Vpbe * Vpbe_ctor(Valist *alist, int ionNum, double *ionConc, double *ionRadii, double *ionQ, double T, double soluteDiel, double solventDiel, double solventRadius, int focusFlag, double sdens, double z_mem, double L, double membraneDiel, double V)
Construct Vpbe object.
VPUBLIC int writematMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out operator matrix from MG calculation to file.
VPUBLIC int setPartMG(NOsh *nosh, MGparm *mgparm, Vpmg *pmg)
Set MG partitions for calculating observables and performing I/O.
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.
Header file for front end auxiliary routines.
VPUBLIC int Vpmg_dbForce(Vpmg *thee, double *dbForce, int atomID, Vsurf_Meth srfm)
Calculate the dielectric boundary forces on the specified atom in units of k_B T/AA.
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
VPUBLIC void killChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded charge maps.
#define APBS_TIMER_SETUP
APBS setup timer ID.
VPUBLIC Vclist * Vclist_ctor(Valist *alist, double max_radius, int npts[VAPBS_DIM], Vclist_DomainMode mode, double lower_corner[VAPBS_DIM], double upper_corner[VAPBS_DIM])
Construct the cell list object.
VPUBLIC void killForce(Vmem *mem, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Free memory from MG force calculation.
APOLparm_calcForce calcforce
char dielXpath[NOSH_MAXMOL][VMAX_ARGLEN]
VPUBLIC int loadKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the kappa maps given in NOsh into grid objects.
VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius)
Return the total solvent accessible surface area (SASA)
VPUBLIC int loadDielMaps(NOsh *nosh, Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL])
Load the dielectric maps given in NOsh into grid objects.
VPUBLIC int Vacc_wcaForceAtom(Vacc *thee, APOLparm *apolparm, Vclist *clist, Vatom *atom, double *force)
Return the WCA integral force.
VPUBLIC Vrc_Codes Vfetk_loadMesh(Vfetk *thee, double center[3], double length[3], Vfetk_MeshLoad meshType, Vio *sock)
Loads a mesh into the Vfetk (and associated) object(s).
VPUBLIC int initMG(int icalc, NOsh *nosh, MGparm *mgparm, PBEparm *pbeparm, double realCenter[3], Vpbe *pbe[NOSH_MAXCALC], Valist *alist[NOSH_MAXMOL], Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL], Vgrid *kappaMap[NOSH_MAXMOL], Vgrid *chargeMap[NOSH_MAXMOL], Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC], Vgrid *potMap[NOSH_MAXMOL])
Initialize an MG calculation.
#define APBS_TIMER_SOLVER
APBS solver timer ID.
VPUBLIC void storeAtomEnergy(Vpmg *pmg, int icalc, double **atomEnergy, int *nenergy)
Store energy in arrays for future use.
VPUBLIC void killMG(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC])
Kill structures initialized during an MG calculation.
AtomData sub-class; stores atom data.
VPUBLIC int printApolForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data.
VPUBLIC void printFEPARM(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Print out FE-specific params loaded from input.
VPUBLIC void Vfetk_dtor(Vfetk **thee)
Object destructor.
VPUBLIC double Vpmg_qfAtomEnergy(Vpmg *thee, Vatom *atom)
Get the per-atom "fixed charge" contribution to the electrostatic energy.
VPUBLIC void killPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded potential maps.
VPUBLIC void Vfetk_setParameters(Vfetk *thee, PBEparm *pbeparm, FEMparm *feparm)
Set the parameter objects.
VPUBLIC void killFE(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vfetk *fetk[NOSH_MAXCALC], Gem *gm[NOSH_MAXMOL])
Kill structures initialized during an FE calculation.
Vdata_Format potfmt[NOSH_MAXMOL]
VPUBLIC Vrc_Codes Valist_readXML(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from an XML file.
VPUBLIC void Vpmg_dtor(Vpmg **thee)
Object destructor.
VPUBLIC int printElecEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data.
Vdata_Format kappafmt[NOSH_MAXMOL]
VPUBLIC int partFE(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Partition mesh (if applicable)
Parameter structure for PBSAM-specific variables from input files.
VPUBLIC void killKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded kappa maps.
APOLparm_calcEnergy calcenergy
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius)
Return the total solvent accessible volume (SAV)
enum eMGparm_CalcType MGparm_CalcType
Declare MGparm_CalcType type.
VPUBLIC void Vacc_atomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible area.
VPUBLIC int Vparam_readFlatFile(Vparam *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read a flat-file format parameter database.
VPUBLIC int Vpmg_fillArray(Vpmg *thee, double *vec, Vdata_Type type, double parm, Vhal_PBEType pbetype, PBEparm *pbeparm)
Fill the specified array with accessibility values.
VPUBLIC int Vstring_strcasecmp(const char *s1, const char *s2)
Case-insensitive string comparison (BSD standard)
VPUBLIC int energyFE(NOsh *nosh, int icalc, Vfetk *fetk[NOSH_MAXCALC], int *nenergy, double *totEnergy, double *qfEnergy, double *qmEnergy, double *dielEnergy)
Calculate electrostatic energies from FE solution.
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
VPUBLIC int printForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data (deprecated...see printElecForce)
VPUBLIC void printPBEPARM(PBEparm *pbeparm)
Print out generic PBE params loaded from input.
#define Vunit_Na
Avogadro's number.
VPUBLIC int writedataMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out observables from MG calculation to file.
Parameter structure for BEM-specific variables from input files.
VPUBLIC Vparam_AtomData * Vparam_getAtomData(Vparam *thee, char resName[VMAX_ARGLEN], char atomName[VMAX_ARGLEN])
Get atom data.
Parameter structure for FEM-specific variables from input files.
Structure to hold atomic forces.
Reads and assigns charge/radii parameters.
int printcalc[NOSH_MAXPRINT][NOSH_MAXPOP]
VPUBLIC double Vpmg_dielEnergy(Vpmg *thee, int extFlag)
Get the "polarization" contribution to the electrostatic energy.
VPUBLIC Vpmg * Vpmg_ctor(Vpmgp *pmgp, Vpbe *pbe, int focusFlag, Vpmg *pmgOLD, MGparm *mgparm, PBEparm_calcEnergy energyFlag)
Constructor for the Vpmg class (allocates new memory)
Contains public data members for Vpmgp class/module.
#define Vunit_kb
Boltzmann constant.
PBEparm_calcEnergy calcenergy
VPUBLIC int printEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data (deprecated...see printElecEnergy)
VPUBLIC double returnEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Access net local energy.
#define APBS_TIMER_ENERGY
APBS energy timer ID.
VPUBLIC int loadPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the potential maps given in NOsh into grid objects.
NOsh_MolFormat molfmt[NOSH_MAXMOL]
char chargepath[NOSH_MAXMOL][VMAX_ARGLEN]
Parameter structure for PBE variables from input files.
VPUBLIC int Vpmg_fillco(Vpmg *thee, Vsurf_Meth surfMeth, double splineWin, Vchrg_Meth chargeMeth, int useDielXMap, Vgrid *dielXMap, int useDielYMap, Vgrid *dielYMap, int useDielZMap, Vgrid *dielZMap, int useKappaMap, Vgrid *kappaMap, int usePotMap, Vgrid *potMap, int useChargeMap, Vgrid *chargeMap)
Fill the coefficient arrays prior to solving the equation.
VPUBLIC void killMolecules(NOsh *nosh, Valist *alist[NOSH_MAXMOL])
Destroy the loaded molecules.
PBEparm_calcForce calcforce
int printnarg[NOSH_MAXPRINT]
VPUBLIC int forceMG(Vmem *mem, NOsh *nosh, PBEparm *pbeparm, MGparm *mgparm, Vpmg *pmg, int *nforce, AtomForce **atomForce, Valist *alist[NOSH_MAXMOL])
Calculate forces from MG solution.
Surface object list of per-atom surface points.
VPUBLIC Vrc_Codes Valist_readPDB(Valist *thee, Vparam *param, Vio *sock)
Fill atom list with information from a PDB file.
char dielZpath[NOSH_MAXMOL][VMAX_ARGLEN]
FEMparm_EstType akeySOLVE
Parameter structure for MG-specific variables from input files.
Vdata_Format chargefmt[NOSH_MAXMOL]
VPUBLIC void Vpbe_dtor(Vpbe **thee)
Object destructor.
VPUBLIC void killDielMaps(NOsh *nosh, Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL])
Destroy the loaded dielectric.
VPUBLIC Vrc_Codes initFE(int icalc, NOsh *nosh, FEMparm *feparm, PBEparm *pbeparm, Vpbe *pbe[NOSH_MAXCALC], Valist *alist[NOSH_MAXMOL], Vfetk *fetk[NOSH_MAXCALC])
Initialize FE solver objects.
VPUBLIC int printElecForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data.
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
VPUBLIC void printMGPARM(MGparm *mgparm, double realCenter[3])
Print out MG-specific params loaded from input.
VPUBLIC int printApolEnergy(NOsh *nosh, int iprint)
Combine and pretty-print energy data.
Contains public data members for Vatom class/module.
char writematstem[VMAX_ARGLEN]
VPUBLIC int Vgrid_readDX(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read in data in OpenDX grid format.
Container class for list of atom objects.
VPUBLIC Vpmgp * Vpmgp_ctor(MGparm *mgparm)
Construct PMG parameter object and initialize to default values.
Class for parsing fixed format input files.
char meshpath[NOSH_MAXMOL][VMAX_ARGLEN]
VPUBLIC int Vparam_readXMLFile(Vparam *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read an XML format parameter database.
int elec2calc[NOSH_MAXCALC]
Contains public data members for Vfetk class/module.
VPUBLIC int Vpmg_ibForce(Vpmg *thee, double *force, int atomID, Vsurf_Meth srfm)
Calculate the osmotic pressure on the specified atom in units of k_B T/AA.
VPUBLIC int loadMolecules(NOsh *nosh, Vparam *param, Valist *alist[NOSH_MAXMOL])
Load the molecules given in NOsh into atom lists.
VPUBLIC void Vgrid_writeGZ(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out OpenDX data in GZIP format.
#define NOSH_MAXCALC
Maximum number of calculations in a run.
VPUBLIC int forceAPOL(Vacc *acc, Vmem *mem, APOLparm *apolparm, int *nforce, AtomForce **atomForce, Valist *alist, Vclist *clist)
Calculate non-polar forces.
Vdata_Type writetype[PBEPARM_MAXWRITE]
VPUBLIC int Vpmg_qfForce(Vpmg *thee, double *force, int atomID, Vchrg_Meth chgm)
Calculate the "charge-field" force on the specified atom in units of k_B T/AA.
VPUBLIC int solveFE(int icalc, PBEparm *pbeparm, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Solve-estimate-refine.
VPUBLIC Vrc_Codes Valist_readPQR(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from a PQR file.
VPUBLIC double Vfetk_energy(Vfetk *thee, int color, int nonlin)
Return the total electrostatic energy.
enum eVfetk_MeshLoad Vfetk_MeshLoad
Declare FEMparm_GuessType type.
VPUBLIC void startVio()
Wrapper to start MALOC Vio layer.
char parmpath[VMAX_ARGLEN]
VPUBLIC int Vpmg_solve(Vpmg *thee)
Solve the PBE using PMG.
VPUBLIC double Vpmg_energy(Vpmg *thee, int extFlag)
Get the total electrostatic energy.
VPUBLIC double Vpbe_getDeblen(Vpbe *thee)
Get Debye-Huckel screening length.
int printop[NOSH_MAXPRINT][NOSH_MAXPOP]
char molpath[NOSH_MAXMOL][VMAX_ARGLEN]
VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type)
Fill an array with the specified data.
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
char writestem[PBEPARM_MAXWRITE][VMAX_ARGLEN]
VPUBLIC void Valist_dtor(Valist **thee)
Destroys atom list object.
char apolname[NOSH_MAXCALC][VMAX_ARGLEN]
VPUBLIC void Vacc_atomdSAV(Vacc *thee, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible volume.
VPUBLIC int initAPOL(NOsh *nosh, Vmem *mem, Vparam *param, APOLparm *apolparm, int *nforce, AtomForce **atomForce, Valist *alist)
Upperlevel routine to the non-polar energy and force routines.
VPUBLIC int Vfetk_write(Vfetk *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, Bvec *vec, Vdata_Format format)
Write out data.
VPUBLIC Valist * Valist_ctor()
Construct the atom list object.
#define APBS_TIMER_FORCE
APBS force timer ID.
Parameter structure for APOL-specific variables from input files.
VPUBLIC int writedataXML(NOsh *nosh, Vcom *com, const char *fname, double totEnergy[NOSH_MAXCALC], double qfEnergy[NOSH_MAXCALC], double qmEnergy[NOSH_MAXCALC], double dielEnergy[NOSH_MAXCALC], int nenergy[NOSH_MAXCALC], double *atomEnergy[NOSH_MAXCALC], int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Write out information to an XML file.
VPUBLIC VaccSurf * Vacc_atomSurf(Vacc *thee, Vatom *atom, VaccSurf *ref, double prad)
Set up an array of points corresponding to the SAS due to a particular atom.
VPUBLIC int energyMG(NOsh *nosh, int icalc, Vpmg *pmg, int *nenergy, double *totEnergy, double *qfEnergy, double *qmEnergy, double *dielEnergy)
Calculate electrostatic energies from MG solution.
VPUBLIC int preRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Pre-refine mesh before solve.
VPUBLIC int writedataFE(int rank, NOsh *nosh, PBEparm *pbeparm, Vfetk *fetk)
Write FEM data to files.
VPUBLIC Vgrid * Vgrid_ctor(int nx, int ny, int nz, double hx, double hy, double hzed, double xmin, double ymin, double zmin, double *data)
Construct Vgrid object with values obtained from Vpmg_readDX (for example)
VPUBLIC void killEnergy()
Kill arrays allocated for energies.
char dielYpath[NOSH_MAXMOL][VMAX_ARGLEN]
char kappapath[NOSH_MAXMOL][VMAX_ARGLEN]