Fixes issue with vsiten and Verlet buffers
[gromacs.git] / src / tools / gmx_sas.c
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2007, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
11 *
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
16 *
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
34 *
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
37 */
38 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif
41 #include <math.h>
42 #include <stdlib.h>
43
44 #include "sysstuff.h"
45 #include <string.h>
46 #include "typedefs.h"
47 #include "smalloc.h"
48 #include "macros.h"
49 #include "vec.h"
50 #include "xvgr.h"
51 #include "pbc.h"
52 #include "copyrite.h"
53 #include "futil.h"
54 #include "statutil.h"
55 #include "index.h"
56 #include "nsc.h"
57 #include "pdbio.h"
58 #include "confio.h"
59 #include "rmpbc.h"
60 #include "names.h"
61 #include "atomprop.h"
62 #include "physics.h"
63 #include "tpxio.h"
64 #include "gmx_ana.h"
65
66
67 typedef struct {
68 atom_id aa, ab;
69 real d2a, d2b;
70 } t_conect;
71
72 void add_rec(t_conect c[], atom_id i, atom_id j, real d2)
73 {
74 if (c[i].aa == NO_ATID)
75 {
76 c[i].aa = j;
77 c[i].d2a = d2;
78 }
79 else if (c[i].ab == NO_ATID)
80 {
81 c[i].ab = j;
82 c[i].d2b = d2;
83 }
84 else if (d2 < c[i].d2a)
85 {
86 c[i].aa = j;
87 c[i].d2a = d2;
88 }
89 else if (d2 < c[i].d2b)
90 {
91 c[i].ab = j;
92 c[i].d2b = d2;
93 }
94 /* Swap them if necessary: a must be larger than b */
95 if (c[i].d2a < c[i].d2b)
96 {
97 j = c[i].ab;
98 c[i].ab = c[i].aa;
99 c[i].aa = j;
100 d2 = c[i].d2b;
101 c[i].d2b = c[i].d2a;
102 c[i].d2a = d2;
103 }
104 }
105
106 void do_conect(const char *fn, int n, rvec x[])
107 {
108 FILE *fp;
109 int i, j;
110 t_conect *c;
111 rvec dx;
112 real d2;
113
114 fprintf(stderr, "Building CONECT records\n");
115 snew(c, n);
116 for (i = 0; (i < n); i++)
117 {
118 c[i].aa = c[i].ab = NO_ATID;
119 }
120
121 for (i = 0; (i < n); i++)
122 {
123 for (j = i+1; (j < n); j++)
124 {
125 rvec_sub(x[i], x[j], dx);
126 d2 = iprod(dx, dx);
127 add_rec(c, i, j, d2);
128 add_rec(c, j, i, d2);
129 }
130 }
131 fp = ffopen(fn, "a");
132 for (i = 0; (i < n); i++)
133 {
134 if ((c[i].aa == NO_ATID) || (c[i].ab == NO_ATID))
135 {
136 fprintf(stderr, "Warning dot %d has no conections\n", i+1);
137 }
138 fprintf(fp, "CONECT%5d%5d%5d\n", i+1, c[i].aa+1, c[i].ab+1);
139 }
140 ffclose(fp);
141 sfree(c);
142 }
143
144 void connelly_plot(const char *fn, int ndots, real dots[], rvec x[], t_atoms *atoms,
145 t_symtab *symtab, int ePBC, matrix box, gmx_bool bSave)
146 {
147 static const char *atomnm = "DOT";
148 static const char *resnm = "DOT";
149 static const char *title = "Connely Dot Surface Generated by g_sas";
150
151 int i, i0, r0, ii0, k;
152 rvec *xnew;
153 t_atoms aaa;
154
155 if (bSave)
156 {
157 i0 = atoms->nr;
158 r0 = atoms->nres;
159 srenew(atoms->atom, atoms->nr+ndots);
160 srenew(atoms->atomname, atoms->nr+ndots);
161 srenew(atoms->resinfo, r0+1);
162 atoms->atom[i0].resind = r0;
163 t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0+1, ' ', 0, ' ');
164 srenew(atoms->pdbinfo, atoms->nr+ndots);
165 snew(xnew, atoms->nr+ndots);
166 for (i = 0; (i < atoms->nr); i++)
167 {
168 copy_rvec(x[i], xnew[i]);
169 }
170 for (i = k = 0; (i < ndots); i++)
171 {
172 ii0 = i0+i;
173 atoms->atomname[ii0] = put_symtab(symtab, atomnm);
174 atoms->pdbinfo[ii0].type = epdbATOM;
175 atoms->pdbinfo[ii0].atomnr = ii0;
176 atoms->atom[ii0].resind = r0;
177 xnew[ii0][XX] = dots[k++];
178 xnew[ii0][YY] = dots[k++];
179 xnew[ii0][ZZ] = dots[k++];
180 atoms->pdbinfo[ii0].bfac = 0.0;
181 atoms->pdbinfo[ii0].occup = 0.0;
182 }
183 atoms->nr = i0+ndots;
184 atoms->nres = r0+1;
185 write_sto_conf(fn, title, atoms, xnew, NULL, ePBC, box);
186 atoms->nres = r0;
187 atoms->nr = i0;
188 }
189 else
190 {
191 init_t_atoms(&aaa, ndots, TRUE);
192 aaa.atom[0].resind = 0;
193 t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
194 snew(xnew, ndots);
195 for (i = k = 0; (i < ndots); i++)
196 {
197 ii0 = i;
198 aaa.atomname[ii0] = put_symtab(symtab, atomnm);
199 aaa.pdbinfo[ii0].type = epdbATOM;
200 aaa.pdbinfo[ii0].atomnr = ii0;
201 aaa.atom[ii0].resind = 0;
202 xnew[ii0][XX] = dots[k++];
203 xnew[ii0][YY] = dots[k++];
204 xnew[ii0][ZZ] = dots[k++];
205 aaa.pdbinfo[ii0].bfac = 0.0;
206 aaa.pdbinfo[ii0].occup = 0.0;
207 }
208 aaa.nr = ndots;
209 write_sto_conf(fn, title, &aaa, xnew, NULL, ePBC, box);
210 do_conect(fn, ndots, xnew);
211 free_t_atoms(&aaa, FALSE);
212 }
213 sfree(xnew);
214 }
215
216 real calc_radius(char *atom)
217 {
218 real r;
219
220 switch (atom[0])
221 {
222 case 'C':
223 r = 0.16;
224 break;
225 case 'O':
226 r = 0.13;
227 break;
228 case 'N':
229 r = 0.14;
230 break;
231 case 'S':
232 r = 0.2;
233 break;
234 case 'H':
235 r = 0.1;
236 break;
237 default:
238 r = 1e-3;
239 }
240 return r;
241 }
242
243 void sas_plot(int nfile, t_filenm fnm[], real solsize, int ndots,
244 real qcut, gmx_bool bSave, real minarea, gmx_bool bPBC,
245 real dgs_default, gmx_bool bFindex, const output_env_t oenv)
246 {
247 FILE *fp, *fp2, *fp3 = NULL, *vp;
248 const char *flegend[] = {
249 "Hydrophobic", "Hydrophilic",
250 "Total", "D Gsolv"
251 };
252 const char *vlegend[] = { "Volume (nm\\S3\\N)", "Density (g/l)" };
253 const char *or_and_oa_legend[] = { "Average (nm\\S2\\N)", "Standard deviation (nm\\S2\\N)" };
254 const char *vfile;
255 real t;
256 gmx_atomprop_t aps = NULL;
257 gmx_rmpbc_t gpbc = NULL;
258 t_trxstatus *status;
259 int ndefault;
260 int i, j, ii, nfr, natoms, flag, nsurfacedots, res;
261 rvec *xtop, *x;
262 matrix topbox, box;
263 t_topology top;
264 char title[STRLEN];
265 int ePBC;
266 gmx_bool bTop;
267 t_atoms *atoms;
268 gmx_bool *bOut, *bPhobic;
269 gmx_bool bConnelly;
270 gmx_bool bResAt, bITP, bDGsol;
271 real *radius, *dgs_factor = NULL, *area = NULL, *surfacedots = NULL;
272 real at_area, *atom_area = NULL, *atom_area2 = NULL;
273 real *res_a = NULL, *res_area = NULL, *res_area2 = NULL;
274 real totarea, totvolume, totmass = 0, density, harea, tarea, fluc2;
275 atom_id **index, *findex;
276 int *nx, nphobic, npcheck, retval;
277 char **grpname, *fgrpname;
278 real dgsolv;
279
280 bITP = opt2bSet("-i", nfile, fnm);
281 bResAt = opt2bSet("-or", nfile, fnm) || opt2bSet("-oa", nfile, fnm) || bITP;
282
283 bTop = read_tps_conf(ftp2fn(efTPS, nfile, fnm), title, &top, &ePBC,
284 &xtop, NULL, topbox, FALSE);
285 atoms = &(top.atoms);
286
287 if (!bTop)
288 {
289 fprintf(stderr, "No tpr file, will not compute Delta G of solvation\n");
290 bDGsol = FALSE;
291 }
292 else
293 {
294 bDGsol = strcmp(*(atoms->atomtype[0]), "?") != 0;
295 if (!bDGsol)
296 {
297 fprintf(stderr, "Warning: your tpr file is too old, will not compute "
298 "Delta G of solvation\n");
299 }
300 else
301 {
302 printf("In case you use free energy of solvation predictions:\n");
303 please_cite(stdout, "Eisenberg86a");
304 }
305 }
306
307 aps = gmx_atomprop_init();
308
309 if ((natoms = read_first_x(oenv, &status, ftp2fn(efTRX, nfile, fnm),
310 &t, &x, box)) == 0)
311 {
312 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
313 }
314
315 if ((ePBC != epbcXYZ) || (TRICLINIC(box)))
316 {
317 fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
318 "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
319 "will certainly crash the analysis.\n\n");
320 }
321 snew(nx, 2);
322 snew(index, 2);
323 snew(grpname, 2);
324 fprintf(stderr, "Select a group for calculation of surface and a group for output:\n");
325 get_index(atoms, ftp2fn_null(efNDX, nfile, fnm), 2, nx, index, grpname);
326
327 if (bFindex)
328 {
329 fprintf(stderr, "Select a group of hydrophobic atoms:\n");
330 get_index(atoms, ftp2fn_null(efNDX, nfile, fnm), 1, &nphobic, &findex, &fgrpname);
331 }
332 snew(bOut, natoms);
333 for (i = 0; i < nx[1]; i++)
334 {
335 bOut[index[1][i]] = TRUE;
336 }
337
338 /* Now compute atomic readii including solvent probe size */
339 snew(radius, natoms);
340 snew(bPhobic, nx[0]);
341 if (bResAt)
342 {
343 snew(atom_area, nx[0]);
344 snew(atom_area2, nx[0]);
345 snew(res_a, atoms->nres);
346 snew(res_area, atoms->nres);
347 snew(res_area2, atoms->nres);
348 }
349 if (bDGsol)
350 {
351 snew(dgs_factor, nx[0]);
352 }
353
354 /* Get a Van der Waals radius for each atom */
355 ndefault = 0;
356 for (i = 0; (i < natoms); i++)
357 {
358 if (!gmx_atomprop_query(aps, epropVDW,
359 *(atoms->resinfo[atoms->atom[i].resind].name),
360 *(atoms->atomname[i]), &radius[i]))
361 {
362 ndefault++;
363 }
364 /* radius[i] = calc_radius(*(top->atoms.atomname[i])); */
365 radius[i] += solsize;
366 }
367 if (ndefault > 0)
368 {
369 fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
370 }
371 /* Determine which atom is counted as hydrophobic */
372 if (bFindex)
373 {
374 npcheck = 0;
375 for (i = 0; (i < nx[0]); i++)
376 {
377 ii = index[0][i];
378 for (j = 0; (j < nphobic); j++)
379 {
380 if (findex[j] == ii)
381 {
382 bPhobic[i] = TRUE;
383 if (bOut[ii])
384 {
385 npcheck++;
386 }
387 }
388 }
389 }
390 if (npcheck != nphobic)
391 {
392 gmx_fatal(FARGS, "Consistency check failed: not all %d atoms in the hydrophobic index\n"
393 "found in the normal index selection (%d atoms)", nphobic, npcheck);
394 }
395 }
396 else
397 {
398 nphobic = 0;
399 }
400
401 for (i = 0; (i < nx[0]); i++)
402 {
403 ii = index[0][i];
404 if (!bFindex)
405 {
406 bPhobic[i] = fabs(atoms->atom[ii].q) <= qcut;
407 if (bPhobic[i] && bOut[ii])
408 {
409 nphobic++;
410 }
411 }
412 if (bDGsol)
413 {
414 if (!gmx_atomprop_query(aps, epropDGsol,
415 *(atoms->resinfo[atoms->atom[ii].resind].name),
416 *(atoms->atomtype[ii]), &(dgs_factor[i])))
417 {
418 dgs_factor[i] = dgs_default;
419 }
420 }
421 if (debug)
422 {
423 fprintf(debug, "Atom %5d %5s-%5s: q= %6.3f, r= %6.3f, dgsol= %6.3f, hydrophobic= %s\n",
424 ii+1, *(atoms->resinfo[atoms->atom[ii].resind].name),
425 *(atoms->atomname[ii]),
426 atoms->atom[ii].q, radius[ii]-solsize, dgs_factor[i],
427 EBOOL(bPhobic[i]));
428 }
429 }
430 fprintf(stderr, "%d out of %d atoms were classified as hydrophobic\n",
431 nphobic, nx[1]);
432
433 fp = xvgropen(opt2fn("-o", nfile, fnm), "Solvent Accessible Surface", "Time (ps)",
434 "Area (nm\\S2\\N)", oenv);
435 xvgr_legend(fp, asize(flegend) - (bDGsol ? 0 : 1), flegend, oenv);
436 vfile = opt2fn_null("-tv", nfile, fnm);
437 if (vfile)
438 {
439 if (!bTop)
440 {
441 gmx_fatal(FARGS, "Need a tpr file for option -tv");
442 }
443 vp = xvgropen(vfile, "Volume and Density", "Time (ps)", "", oenv);
444 xvgr_legend(vp, asize(vlegend), vlegend, oenv);
445 totmass = 0;
446 ndefault = 0;
447 for (i = 0; (i < nx[0]); i++)
448 {
449 real mm;
450 ii = index[0][i];
451 /*
452 if (!query_atomprop(atomprop,epropMass,
453 *(top->atoms.resname[top->atoms.atom[ii].resnr]),
454 *(top->atoms.atomname[ii]),&mm))
455 ndefault++;
456 totmass += mm;
457 */
458 totmass += atoms->atom[ii].m;
459 }
460 if (ndefault)
461 {
462 fprintf(stderr, "WARNING: Using %d default masses for density calculation, which most likely are inaccurate\n", ndefault);
463 }
464 }
465 else
466 {
467 vp = NULL;
468 }
469
470 gmx_atomprop_destroy(aps);
471
472 if (bPBC)
473 {
474 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms, box);
475 }
476
477 nfr = 0;
478 do
479 {
480 if (bPBC)
481 {
482 gmx_rmpbc(gpbc, natoms, box, x);
483 }
484
485 bConnelly = (nfr == 0 && opt2bSet("-q", nfile, fnm));
486 if (bConnelly)
487 {
488 if (!bTop)
489 {
490 gmx_fatal(FARGS, "Need a tpr file for Connelly plot");
491 }
492 flag = FLAG_ATOM_AREA | FLAG_DOTS;
493 }
494 else
495 {
496 flag = FLAG_ATOM_AREA;
497 }
498 if (vp)
499 {
500 flag = flag | FLAG_VOLUME;
501 }
502
503 if (debug)
504 {
505 write_sto_conf("check.pdb", "pbc check", atoms, x, NULL, ePBC, box);
506 }
507
508 retval = nsc_dclm_pbc(x, radius, nx[0], ndots, flag, &totarea,
509 &area, &totvolume, &surfacedots, &nsurfacedots,
510 index[0], ePBC, bPBC ? box : NULL);
511 if (retval)
512 {
513 gmx_fatal(FARGS, "Something wrong in nsc_dclm_pbc");
514 }
515
516 if (bConnelly)
517 {
518 connelly_plot(ftp2fn(efPDB, nfile, fnm),
519 nsurfacedots, surfacedots, x, atoms,
520 &(top.symtab), ePBC, box, bSave);
521 }
522 harea = 0;
523 tarea = 0;
524 dgsolv = 0;
525 if (bResAt)
526 {
527 for (i = 0; i < atoms->nres; i++)
528 {
529 res_a[i] = 0;
530 }
531 }
532 for (i = 0; (i < nx[0]); i++)
533 {
534 ii = index[0][i];
535 if (bOut[ii])
536 {
537 at_area = area[i];
538 if (bResAt)
539 {
540 atom_area[i] += at_area;
541 atom_area2[i] += sqr(at_area);
542 res_a[atoms->atom[ii].resind] += at_area;
543 }
544 tarea += at_area;
545 if (bDGsol)
546 {
547 dgsolv += at_area*dgs_factor[i];
548 }
549 if (bPhobic[i])
550 {
551 harea += at_area;
552 }
553 }
554 }
555 if (bResAt)
556 {
557 for (i = 0; i < atoms->nres; i++)
558 {
559 res_area[i] += res_a[i];
560 res_area2[i] += sqr(res_a[i]);
561 }
562 }
563 fprintf(fp, "%10g %10g %10g %10g", t, harea, tarea-harea, tarea);
564 if (bDGsol)
565 {
566 fprintf(fp, " %10g\n", dgsolv);
567 }
568 else
569 {
570 fprintf(fp, "\n");
571 }
572
573 /* Print volume */
574 if (vp)
575 {
576 density = totmass*AMU/(totvolume*NANO*NANO*NANO);
577 fprintf(vp, "%12.5e %12.5e %12.5e\n", t, totvolume, density);
578 }
579 if (area)
580 {
581 sfree(area);
582 area = NULL;
583 }
584 if (surfacedots)
585 {
586 sfree(surfacedots);
587 surfacedots = NULL;
588 }
589 nfr++;
590 }
591 while (read_next_x(oenv, status, &t, natoms, x, box));
592
593 if (bPBC)
594 {
595 gmx_rmpbc_done(gpbc);
596 }
597
598 fprintf(stderr, "\n");
599 close_trj(status);
600 ffclose(fp);
601 if (vp)
602 {
603 ffclose(vp);
604 }
605
606 /* if necessary, print areas per atom to file too: */
607 if (bResAt)
608 {
609 for (i = 0; i < atoms->nres; i++)
610 {
611 res_area[i] /= nfr;
612 res_area2[i] /= nfr;
613 }
614 for (i = 0; i < nx[0]; i++)
615 {
616 atom_area[i] /= nfr;
617 atom_area2[i] /= nfr;
618 }
619 fprintf(stderr, "Printing out areas per atom\n");
620 fp = xvgropen(opt2fn("-or", nfile, fnm), "Area per residue over the trajectory", "Residue",
621 "Area (nm\\S2\\N)", oenv);
622 xvgr_legend(fp, asize(or_and_oa_legend), or_and_oa_legend, oenv);
623 fp2 = xvgropen(opt2fn("-oa", nfile, fnm), "Area per atom over the trajectory", "Atom #",
624 "Area (nm\\S2\\N)", oenv);
625 xvgr_legend(fp2, asize(or_and_oa_legend), or_and_oa_legend, oenv);
626 if (bITP)
627 {
628 fp3 = ftp2FILE(efITP, nfile, fnm, "w");
629 fprintf(fp3, "[ position_restraints ]\n"
630 "#define FCX 1000\n"
631 "#define FCY 1000\n"
632 "#define FCZ 1000\n"
633 "; Atom Type fx fy fz\n");
634 }
635 for (i = 0; i < nx[0]; i++)
636 {
637 ii = index[0][i];
638 res = atoms->atom[ii].resind;
639 if (i == nx[0]-1 || res != atoms->atom[index[0][i+1]].resind)
640 {
641 fluc2 = res_area2[res]-sqr(res_area[res]);
642 if (fluc2 < 0)
643 {
644 fluc2 = 0;
645 }
646 fprintf(fp, "%10d %10g %10g\n",
647 atoms->resinfo[res].nr, res_area[res], sqrt(fluc2));
648 }
649 fluc2 = atom_area2[i]-sqr(atom_area[i]);
650 if (fluc2 < 0)
651 {
652 fluc2 = 0;
653 }
654 fprintf(fp2, "%d %g %g\n", index[0][i]+1, atom_area[i], sqrt(fluc2));
655 if (bITP && (atom_area[i] > minarea))
656 {
657 fprintf(fp3, "%5d 1 FCX FCX FCZ\n", ii+1);
658 }
659 }
660 if (bITP)
661 {
662 ffclose(fp3);
663 }
664 ffclose(fp);
665 }
666
667 /* Be a good citizen, keep our memory free! */
668 sfree(x);
669 sfree(nx);
670 for (i = 0; i < 2; i++)
671 {
672 sfree(index[i]);
673 sfree(grpname[i]);
674 }
675 sfree(bOut);
676 sfree(radius);
677 sfree(bPhobic);
678
679 if (bResAt)
680 {
681 sfree(atom_area);
682 sfree(atom_area2);
683 sfree(res_a);
684 sfree(res_area);
685 sfree(res_area2);
686 }
687 if (bDGsol)
688 {
689 sfree(dgs_factor);
690 }
691 }
692
693 int gmx_sas(int argc, char *argv[])
694 {
695 const char *desc[] = {
696 "[TT]g_sas[tt] computes hydrophobic, hydrophilic and total solvent",
697 "accessible surface area. See Eisenhaber F, Lijnzaad P, Argos P,",
698 "Sander C, & Scharf M (1995) J. Comput. Chem. 16, 273-284.",
699 "As a side effect, the Connolly surface can be generated as well in",
700 "a [TT].pdb[tt] file where the nodes are represented as atoms and the",
701 "vertice connecting the nearest nodes as CONECT records.",
702 "The program will ask for a group for the surface calculation",
703 "and a group for the output. The calculation group should always",
704 "consists of all the non-solvent atoms in the system.",
705 "The output group can be the whole or part of the calculation group.",
706 "The average and standard deviation of the area over the trajectory can be plotted",
707 "per residue and atom as well (options [TT]-or[tt] and [TT]-oa[tt]).",
708 "In combination with the latter option an [TT].itp[tt] file can be",
709 "generated (option [TT]-i[tt])",
710 "which can be used to restrain surface atoms.[PAR]",
711
712 "By default, periodic boundary conditions are taken into account,",
713 "this can be turned off using the [TT]-nopbc[tt] option.[PAR]",
714
715 "With the [TT]-tv[tt] option the total volume and density of the",
716 "molecule can be computed.",
717 "Please consider whether the normal probe radius is appropriate",
718 "in this case or whether you would rather use e.g. 0. It is good",
719 "to keep in mind that the results for volume and density are very",
720 "approximate. For example, in ice Ih, one can easily fit water molecules in the",
721 "pores which would yield a volume that is too low, and surface area and density",
722 "that are both too high."
723 };
724
725 output_env_t oenv;
726 static real solsize = 0.14;
727 static int ndots = 24;
728 static real qcut = 0.2;
729 static real minarea = 0.5, dgs_default = 0;
730 static gmx_bool bSave = TRUE, bPBC = TRUE, bFindex = FALSE;
731 t_pargs pa[] = {
732 { "-probe", FALSE, etREAL, {&solsize},
733 "Radius of the solvent probe (nm)" },
734 { "-ndots", FALSE, etINT, {&ndots},
735 "Number of dots per sphere, more dots means more accuracy" },
736 { "-qmax", FALSE, etREAL, {&qcut},
737 "The maximum charge (e, absolute value) of a hydrophobic atom" },
738 { "-f_index", FALSE, etBOOL, {&bFindex},
739 "Determine from a group in the index file what are the hydrophobic atoms rather than from the charge" },
740 { "-minarea", FALSE, etREAL, {&minarea},
741 "The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)" },
742 { "-pbc", FALSE, etBOOL, {&bPBC},
743 "Take periodicity into account" },
744 { "-prot", FALSE, etBOOL, {&bSave},
745 "Output the protein to the Connelly [TT].pdb[tt] file too" },
746 { "-dgs", FALSE, etREAL, {&dgs_default},
747 "Default value for solvation free energy per area (kJ/mol/nm^2)" }
748 };
749 t_filenm fnm[] = {
750 { efTRX, "-f", NULL, ffREAD },
751 { efTPS, "-s", NULL, ffREAD },
752 { efXVG, "-o", "area", ffWRITE },
753 { efXVG, "-or", "resarea", ffOPTWR },
754 { efXVG, "-oa", "atomarea", ffOPTWR },
755 { efXVG, "-tv", "volume", ffOPTWR },
756 { efPDB, "-q", "connelly", ffOPTWR },
757 { efNDX, "-n", "index", ffOPTRD },
758 { efITP, "-i", "surfat", ffOPTWR }
759 };
760 #define NFILE asize(fnm)
761
762 CopyRight(stderr, argv[0]);
763 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
764 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
765 if (solsize < 0)
766 {
767 solsize = 1e-3;
768 fprintf(stderr, "Probe size too small, setting it to %g\n", solsize);
769 }
770 if (ndots < 20)
771 {
772 ndots = 20;
773 fprintf(stderr, "Ndots too small, setting it to %d\n", ndots);
774 }
775
776 please_cite(stderr, "Eisenhaber95");
777
778 sas_plot(NFILE, fnm, solsize, ndots, qcut, bSave, minarea, bPBC, dgs_default, bFindex,
779 oenv);
780
781 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
782 do_view(oenv, opt2fn_null("-or", NFILE, fnm), "-nxy");
783 do_view(oenv, opt2fn_null("-oa", NFILE, fnm), "-nxy");
784
785 thanx(stderr);
786
787 return 0;
788 }