*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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/smalloc.h"
#include "gromacs/utility/sysinfo.h"
-int gmx_covar(int argc, char *argv[])
+int gmx_covar(int argc, char* argv[])
{
- const char *desc[] = {
+ const char* desc[] = {
"[THISMODULE] calculates and diagonalizes the (mass-weighted)",
"covariance matrix.",
"All structures are fitted to the structure in the structure file.",
"should consider carefully whether a reduced set of atoms will meet",
"your needs for lower costs."
};
- static gmx_bool bFit = TRUE, bRef = FALSE, bM = FALSE, bPBC = TRUE;
- static int end = -1;
- t_pargs pa[] = {
- { "-fit", FALSE, etBOOL, {&bFit},
- "Fit to a reference structure"},
- { "-ref", FALSE, etBOOL, {&bRef},
- "Use the deviation from the conformation in the structure file instead of from the average" },
- { "-mwa", FALSE, etBOOL, {&bM},
- "Mass-weighted covariance analysis"},
- { "-last", FALSE, etINT, {&end},
- "Last eigenvector to write away (-1 is till the last)" },
- { "-pbc", FALSE, etBOOL, {&bPBC},
- "Apply corrections for periodic boundary conditions" }
+ static gmx_bool bFit = TRUE, bRef = FALSE, bM = FALSE, bPBC = TRUE;
+ static int end = -1;
+ t_pargs pa[] = {
+ { "-fit", FALSE, etBOOL, { &bFit }, "Fit to a reference structure" },
+ { "-ref",
+ FALSE,
+ etBOOL,
+ { &bRef },
+ "Use the deviation from the conformation in the structure file instead of from the "
+ "average" },
+ { "-mwa", FALSE, etBOOL, { &bM }, "Mass-weighted covariance analysis" },
+ { "-last", FALSE, etINT, { &end }, "Last eigenvector to write away (-1 is till the last)" },
+ { "-pbc", FALSE, etBOOL, { &bPBC }, "Apply corrections for periodic boundary conditions" }
};
- FILE *out = nullptr; /* initialization makes all compilers happy */
- t_trxstatus *status;
+ FILE* out = nullptr; /* initialization makes all compilers happy */
+ t_trxstatus* status;
t_topology top;
int ePBC;
- t_atoms *atoms;
- rvec *x, *xread, *xref, *xav, *xproj;
+ t_atoms* atoms;
+ rvec * x, *xread, *xref, *xav, *xproj;
matrix box, zerobox;
- real *sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
+ real * sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
real t, tstart, tend, **mat2;
real xj, *w_rls = nullptr;
real min, max, *axis;
int natoms, nat, nframes0, nframes, nlevels;
int64_t ndim, i, j, k, l;
int WriteXref;
- const char *fitfile, *trxfile, *ndxfile;
- const char *eigvalfile, *eigvecfile, *averfile, *logfile;
- const char *asciifile, *xpmfile, *xpmafile;
+ const char * fitfile, *trxfile, *ndxfile;
+ const char * eigvalfile, *eigvecfile, *averfile, *logfile;
+ const char * asciifile, *xpmfile, *xpmafile;
char str[STRLEN], *fitname, *ananame;
int d, dj, nfit;
- int *index, *ifit;
+ int * index, *ifit;
gmx_bool bDiffMass1, bDiffMass2;
t_rgb rlo, rmi, rhi;
- real *eigenvectors;
- gmx_output_env_t *oenv;
+ real* eigenvectors;
+ 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 },
- { efXVG, nullptr, "eigenval", ffWRITE },
- { efTRN, "-v", "eigenvec", ffWRITE },
- { efSTO, "-av", "average.pdb", ffWRITE },
- { efLOG, nullptr, "covar", ffWRITE },
- { efDAT, "-ascii", "covar", ffOPTWR },
- { efXPM, "-xpm", "covar", ffOPTWR },
- { efXPM, "-xpma", "covara", ffOPTWR }
+ t_filenm fnm[] = {
+ { efTRX, "-f", nullptr, ffREAD }, { efTPS, nullptr, nullptr, ffREAD },
+ { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, nullptr, "eigenval", ffWRITE },
+ { efTRN, "-v", "eigenvec", ffWRITE }, { efSTO, "-av", "average.pdb", ffWRITE },
+ { efLOG, nullptr, "covar", ffWRITE }, { efDAT, "-ascii", "covar", ffOPTWR },
+ { efXPM, "-xpm", "covar", ffOPTWR }, { efXPM, "-xpma", "covara", ffOPTWR }
};
#define NFILE asize(fnm)
- if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
+ if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, asize(pa), pa,
+ asize(desc), desc, 0, nullptr, &oenv))
{
return 0;
}
w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
if (i)
{
- bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i-1]]);
+ bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i - 1]]);
}
}
}
sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
if (i)
{
- bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i-1]);
+ bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i - 1]);
}
}
else
}
if (!bDiffMass1)
{
- fprintf(stderr, "\n"
+ fprintf(stderr,
+ "\n"
"Note: the fit and analysis group are identical,\n"
" while the fit is mass weighted and the analysis is not.\n"
" Making the fit non mass weighted.\n\n");
snew(x, natoms);
snew(xav, natoms);
- ndim = natoms*DIM;
+ ndim = natoms * DIM;
if (std::sqrt(static_cast<real>(INT64_MAX)) < static_cast<real>(ndim))
{
gmx_fatal(FARGS, "Number of degrees of freedoms to large for matrix.\n");
}
- snew(mat, ndim*ndim);
+ snew(mat, ndim * ndim);
fprintf(stderr, "Calculating the average structure ...\n");
nframes0 = 0;
nat = read_first_x(oenv, &status, trxfile, &t, &xread, box);
if (nat != atoms->nr)
{
- fprintf(stderr, "\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n", natoms, nat);
+ fprintf(stderr, "\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",
+ natoms, nat);
}
do
{
{
rvec_inc(xav[i], xread[index[i]]);
}
- }
- while (read_next_x(oenv, status, &t, xread, box));
+ } while (read_next_x(oenv, status, &t, xread, box));
close_trx(status);
- inv_nframes = 1.0/nframes0;
+ inv_nframes = 1.0 / nframes0;
for (i = 0; i < natoms; i++)
{
for (d = 0; d < DIM; d++)
{
- xav[i][d] *= inv_nframes;
+ xav[i][d] *= inv_nframes;
xread[index[i]][d] = xav[i][d];
}
}
- write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure",
- atoms, xread, nullptr, epbcNONE, zerobox, natoms, index);
+ write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure", atoms, xread, nullptr,
+ epbcNONE, zerobox, natoms, index);
sfree(xread);
- fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", static_cast<int>(ndim), static_cast<int>(ndim));
+ fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", static_cast<int>(ndim),
+ static_cast<int>(ndim));
nframes = 0;
nat = read_first_x(oenv, &status, trxfile, &t, &xread, box);
tstart = t;
{
for (dj = 0; dj < DIM; dj++)
{
- k = ndim*(DIM*j+dj);
+ k = ndim * (DIM * j + dj);
xj = x[j][dj];
for (i = j; i < natoms; i++)
{
- l = k+DIM*i;
+ l = k + DIM * i;
for (d = 0; d < DIM; d++)
{
- mat[l+d] += x[i][d]*xj;
+ mat[l + d] += x[i][d] * xj;
}
}
}
}
- }
- while (read_next_x(oenv, status, &t, xread, box) &&
- (bRef || nframes < nframes0));
+ } while (read_next_x(oenv, status, &t, xread, box) && (bRef || nframes < nframes0));
close_trx(status);
gmx_rmpbc_done(gpbc);
}
/* correct the covariance matrix for the mass */
- inv_nframes = 1.0/nframes;
+ inv_nframes = 1.0 / nframes;
for (j = 0; j < natoms; j++)
{
for (dj = 0; dj < DIM; dj++)
{
for (i = j; i < natoms; i++)
{
- k = ndim*(DIM*j+dj)+DIM*i;
+ k = ndim * (DIM * j + dj) + DIM * i;
for (d = 0; d < DIM; d++)
{
- mat[k+d] = mat[k+d]*inv_nframes*sqrtm[i]*sqrtm[j];
+ mat[k + d] = mat[k + d] * inv_nframes * sqrtm[i] * sqrtm[j];
}
}
}
{
for (i = j; i < ndim; i++)
{
- mat[ndim*i+j] = mat[ndim*j+i];
+ mat[ndim * i + j] = mat[ndim * j + i];
}
}
trace = 0;
for (i = 0; i < ndim; i++)
{
- trace += mat[i*ndim+i];
+ trace += mat[i * ndim + i];
}
- fprintf(stderr, "\nTrace of the covariance matrix: %g (%snm^2)\n",
- trace, bM ? "u " : "");
+ fprintf(stderr, "\nTrace of the covariance matrix: %g (%snm^2)\n", trace, bM ? "u " : "");
if (asciifile)
{
{
for (i = 0; i < ndim; i += 3)
{
- fprintf(out, "%g %g %g\n",
- mat[ndim*j+i], mat[ndim*j+i+1], mat[ndim*j+i+2]);
+ fprintf(out, "%g %g %g\n", mat[ndim * j + i], mat[ndim * j + i + 1],
+ mat[ndim * j + i + 2]);
}
}
gmx_ffclose(out);
snew(mat2, ndim);
for (j = 0; j < ndim; j++)
{
- mat2[j] = &(mat[ndim*j]);
+ mat2[j] = &(mat[ndim * j]);
for (i = 0; i <= j; i++)
{
if (mat2[j][i] < min)
snew(axis, ndim);
for (i = 0; i < ndim; i++)
{
- axis[i] = i+1;
+ axis[i] = i + 1;
}
- rlo.r = 0; rlo.g = 0; rlo.b = 1;
- rmi.r = 1; rmi.g = 1; rmi.b = 1;
- rhi.r = 1; rhi.g = 0; rhi.b = 0;
+ rlo.r = 0;
+ rlo.g = 0;
+ rlo.b = 1;
+ rmi.r = 1;
+ rmi.g = 1;
+ rmi.b = 1;
+ rhi.r = 1;
+ rhi.g = 0;
+ rhi.b = 0;
out = gmx_ffopen(xpmfile, "w");
nlevels = 80;
- write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
- "dim", "dim", ndim, ndim, axis, axis,
- mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
+ write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2", "dim", "dim", ndim, ndim, axis,
+ axis, mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
gmx_ffclose(out);
sfree(axis);
sfree(mat2);
{
min = 0;
max = 0;
- snew(mat2, ndim/DIM);
- for (i = 0; i < ndim/DIM; i++)
+ snew(mat2, ndim / DIM);
+ for (i = 0; i < ndim / DIM; i++)
{
- snew(mat2[i], ndim/DIM);
+ snew(mat2[i], ndim / DIM);
}
- for (j = 0; j < ndim/DIM; j++)
+ for (j = 0; j < ndim / DIM; j++)
{
for (i = 0; i <= j; i++)
{
mat2[j][i] = 0;
for (d = 0; d < DIM; d++)
{
- mat2[j][i] += mat[ndim*(DIM*j+d)+DIM*i+d];
+ mat2[j][i] += mat[ndim * (DIM * j + d) + DIM * i + d];
}
if (mat2[j][i] < min)
{
mat2[i][j] = mat2[j][i];
}
}
- snew(axis, ndim/DIM);
- for (i = 0; i < ndim/DIM; i++)
+ snew(axis, ndim / DIM);
+ for (i = 0; i < ndim / DIM; i++)
{
- axis[i] = i+1;
+ axis[i] = i + 1;
}
- rlo.r = 0; rlo.g = 0; rlo.b = 1;
- rmi.r = 1; rmi.g = 1; rmi.b = 1;
- rhi.r = 1; rhi.g = 0; rhi.b = 0;
+ rlo.r = 0;
+ rlo.g = 0;
+ rlo.b = 1;
+ rmi.r = 1;
+ rmi.g = 1;
+ rmi.b = 1;
+ rhi.r = 1;
+ rhi.g = 0;
+ rhi.b = 0;
out = gmx_ffopen(xpmafile, "w");
nlevels = 80;
- write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
- "atom", "atom", ndim/DIM, ndim/DIM, axis, axis,
- mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
+ write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2", "atom", "atom", ndim / DIM,
+ ndim / DIM, axis, axis, mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
gmx_ffclose(out);
sfree(axis);
- for (i = 0; i < ndim/DIM; i++)
+ for (i = 0; i < ndim / DIM; i++)
{
sfree(mat2[i]);
}
/* call diagonalization routine */
snew(eigenvalues, ndim);
- snew(eigenvectors, ndim*ndim);
+ snew(eigenvectors, ndim * ndim);
- std::memcpy(eigenvectors, mat, ndim*ndim*sizeof(real));
+ std::memcpy(eigenvectors, mat, ndim * ndim * sizeof(real));
fprintf(stderr, "\nDiagonalizing ...\n");
fflush(stderr);
eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
{
sum += eigenvalues[i];
}
- fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n",
- sum, bM ? "u " : "");
- if (std::abs(trace-sum) > 0.01*trace)
+ fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n", sum, bM ? "u " : "");
+ if (std::abs(trace - sum) > 0.01 * trace)
{
- fprintf(stderr, "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
+ fprintf(stderr,
+ "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
}
/* Set 'end', the maximum eigenvector and -value index used for output */
if (end == -1)
{
- if (nframes-1 < ndim)
+ if (nframes - 1 < ndim)
{
- end = nframes-1;
- fprintf(stderr, "\nWARNING: there are fewer frames in your trajectory than there are\n");
+ end = nframes - 1;
+ fprintf(stderr,
+ "\nWARNING: there are fewer frames in your trajectory than there are\n");
fprintf(stderr, "degrees of freedom in your system. Only generating the first\n");
fprintf(stderr, "%d out of %d eigenvectors and eigenvalues.\n", end, static_cast<int>(ndim));
}
fprintf(stderr, "\nWriting eigenvalues to %s\n", eigvalfile);
sprintf(str, "(%snm\\S2\\N)", bM ? "u " : "");
- out = xvgropen(eigvalfile,
- "Eigenvalues of the covariance matrix",
- "Eigenvector index", str, oenv);
+ out = xvgropen(eigvalfile, "Eigenvalues of the covariance matrix", "Eigenvector index", str, oenv);
for (i = 0; (i < end); i++)
{
- fprintf (out, "%10d %g\n", static_cast<int>(i+1), eigenvalues[ndim-1-i]);
+ fprintf(out, "%10d %g\n", static_cast<int>(i + 1), eigenvalues[ndim - 1 - i]);
}
xvgrclose(out);
WriteXref = eWXR_NOFIT;
}
- write_eigenvectors(eigvecfile, natoms, mat, TRUE, 1, end,
- WriteXref, x, bDiffMass1, xproj, bM, eigenvalues);
+ write_eigenvectors(eigvecfile, natoms, mat, TRUE, 1, end, WriteXref, x, bDiffMass1, xproj, bM,
+ eigenvalues);
out = gmx_ffopen(logfile, "w");
fprintf(out, "Working directory: %s\n\n", str);
fprintf(out, "Read %d frames from %s (time %g to %g %s)\n", nframes, trxfile,
- output_env_conv_time(oenv, tstart), output_env_conv_time(oenv, tend), output_env_get_time_unit(oenv).c_str());
+ output_env_conv_time(oenv, tstart), output_env_conv_time(oenv, tend),
+ output_env_get_time_unit(oenv).c_str());
if (bFit)
{
fprintf(out, "Read reference structure for fit from %s\n", fitfile);
{
fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
}
- fprintf(out, "Diagonalized the %dx%d covariance matrix\n", static_cast<int>(ndim), static_cast<int>(ndim));
- fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n",
- trace);
- fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n",
- sum);
+ fprintf(out, "Diagonalized the %dx%d covariance matrix\n", static_cast<int>(ndim),
+ static_cast<int>(ndim));
+ fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n", trace);
+ fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n", sum);
fprintf(out, "Wrote %d eigenvalues to %s\n", static_cast<int>(end), eigvalfile);
if (WriteXref == eWXR_YES)
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff