Update instructions in containers.rst
[gromacs.git] / src / gromacs / gmxana / gmx_do_dssp.cpp
blob1cfcde6f3876650d12fe9f2c2f8951bbd8cacf8f
1 /*
2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
38 #include "gmxpre.h"
40 #include "config.h"
42 #include <cstdlib>
43 #include <cstring>
45 #include <algorithm>
46 #include <numeric>
48 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/commandline/viewit.h"
50 #include "gromacs/fileio/confio.h"
51 #include "gromacs/fileio/matio.h"
52 #include "gromacs/fileio/pdbio.h"
53 #include "gromacs/fileio/trxio.h"
54 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/gmxana/gmx_ana.h"
56 #include "gromacs/gmxana/gstat.h"
57 #include "gromacs/pbcutil/rmpbc.h"
58 #include "gromacs/topology/index.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/utility/arraysize.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/dir_separator.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/smalloc.h"
66 #include "gromacs/utility/strdb.h"
68 #if GMX_NATIVE_WINDOWS
69 # define NULL_DEVICE "nul"
70 # define popen _popen
71 # define pclose _pclose
72 #else
73 # define NULL_DEVICE "/dev/null"
74 #endif
76 struct DsspInputStrings
78 std::string dptr;
79 std::string pdbfile;
80 std::string tmpfile;
83 static const char* prepareToRedirectStdout(bool bVerbose)
85 return bVerbose ? "" : "2>" NULL_DEVICE;
88 static void printDsspResult(char* dssp, const DsspInputStrings& strings, const std::string& redirectionString)
90 #if HAVE_PIPES || GMX_NATIVE_WINDOWS
91 sprintf(dssp, "%s -i %s %s", strings.dptr.c_str(), strings.pdbfile.c_str(), redirectionString.c_str());
92 #else
93 sprintf(dssp, "%s -i %s -o %s > %s %s", strings.dptr.c_str(), strings.pdbfile.c_str(),
94 strings.tmpfile.c_str(), NULL_DEVICE, redirectionString.c_str());
95 #endif
99 static int strip_dssp(FILE* tapeout,
100 int nres,
101 const gmx_bool bPhobres[],
102 real t,
103 real* acc,
104 FILE* fTArea,
105 t_matrix* mat,
106 int average_area[],
107 const gmx_output_env_t* oenv)
109 static gmx_bool bFirst = TRUE;
110 static char* ssbuf;
111 char buf[STRLEN + 1];
112 char SSTP;
113 int nr, iacc, nresidues;
114 int naccf, naccb; /* Count hydrophobic and hydrophilic residues */
115 real iaccf, iaccb;
118 /* Skip header */
121 fgets2(buf, STRLEN, tapeout);
122 } while (std::strstr(buf, "KAPPA") == nullptr);
123 if (bFirst)
125 /* Since we also have empty lines in the dssp output (temp) file,
126 * and some line content is saved to the ssbuf variable,
127 * we need more memory than just nres elements. To be shure,
128 * we allocate 2*nres-1, since for each chain there is a
129 * separating line in the temp file. (At most each residue
130 * could have been defined as a separate chain.) */
131 snew(ssbuf, 2 * nres - 1);
134 iaccb = iaccf = 0;
135 nresidues = 0;
136 naccf = 0;
137 naccb = 0;
138 for (nr = 0; (fgets2(buf, STRLEN, tapeout) != nullptr); nr++)
140 if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
142 SSTP = '='; /* Chain separator sign '=' */
144 else
146 SSTP = buf[16] == ' ' ? '~' : buf[16];
148 ssbuf[nr] = SSTP;
150 buf[39] = '\0';
152 /* Only calculate solvent accessible area if needed */
153 if ((nullptr != acc) && (buf[13] != '!'))
155 sscanf(&(buf[34]), "%d", &iacc);
156 acc[nr] = iacc;
157 /* average_area and bPhobres are counted from 0...nres-1 */
158 average_area[nresidues] += iacc;
159 if (bPhobres[nresidues])
161 naccb++;
162 iaccb += iacc;
164 else
166 naccf++;
167 iaccf += iacc;
169 /* Keep track of the residue number (do not count chain separator lines '!') */
170 nresidues++;
173 ssbuf[nr] = '\0';
175 if (bFirst)
177 if (nullptr != acc)
179 fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n",
180 naccb, naccf);
183 mat->title = "Secondary structure";
184 mat->legend = "";
185 mat->label_x = output_env_get_time_label(oenv);
186 mat->label_y = "Residue";
187 mat->bDiscrete = true;
188 mat->ny = nr;
189 mat->axis_y.resize(nr);
190 std::iota(mat->axis_y.begin(), mat->axis_y.end(), 1);
191 mat->axis_x.resize(0);
192 mat->matrix.resize(1, 1);
193 bFirst = false;
195 mat->axis_x.push_back(t);
196 mat->matrix.resize(++(mat->nx), nr);
197 auto columnIndex = mat->nx - 1;
198 for (int i = 0; i < nr; i++)
200 t_xpmelmt c = { ssbuf[i], 0 };
201 mat->matrix(columnIndex, i) = std::max(static_cast<t_matelmt>(0), searchcmap(mat->map, c));
204 if (fTArea)
206 fprintf(fTArea, "%10g %10g %10g\n", t, 0.01 * iaccb, 0.01 * iaccf);
209 /* Return the number of lines found in the dssp file (i.e. number
210 * of redidues plus chain separator lines).
211 * This is the number of y elements needed for the area xpm file */
212 return nr;
215 static gmx_bool* bPhobics(t_atoms* atoms)
217 int j, i, nb;
218 char** cb;
219 gmx_bool* bb;
220 int n_surf;
221 char surffn[] = "surface.dat";
222 char ** surf_res, **surf_lines;
225 nb = get_lines("phbres.dat", &cb);
226 snew(bb, atoms->nres);
228 n_surf = get_lines(surffn, &surf_lines);
229 snew(surf_res, n_surf);
230 for (i = 0; (i < n_surf); i++)
232 snew(surf_res[i], 5);
233 sscanf(surf_lines[i], "%s", surf_res[i]);
237 for (i = 0, j = 0; (i < atoms->nres); i++)
239 if (-1 != search_str(n_surf, surf_res, *atoms->resinfo[i].name))
241 bb[j++] = (-1 != search_str(nb, cb, *atoms->resinfo[i].name));
245 if (i != j)
247 fprintf(stderr,
248 "Not all residues were recognized (%d from %d), the result may be inaccurate!\n", j, i);
251 for (i = 0; (i < n_surf); i++)
253 sfree(surf_res[i]);
255 sfree(surf_res);
257 return bb;
260 static void check_oo(t_atoms* atoms)
262 char* OOO = gmx_strdup("O");
264 for (int i = 0; (i < atoms->nr); i++)
266 if ((std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
267 || (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
268 || (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
269 || (std::strcmp(*(atoms->atomname[i]), "OT1") == 0))
271 *atoms->atomname[i] = OOO;
276 static void norm_acc(t_atoms* atoms, int nres, const real av_area[], real norm_av_area[])
278 int i, n, n_surf;
280 char surffn[] = "surface.dat";
281 char ** surf_res, **surf_lines;
282 double* surf;
284 n_surf = get_lines(surffn, &surf_lines);
285 snew(surf, n_surf);
286 snew(surf_res, n_surf);
287 for (i = 0; (i < n_surf); i++)
289 snew(surf_res[i], 5);
290 sscanf(surf_lines[i], "%s %lf", surf_res[i], &surf[i]);
293 for (i = 0; (i < nres); i++)
295 n = search_str(n_surf, surf_res, *atoms->resinfo[i].name);
296 if (n != -1)
298 norm_av_area[i] = av_area[i] / surf[n];
300 else
302 fprintf(stderr, "Residue %s not found in surface database (%s)\n",
303 *atoms->resinfo[i].name, surffn);
308 static void prune_ss_legend(t_matrix* mat)
310 std::vector<bool> isPresent(mat->map.size());
311 std::vector<int> newnum(mat->map.size());
312 std::vector<t_mapping> newmap;
314 for (int f = 0; f < mat->nx; f++)
316 for (int r = 0; r < mat->ny; r++)
318 isPresent[mat->matrix(f, r)] = true;
322 for (size_t i = 0; i < mat->map.size(); i++)
324 newnum[i] = -1;
325 if (isPresent[i])
327 newnum[i] = newmap.size();
328 newmap.emplace_back(mat->map[i]);
331 if (newmap.size() != mat->map.size())
333 std::swap(mat->map, newmap);
334 for (int f = 0; f < mat->nx; f++)
336 for (int r = 0; r < mat->ny; r++)
338 mat->matrix(f, r) = newnum[mat->matrix(f, r)];
344 static void write_sas_mat(const char* fn, real** accr, int nframe, int nres, t_matrix* mat)
346 real lo, hi;
347 int i, j, nlev;
348 t_rgb rlo = { 1, 1, 1 }, rhi = { 0, 0, 0 };
349 FILE* fp;
351 if (fn)
353 hi = lo = accr[0][0];
354 for (i = 0; i < nframe; i++)
356 for (j = 0; j < nres; j++)
358 lo = std::min(lo, accr[i][j]);
359 hi = std::max(hi, accr[i][j]);
362 fp = gmx_ffopen(fn, "w");
363 nlev = static_cast<int>(hi - lo + 1);
364 write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)", "Time", "Residue Index",
365 nframe, nres, mat->axis_x.data(), mat->axis_y.data(), accr, lo, hi, rlo, rhi, &nlev);
366 gmx_ffclose(fp);
370 static void analyse_ss(const char* outfile, t_matrix* mat, const char* ss_string, const gmx_output_env_t* oenv)
372 FILE* fp;
373 int ss_count, total_count;
375 gmx::ArrayRef<t_mapping> map = mat->map;
376 std::vector<int> count(map.size());
377 std::vector<int> total(map.size(), 0);
378 // This copying would not be necessary if xvgr_legend could take a
379 // view of string views
380 std::vector<std::string> leg;
381 leg.reserve(map.size() + 1);
382 leg.emplace_back("Structure");
383 for (const auto& m : map)
385 leg.emplace_back(m.desc);
388 fp = xvgropen(outfile, "Secondary Structure", output_env_get_xvgr_tlabel(oenv),
389 "Number of Residues", oenv);
390 if (output_env_get_print_xvgr_codes(oenv))
392 fprintf(fp, "@ subtitle \"Structure = ");
394 for (size_t s = 0; s < std::strlen(ss_string); s++)
396 if (s > 0)
398 fprintf(fp, " + ");
400 for (const auto& m : map)
402 if (ss_string[s] == m.code.c1)
404 fprintf(fp, "%s", m.desc);
408 fprintf(fp, "\"\n");
409 xvgrLegend(fp, leg, oenv);
411 total_count = 0;
412 for (int f = 0; f < mat->nx; f++)
414 ss_count = 0;
415 for (auto& c : count)
417 c = 0;
419 for (int r = 0; r < mat->ny; r++)
421 count[mat->matrix(f, r)]++;
422 total[mat->matrix(f, r)]++;
424 for (gmx::index s = 0; s != gmx::ssize(map); ++s)
426 if (std::strchr(ss_string, map[s].code.c1))
428 ss_count += count[s];
429 total_count += count[s];
432 fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
433 for (const auto& c : count)
435 fprintf(fp, " %5d", c);
437 fprintf(fp, "\n");
439 /* now print column totals */
440 fprintf(fp, "%-8s %5d", "# Totals", total_count);
441 for (const auto& t : total)
443 fprintf(fp, " %5d", t);
445 fprintf(fp, "\n");
447 /* now print probabilities */
448 fprintf(fp, "%-8s %5.2f", "# SS pr.", total_count / static_cast<real>(mat->nx * mat->ny));
449 for (const auto& t : total)
451 fprintf(fp, " %5.2f", t / static_cast<real>(mat->nx * mat->ny));
453 fprintf(fp, "\n");
455 xvgrclose(fp);
458 int gmx_do_dssp(int argc, char* argv[])
460 const char* desc[] = {
461 "[THISMODULE] ", "reads a trajectory file and computes the secondary structure for",
462 "each time frame ", "calling the dssp program. If you do not have the dssp program,",
463 "get it from http://swift.cmbi.ru.nl/gv/dssp. [THISMODULE] assumes ",
464 "that the dssp executable is located in ",
465 // NOLINTNEXTLINE(bugprone-suspicious-missing-comma)
466 "[TT]" GMX_DSSP_PROGRAM_PATH "[tt]. If this is not the case, then you should",
467 "set an environment variable [TT]DSSP[tt] pointing to the dssp", "executable, e.g.: [PAR]",
468 "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
469 "Since version 2.0.0, dssp is invoked with a syntax that differs",
470 "from earlier versions. If you have an older version of dssp,",
471 "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
472 "By default, do_dssp uses the syntax introduced with version 2.0.0.",
473 "Even newer versions (which at the time of writing are not yet released)",
474 "are assumed to have the same syntax as 2.0.0.[PAR]",
475 "The structure assignment for each residue and time is written to an",
476 "[REF].xpm[ref] matrix file. This file can be visualized with for instance",
477 "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
478 "Individual chains are separated by light grey lines in the [REF].xpm[ref] and",
479 "postscript files.", "The number of residues with each secondary structure type and the",
480 "total secondary structure ([TT]-sss[tt]) count as a function of",
481 "time are also written to file ([TT]-sc[tt]).[PAR]",
482 "Solvent accessible surface (SAS) per residue can be calculated, both in",
483 "absolute values (A^2) and in fractions of the maximal accessible",
484 "surface of a residue. The maximal accessible surface is defined as",
485 "the accessible surface of a residue in a chain of glycines.",
486 "[BB]Note[bb] that the program [gmx-sas] can also compute SAS",
487 "and that is more efficient.[PAR]",
488 "Finally, this program can dump the secondary structure in a special file",
489 "[TT]ssdump.dat[tt] for usage in the program [gmx-chi]. Together",
490 "these two programs can be used to analyze dihedral properties as a",
491 "function of secondary structure type."
493 static gmx_bool bVerbose;
494 static const char* ss_string = "HEBT";
495 static int dsspVersion = 2;
496 t_pargs pa[] = {
497 { "-v", FALSE, etBOOL, { &bVerbose }, "HIDDENGenerate miles of useless information" },
498 { "-sss", FALSE, etSTR, { &ss_string }, "Secondary structures for structure count" },
499 { "-ver",
500 FALSE,
501 etINT,
502 { &dsspVersion },
503 "DSSP major version. Syntax changed with version 2" }
506 t_trxstatus* status;
507 FILE * tapein, *tapeout;
508 FILE * ss, *acc, *fTArea, *tmpf;
509 const char * fnSCount, *fnArea, *fnTArea, *fnAArea;
510 const char* leg[] = { "Phobic", "Phylic" };
511 t_topology top;
512 PbcType pbcType;
513 t_atoms* atoms;
514 t_matrix mat;
515 int nres, nr0, naccr, nres_plus_separators;
516 gmx_bool * bPhbres, bDoAccSurf;
517 real t;
518 int natoms, nframe = 0;
519 matrix box = { { 0 } };
520 int gnx;
521 char* grpnm;
522 int* index;
523 rvec * xp, *x;
524 int* average_area;
525 real ** accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
526 char pdbfile[32], tmpfile[32];
527 char dssp[256];
528 const char* dptr;
529 gmx_output_env_t* oenv;
530 gmx_rmpbc_t gpbc = nullptr;
532 t_filenm fnm[] = {
533 { efTRX, "-f", nullptr, ffREAD }, { efTPS, nullptr, nullptr, ffREAD },
534 { efNDX, nullptr, nullptr, ffOPTRD }, { efDAT, "-ssdump", "ssdump", ffOPTWR },
535 { efMAP, "-map", "ss", ffLIBRD }, { efXPM, "-o", "ss", ffWRITE },
536 { efXVG, "-sc", "scount", ffWRITE }, { efXPM, "-a", "area", ffOPTWR },
537 { efXVG, "-ta", "totarea", ffOPTWR }, { efXVG, "-aa", "averarea", ffOPTWR }
539 #define NFILE asize(fnm)
541 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT, NFILE, fnm,
542 asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
544 return 0;
546 fnSCount = opt2fn("-sc", NFILE, fnm);
547 fnArea = opt2fn_null("-a", NFILE, fnm);
548 fnTArea = opt2fn_null("-ta", NFILE, fnm);
549 fnAArea = opt2fn_null("-aa", NFILE, fnm);
550 bDoAccSurf = ((fnArea != nullptr) || (fnTArea != nullptr) || (fnAArea != nullptr));
552 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xp, nullptr, box, FALSE);
553 atoms = &(top.atoms);
554 check_oo(atoms);
555 bPhbres = bPhobics(atoms);
557 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpnm);
558 nres = 0;
559 nr0 = -1;
560 for (int i = 0; (i < gnx); i++)
562 if (atoms->atom[index[i]].resind != nr0)
564 nr0 = atoms->atom[index[i]].resind;
565 nres++;
568 fprintf(stderr, "There are %d residues in your selected group\n", nres);
570 std::strcpy(pdbfile, "ddXXXXXX");
571 gmx_tmpnam(pdbfile);
572 if ((tmpf = fopen(pdbfile, "w")) == nullptr)
574 sprintf(pdbfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
575 gmx_tmpnam(pdbfile);
576 if ((tmpf = fopen(pdbfile, "w")) == nullptr)
578 gmx_fatal(FARGS, "Can not open tmp file %s", pdbfile);
581 fclose(tmpf);
583 std::strcpy(tmpfile, "ddXXXXXX");
584 gmx_tmpnam(tmpfile);
585 if ((tmpf = fopen(tmpfile, "w")) == nullptr)
587 sprintf(tmpfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
588 gmx_tmpnam(tmpfile);
589 if ((tmpf = fopen(tmpfile, "w")) == nullptr)
591 gmx_fatal(FARGS, "Can not open tmp file %s", tmpfile);
594 fclose(tmpf);
596 const std::string defpathenv = GMX_DSSP_PROGRAM_PATH;
597 if ((dptr = getenv("DSSP")) == nullptr)
599 dptr = defpathenv.c_str();
601 if (!gmx_fexist(dptr))
603 gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)", dptr);
605 std::string redirectionString;
606 redirectionString = prepareToRedirectStdout(bVerbose);
607 DsspInputStrings dsspStrings;
608 dsspStrings.dptr = dptr;
609 dsspStrings.pdbfile = pdbfile;
610 dsspStrings.tmpfile = tmpfile;
611 if (dsspVersion >= 2)
613 if (dsspVersion > 2)
615 printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by "
616 "do_dssp. Assuming version 2 syntax.\n\n",
617 dsspVersion);
620 printDsspResult(dssp, dsspStrings, redirectionString);
622 else
624 if (bDoAccSurf)
626 dsspStrings.dptr.clear();
628 else
630 dsspStrings.dptr = "-na";
632 printDsspResult(dssp, dsspStrings, redirectionString);
634 fprintf(stderr, "dssp cmd='%s'\n", dssp);
636 if (fnTArea)
638 fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",
639 output_env_get_xvgr_tlabel(oenv), "Area (nm\\S2\\N)", oenv);
640 xvgr_legend(fTArea, 2, leg, oenv);
642 else
644 fTArea = nullptr;
647 mat.map = readcmap(opt2fn("-map", NFILE, fnm));
649 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
650 if (natoms > atoms->nr)
652 gmx_fatal(FARGS, "\nTrajectory does not match topology!");
654 if (gnx > natoms)
656 gmx_fatal(FARGS, "\nTrajectory does not match selected group!");
659 snew(average_area, atoms->nres);
660 snew(av_area, atoms->nres);
661 snew(norm_av_area, atoms->nres);
662 accr = nullptr;
663 naccr = 0;
665 gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
668 t = output_env_conv_time(oenv, t);
669 if (bDoAccSurf && nframe >= naccr)
671 naccr += 10;
672 srenew(accr, naccr);
673 for (int i = naccr - 10; i < naccr; i++)
675 snew(accr[i], 2 * atoms->nres - 1);
678 gmx_rmpbc(gpbc, natoms, box, x);
679 tapein = gmx_ffopen(pdbfile, "w");
680 write_pdbfile_indexed(tapein, nullptr, atoms, x, pbcType, box, ' ', -1, gnx, index, nullptr, FALSE);
681 gmx_ffclose(tapein);
682 /* strip_dssp returns the number of lines found in the dssp file, i.e.
683 * the number of residues plus the separator lines */
685 #if HAVE_PIPES || GMX_NATIVE_WINDOWS
686 if (nullptr == (tapeout = popen(dssp, "r")))
687 #else
688 if (0 != system(dssp) || nullptr == (tapeout = gmx_ffopen(tmpfile, "r")))
689 #endif
691 remove(pdbfile);
692 remove(tmpfile);
693 gmx_fatal(FARGS,
694 "Failed to execute command: %s\n"
695 "Try specifying your dssp version with the -ver option.",
696 dssp);
698 if (bDoAccSurf)
700 accr_ptr = accr[nframe];
702 /* strip_dssp returns the number of lines found in the dssp file, i.e.
703 * the number of residues plus the separator lines */
704 nres_plus_separators =
705 strip_dssp(tapeout, nres, bPhbres, t, accr_ptr, fTArea, &mat, average_area, oenv);
706 #if HAVE_PIPES || GMX_NATIVE_WINDOWS
707 pclose(tapeout);
708 #else
709 gmx_ffclose(tapeout);
710 #endif
711 remove(tmpfile);
712 remove(pdbfile);
713 nframe++;
714 } while (read_next_x(oenv, status, &t, x, box));
715 fprintf(stderr, "\n");
716 close_trx(status);
717 if (fTArea)
719 xvgrclose(fTArea);
721 gmx_rmpbc_done(gpbc);
723 prune_ss_legend(&mat);
725 ss = opt2FILE("-o", NFILE, fnm, "w");
726 mat.flags = 0;
727 write_xpm_m(ss, mat);
728 gmx_ffclose(ss);
730 if (opt2bSet("-ssdump", NFILE, fnm))
732 ss = opt2FILE("-ssdump", NFILE, fnm, "w");
733 fprintf(ss, "%d\n", nres);
734 for (gmx::index j = 0; j != mat.matrix.extent(0); ++j)
736 auto row = mat.matrix.asView()[j];
737 for (gmx::index i = 0; i != row.extent(0); ++i)
739 fputc(mat.map[row[i]].code.c1, ss);
741 fputc('\n', ss);
743 gmx_ffclose(ss);
745 analyse_ss(fnSCount, &mat, ss_string, oenv);
747 if (bDoAccSurf)
749 write_sas_mat(fnArea, accr, nframe, nres_plus_separators, &mat);
751 for (int i = 0; i < atoms->nres; i++)
753 av_area[i] = (average_area[i] / static_cast<real>(nframe));
756 norm_acc(atoms, nres, av_area, norm_av_area);
758 if (fnAArea)
760 acc = xvgropen(fnAArea, "Average Accessible Area", "Residue", "A\\S2", oenv);
761 for (int i = 0; (i < nres); i++)
763 fprintf(acc, "%5d %10g %10g\n", i + 1, av_area[i], norm_av_area[i]);
765 xvgrclose(acc);
769 view_all(oenv, NFILE, fnm);
771 return 0;