*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2013,2014,2015,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/utility/strdb.h"
#if GMX_NATIVE_WINDOWS
- #define NULL_DEVICE "nul"
- #define popen _popen
- #define pclose _pclose
+# define NULL_DEVICE "nul"
+# define popen _popen
+# define pclose _pclose
#else
- #define NULL_DEVICE "/dev/null"
+# define NULL_DEVICE "/dev/null"
#endif
struct DsspInputStrings
return bVerbose ? "" : "2>" NULL_DEVICE;
}
-static void printDsspResult(char *dssp, const DsspInputStrings &strings,
- const std::string &redirectionString)
+static void printDsspResult(char* dssp, const DsspInputStrings& strings, const std::string& redirectionString)
{
#if HAVE_PIPES || GMX_NATIVE_WINDOWS
- sprintf(dssp, "%s -i %s %s",
- strings.dptr.c_str(), strings.pdbfile.c_str(), redirectionString.c_str());
+ sprintf(dssp, "%s -i %s %s", strings.dptr.c_str(), strings.pdbfile.c_str(), redirectionString.c_str());
#else
- sprintf(dssp, "%s -i %s -o %s > %s %s",
- strings.dptr.c_str(), strings.pdbfile.c_str(), strings.tmpfile.c_str(), NULL_DEVICE, redirectionString.c_str());
+ sprintf(dssp,
+ "%s -i %s -o %s > %s %s",
+ strings.dptr.c_str(),
+ strings.pdbfile.c_str(),
+ strings.tmpfile.c_str(),
+ NULL_DEVICE,
+ redirectionString.c_str());
#endif
}
-static int strip_dssp(FILE *tapeout, int nres,
- const gmx_bool bPhobres[], real t,
- real *acc, FILE *fTArea,
- t_matrix *mat, int average_area[],
- const gmx_output_env_t *oenv)
+static int strip_dssp(FILE* tapeout,
+ int nres,
+ const gmx_bool bPhobres[],
+ real t,
+ real* acc,
+ FILE* fTArea,
+ t_matrix* mat,
+ int average_area[],
+ const gmx_output_env_t* oenv)
{
- static gmx_bool bFirst = TRUE;
- static char *ssbuf;
- char buf[STRLEN+1];
- char SSTP;
- int nr, iacc, nresidues;
- int naccf, naccb; /* Count hydrophobic and hydrophilic residues */
- real iaccf, iaccb;
+ static gmx_bool bFirst = TRUE;
+ static std::string ssbuf;
+ char buf[STRLEN + 1];
+ char SSTP;
+ int nr, iacc, nresidues;
+ int naccf, naccb; /* Count hydrophobic and hydrophilic residues */
+ real iaccf, iaccb;
/* Skip header */
do
{
fgets2(buf, STRLEN, tapeout);
- }
- while (std::strstr(buf, "KAPPA") == nullptr);
+ } while (std::strstr(buf, "KAPPA") == nullptr);
if (bFirst)
{
/* Since we also have empty lines in the dssp output (temp) file,
* we allocate 2*nres-1, since for each chain there is a
* separating line in the temp file. (At most each residue
* could have been defined as a separate chain.) */
- snew(ssbuf, 2*nres-1);
+ ssbuf.resize(2 * nres - 1);
}
- iaccb = iaccf = 0;
- nresidues = 0;
- naccf = 0;
- naccb = 0;
+ iaccb = iaccf = 0;
+ nresidues = 0;
+ naccf = 0;
+ naccb = 0;
for (nr = 0; (fgets2(buf, STRLEN, tapeout) != nullptr); nr++)
{
if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
{
- SSTP = '='; /* Chain separator sign '=' */
+ SSTP = '='; /* Chain separator sign '=' */
}
else
{
mat->axis_y.resize(nr);
std::iota(mat->axis_y.begin(), mat->axis_y.end(), 1);
mat->axis_x.resize(0);
- mat->matrix.resize(0, 0);
- bFirst = false;
+ mat->matrix.resize(1, 1);
+ bFirst = false;
}
mat->axis_x.push_back(t);
- mat->matrix.resize(mat->matrix.extent(0), nr);
- mat->nx = mat->matrix.extent(0);
+ mat->matrix.resize(++(mat->nx), nr);
auto columnIndex = mat->nx - 1;
for (int i = 0; i < nr; i++)
{
- t_xpmelmt c = {ssbuf[i], 0};
+ t_xpmelmt c = { ssbuf[i], 0 };
mat->matrix(columnIndex, i) = std::max(static_cast<t_matelmt>(0), searchcmap(mat->map, c));
}
if (fTArea)
{
- fprintf(fTArea, "%10g %10g %10g\n", t, 0.01*iaccb, 0.01*iaccf);
+ fprintf(fTArea, "%10g %10g %10g\n", t, 0.01 * iaccb, 0.01 * iaccf);
}
/* Return the number of lines found in the dssp file (i.e. number
return nr;
}
-static gmx_bool *bPhobics(t_atoms *atoms)
+static gmx_bool* bPhobics(t_atoms* atoms)
{
- int j, i, nb;
- char **cb;
- gmx_bool *bb;
- int n_surf;
- char surffn[] = "surface.dat";
- char **surf_res, **surf_lines;
+ int j, i, nb;
+ char** cb;
+ gmx_bool* bb;
+ int n_surf;
+ char surffn[] = "surface.dat";
+ char ** surf_res, **surf_lines;
nb = get_lines("phbres.dat", &cb);
for (i = 0, j = 0; (i < atoms->nres); i++)
{
- if (-1 != search_str(n_surf, surf_res, *atoms->resinfo[i].name) )
+ if (-1 != search_str(n_surf, surf_res, *atoms->resinfo[i].name))
{
bb[j++] = (-1 != search_str(nb, cb, *atoms->resinfo[i].name));
}
if (i != j)
{
- fprintf(stderr, "Not all residues were recognized (%d from %d), the result may be inaccurate!\n", j, i);
+ fprintf(stderr,
+ "Not all residues were recognized (%d from %d), the result may be inaccurate!\n",
+ j,
+ i);
}
for (i = 0; (i < n_surf); i++)
return bb;
}
-static void check_oo(t_atoms *atoms)
+static void check_oo(t_atoms* atoms)
{
- char *OOO;
-
- int i;
-
- OOO = gmx_strdup("O");
+ char* OOO = gmx_strdup("O");
- for (i = 0; (i < atoms->nr); i++)
+ for (int i = 0; (i < atoms->nr); i++)
{
- if (std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
- {
- *atoms->atomname[i] = OOO;
- }
- else if (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
- {
- *atoms->atomname[i] = OOO;
- }
- else if (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
+ if ((std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
+ || (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
+ || (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
+ || (std::strcmp(*(atoms->atomname[i]), "OT1") == 0))
{
*atoms->atomname[i] = OOO;
}
}
}
-static void norm_acc(t_atoms *atoms, int nres,
- const real av_area[], real norm_av_area[])
+static void norm_acc(t_atoms* atoms, int nres, const real av_area[], real norm_av_area[])
{
- int i, n, n_surf;
+ int i, n, n_surf;
char surffn[] = "surface.dat";
- char **surf_res, **surf_lines;
- double *surf;
+ char ** surf_res, **surf_lines;
+ double* surf;
n_surf = get_lines(surffn, &surf_lines);
snew(surf, n_surf);
}
else
{
- fprintf(stderr, "Residue %s not found in surface database (%s)\n",
- *atoms->resinfo[i].name, surffn);
+ fprintf(stderr, "Residue %s not found in surface database (%s)\n", *atoms->resinfo[i].name, surffn);
}
}
}
-static void prune_ss_legend(t_matrix *mat)
+static void prune_ss_legend(t_matrix* mat)
{
std::vector<bool> isPresent(mat->map.size());
std::vector<int> newnum(mat->map.size());
}
}
-static void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_matrix *mat)
+static void write_sas_mat(const char* fn, real** accr, int nframe, int nres, t_matrix* mat)
{
real lo, hi;
int i, j, nlev;
- t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
- FILE *fp;
+ t_rgb rlo = { 1, 1, 1 }, rhi = { 0, 0, 0 };
+ FILE* fp;
if (fn)
{
}
}
fp = gmx_ffopen(fn, "w");
- nlev = static_cast<int>(hi-lo+1);
- write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)",
- "Time", "Residue Index", nframe, nres,
- mat->axis_x.data(), mat->axis_y.data(), accr, lo, hi, rlo, rhi, &nlev);
+ nlev = static_cast<int>(hi - lo + 1);
+ write_xpm(fp,
+ 0,
+ "Solvent Accessible Surface",
+ "Surface (A^2)",
+ "Time",
+ "Residue Index",
+ nframe,
+ nres,
+ mat->axis_x.data(),
+ mat->axis_y.data(),
+ accr,
+ lo,
+ hi,
+ rlo,
+ rhi,
+ &nlev);
gmx_ffclose(fp);
}
}
-static void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
- const gmx_output_env_t *oenv)
+static void analyse_ss(const char* outfile, t_matrix* mat, const char* ss_string, const gmx_output_env_t* oenv)
{
- FILE *fp;
- int ss_count, total_count;
+ FILE* fp;
+ int ss_count, total_count;
gmx::ArrayRef<t_mapping> map = mat->map;
std::vector<int> count(map.size());
// This copying would not be necessary if xvgr_legend could take a
// view of string views
std::vector<std::string> leg;
- leg.reserve(map.size()+1);
+ leg.reserve(map.size() + 1);
leg.emplace_back("Structure");
- for (const auto &m : map)
+ for (const auto& m : map)
{
leg.emplace_back(m.desc);
}
- fp = xvgropen(outfile, "Secondary Structure",
- output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
+ fp = xvgropen(
+ outfile, "Secondary Structure", output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
if (output_env_get_print_xvgr_codes(oenv))
{
fprintf(fp, "@ subtitle \"Structure = ");
{
fprintf(fp, " + ");
}
- for (const auto &m : map)
+ for (const auto& m : map)
{
if (ss_string[s] == m.code.c1)
{
for (int f = 0; f < mat->nx; f++)
{
ss_count = 0;
- for (auto &c : count)
+ for (auto& c : count)
{
c = 0;
}
{
if (std::strchr(ss_string, map[s].code.c1))
{
- ss_count += count[s];
+ ss_count += count[s];
total_count += count[s];
}
}
fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
- for (const auto &c : count)
+ for (const auto& c : count)
{
fprintf(fp, " %5d", c);
}
}
/* now print column totals */
fprintf(fp, "%-8s %5d", "# Totals", total_count);
- for (const auto &t : total)
+ for (const auto& t : total)
{
fprintf(fp, " %5d", t);
}
/* now print probabilities */
fprintf(fp, "%-8s %5.2f", "# SS pr.", total_count / static_cast<real>(mat->nx * mat->ny));
- for (const auto &t : total)
+ for (const auto& t : total)
{
fprintf(fp, " %5.2f", t / static_cast<real>(mat->nx * mat->ny));
}
xvgrclose(fp);
}
-int gmx_do_dssp(int argc, char *argv[])
+int gmx_do_dssp(int argc, char* argv[])
{
- const char *desc[] = {
+ const char* desc[] = {
"[THISMODULE] ",
"reads a trajectory file and computes the secondary structure for",
"each time frame ",
"calling the dssp program. If you do not have the dssp program,",
- "get it from http://swift.cmbi.ru.nl/gv/dssp. [THISMODULE] assumes ",
+ "get it from https://swift.cmbi.umcn.nl/gv/dssp. [THISMODULE] assumes ",
"that the dssp executable is located in ",
- "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
- "set an environment variable [TT]DSSP[tt] pointing to the dssp",
+ // NOLINTNEXTLINE(bugprone-suspicious-missing-comma)
+ "[TT]" GMX_DSSP_PROGRAM_PATH "[tt]. If this is not the case, then you should",
+ "set an environment variable [TT]DSSP[tt] pointing to the dssp ",
"executable, e.g.: [PAR]",
"[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
- "Since version 2.0.0, dssp is invoked with a syntax that differs",
- "from earlier versions. If you have an older version of dssp,",
- "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
+ "The dssp program is invoked with a syntax that differs",
+ "depending on version. Version 1, 2 and 4 are supported, and the correct",
+ "invocation format can be selected using the [TT]-ver[tt] option.",
"By default, do_dssp uses the syntax introduced with version 2.0.0.",
- "Even newer versions (which at the time of writing are not yet released)",
- "are assumed to have the same syntax as 2.0.0.[PAR]",
+ "Newer versions might also have executable name [TT]mkdssp[tt] instead",
+ "of [TT]dssp[tt].[PAR]",
"The structure assignment for each residue and time is written to an",
"[REF].xpm[ref] matrix file. This file can be visualized with for instance",
"[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
"function of secondary structure type."
};
static gmx_bool bVerbose;
- static const char *ss_string = "HEBT";
+ static const char* ss_string = "HEBT";
static int dsspVersion = 2;
t_pargs pa[] = {
- { "-v", FALSE, etBOOL, {&bVerbose},
- "HIDDENGenerate miles of useless information" },
- { "-sss", FALSE, etSTR, {&ss_string},
- "Secondary structures for structure count"},
- { "-ver", FALSE, etINT, {&dsspVersion},
- "DSSP major version. Syntax changed with version 2"}
+ { "-v", FALSE, etBOOL, { &bVerbose }, "HIDDENGenerate miles of useless information" },
+ { "-sss", FALSE, etSTR, { &ss_string }, "Secondary structures for structure count" },
+ { "-ver",
+ FALSE,
+ etINT,
+ { &dsspVersion },
+ "DSSP major version. Syntax changed with version 2 and 4." }
};
- t_trxstatus *status;
- FILE *tapein, *tapeout;
- FILE *ss, *acc, *fTArea, *tmpf;
- const char *fnSCount, *fnArea, *fnTArea, *fnAArea;
- const char *leg[] = { "Phobic", "Phylic" };
- t_topology top;
- int ePBC;
- t_atoms *atoms;
- t_matrix mat;
- int nres, nr0, naccr, nres_plus_separators;
- gmx_bool *bPhbres, bDoAccSurf;
- real t;
- int natoms, nframe = 0;
- matrix box = {{0}};
- int gnx;
- char *grpnm;
- int *index;
- rvec *xp, *x;
- int *average_area;
- real **accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
- char pdbfile[32], tmpfile[32];
- char dssp[256];
- const char *dptr;
- gmx_output_env_t *oenv;
- gmx_rmpbc_t gpbc = nullptr;
-
- t_filenm fnm[] = {
- { efTRX, "-f", nullptr, ffREAD },
- { efTPS, nullptr, nullptr, ffREAD },
- { efNDX, nullptr, nullptr, ffOPTRD },
- { efDAT, "-ssdump", "ssdump", ffOPTWR },
- { efMAP, "-map", "ss", ffLIBRD },
- { efXPM, "-o", "ss", ffWRITE },
- { efXVG, "-sc", "scount", ffWRITE },
- { efXPM, "-a", "area", ffOPTWR },
- { efXVG, "-ta", "totarea", ffOPTWR },
- { efXVG, "-aa", "averarea", ffOPTWR }
+ t_trxstatus* status;
+ FILE * tapein, *tapeout;
+ FILE * ss, *acc, *fTArea, *tmpf;
+ const char * fnSCount, *fnArea, *fnTArea, *fnAArea;
+ const char* leg[] = { "Phobic", "Phylic" };
+ t_topology top;
+ PbcType pbcType;
+ t_atoms* atoms;
+ t_matrix mat;
+ int nres, nr0, naccr, nres_plus_separators;
+ gmx_bool * bPhbres, bDoAccSurf;
+ real t;
+ int natoms, nframe = 0;
+ matrix box = { { 0 } };
+ int gnx;
+ char* grpnm;
+ int* index;
+ rvec * xp, *x;
+ int* average_area;
+ real ** accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
+ char pdbfile[32], tmpfile[32];
+ char dssp[256];
+ const char* dptr;
+ gmx_output_env_t* oenv;
+ gmx_rmpbc_t gpbc = nullptr;
+
+ t_filenm fnm[] = {
+ { efTRX, "-f", nullptr, ffREAD }, { efTPS, nullptr, nullptr, ffREAD },
+ { efNDX, nullptr, nullptr, ffOPTRD }, { efDAT, "-ssdump", "ssdump", ffOPTWR },
+ { efMAP, "-map", "ss", ffLIBRD }, { efXPM, "-o", "ss", ffWRITE },
+ { efXVG, "-sc", "scount", ffWRITE }, { efXPM, "-a", "area", ffOPTWR },
+ { efXVG, "-ta", "totarea", ffOPTWR }, { efXVG, "-aa", "averarea", ffOPTWR }
};
#define NFILE asize(fnm)
- if (!parse_common_args(&argc, argv,
+ if (!parse_common_args(&argc,
+ argv,
PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
+ NFILE,
+ fnm,
+ asize(pa),
+ pa,
+ asize(desc),
+ desc,
+ 0,
+ nullptr,
+ &oenv))
{
return 0;
}
fnAArea = opt2fn_null("-aa", NFILE, fnm);
bDoAccSurf = ((fnArea != nullptr) || (fnTArea != nullptr) || (fnAArea != nullptr));
- read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xp, nullptr, box, FALSE);
+ read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xp, nullptr, box, FALSE);
atoms = &(top.atoms);
check_oo(atoms);
bPhbres = bPhobics(atoms);
}
fclose(tmpf);
+ const std::string defpathenv = GMX_DSSP_PROGRAM_PATH;
if ((dptr = getenv("DSSP")) == nullptr)
{
- dptr = "/usr/local/bin/dssp";
+ dptr = defpathenv.c_str();
}
if (!gmx_fexist(dptr))
{
- gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)",
- dptr);
+ gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)", dptr);
}
- std::string redirectionString;
+ std::string redirectionString;
redirectionString = prepareToRedirectStdout(bVerbose);
DsspInputStrings dsspStrings;
dsspStrings.dptr = dptr;
dsspStrings.tmpfile = tmpfile;
if (dsspVersion >= 2)
{
- if (dsspVersion > 2)
+ if (dsspVersion == 4)
{
- printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
- }
+ std::string mkdsspCommandLine = dsspStrings.dptr;
+ mkdsspCommandLine += " --output-format dssp ";
+ mkdsspCommandLine += dsspStrings.pdbfile;
- printDsspResult(dssp, dsspStrings, redirectionString);
+#if not HAVE_PIPES && not GMX_NATIVE_WINDOWS
+ // Without pipe/popen, rely on temporary file for output
+ mkdsspCommandLine += " " + dsspStrings.tmpfile;
+#endif
+
+ GMX_RELEASE_ASSERT(mkdsspCommandLine.size() < 255, "DSSP v4 command line too long");
+ strcpy(dssp, mkdsspCommandLine.c_str());
+ }
+ else if (dsspVersion == 2)
+ {
+ printDsspResult(dssp, dsspStrings, redirectionString);
+ }
+ else
+ {
+ printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by "
+ "do_dssp. Assuming version 2 syntax.\n\n",
+ dsspVersion);
+ }
}
else
{
if (fnTArea)
{
- fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",
- output_env_get_xvgr_tlabel(oenv), "Area (nm\\S2\\N)", oenv);
+ fTArea = xvgropen(fnTArea,
+ "Solvent Accessible Surface Area",
+ output_env_get_xvgr_tlabel(oenv),
+ "Area (nm\\S2\\N)",
+ oenv);
xvgr_legend(fTArea, 2, leg, oenv);
}
else
accr = nullptr;
naccr = 0;
- gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
+ gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
do
{
t = output_env_conv_time(oenv, t);
{
naccr += 10;
srenew(accr, naccr);
- for (int i = naccr-10; i < naccr; i++)
+ for (int i = naccr - 10; i < naccr; i++)
{
- snew(accr[i], 2*atoms->nres-1);
+ snew(accr[i], 2 * atoms->nres - 1);
}
}
gmx_rmpbc(gpbc, natoms, box, x);
tapein = gmx_ffopen(pdbfile, "w");
- write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, FALSE);
+ write_pdbfile_indexed(
+ tapein, nullptr, atoms, x, pbcType, box, 'A', -1, gnx, index, nullptr, FALSE, true);
gmx_ffclose(tapein);
/* strip_dssp returns the number of lines found in the dssp file, i.e.
* the number of residues plus the separator lines */
{
remove(pdbfile);
remove(tmpfile);
- gmx_fatal(FARGS, "Failed to execute command: %s\n"
- "Try specifying your dssp version with the -ver option.", dssp);
+ gmx_fatal(FARGS,
+ "Failed to execute command: %s\n"
+ "Try specifying your dssp version with the -ver option.",
+ dssp);
}
if (bDoAccSurf)
{
}
/* strip_dssp returns the number of lines found in the dssp file, i.e.
* the number of residues plus the separator lines */
- nres_plus_separators = strip_dssp(tapeout, nres, bPhbres, t,
- accr_ptr, fTArea, &mat, average_area, oenv);
+ nres_plus_separators =
+ strip_dssp(tapeout, nres, bPhbres, t, accr_ptr, fTArea, &mat, average_area, oenv);
#if HAVE_PIPES || GMX_NATIVE_WINDOWS
pclose(tapeout);
#else
remove(tmpfile);
remove(pdbfile);
nframe++;
- }
- while (read_next_x(oenv, status, &t, x, box));
+ } while (read_next_x(oenv, status, &t, x, box));
fprintf(stderr, "\n");
close_trx(status);
if (fTArea)
if (fnAArea)
{
- acc = xvgropen(fnAArea, "Average Accessible Area",
- "Residue", "A\\S2", oenv);
+ acc = xvgropen(fnAArea, "Average Accessible Area", "Residue", "A\\S2", oenv);
for (int i = 0; (i < nres); i++)
{
- fprintf(acc, "%5d %10g %10g\n", i+1, av_area[i], norm_av_area[i]);
+ fprintf(acc, "%5d %10g %10g\n", i + 1, av_area[i], norm_av_area[i]);
}
xvgrclose(acc);
}