* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-#include <math.h>
+#include "gmxpre.h"
-/*#define HAVE_NN_LOOPS*/
+#include "config.h"
-#include "gromacs/commandline/pargs.h"
-#include "copyrite.h"
-#include "sysstuff.h"
-#include "txtdump.h"
-#include "physics.h"
-#include "macros.h"
-#include "gmx_fatal.h"
-#include "index.h"
-#include "gromacs/utility/smalloc.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "gstat.h"
-#include "gromacs/utility/cstringutil.h"
-#include "pbc.h"
-#include "correl.h"
-#include "gmx_ana.h"
-#include "geminate.h"
+#include <math.h>
-#include "gromacs/fileio/futil.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/correlationfunctions/autocorr.h"
+#include "gromacs/correlationfunctions/crosscorr.h"
+#include "gromacs/correlationfunctions/expfit.h"
+#include "gromacs/correlationfunctions/integrate.h"
#include "gromacs/fileio/matio.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/geminate.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/viewit.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxomp.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/snprintf.h"
+
+#include "geminate.h"
+
+/*#define HAVE_NN_LOOPS*/
typedef short int t_E;
typedef int t_EEst;
{
int i;
- if (ISDON(datable[id]) || !datable)
+ if (!datable || ISDON(datable[id]))
{
if (ddd->dptr[id] == NOTSET) /* New donor */
{
integral += x1;
}
integral *= dt;
- gmx_ffclose(fp);
+ xvgrclose(fp);
printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
printf("Note that the lifetime obtained in this manner is close to useless\n");
printf("Use the -ac option instead and check the Forward lifetime\n");
static void normalizeACF(real *ct, real *gt, int nhb, int len)
{
- real ct_fac, gt_fac;
+ real ct_fac, gt_fac = 0;
int i;
/* Xu and Berne use the same normalization constant */
ct_fac = 1.0/ct[0];
- gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
+ if (nhb != 0)
+ {
+ gt_fac = 1.0/(real)nhb;
+ }
printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
for (i = 0; i < len; i++)
fprintf(fp, "%10g %10g %10g %10g %10g\n",
hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
}
- gmx_ffclose(fp);
+ xvgrclose(fp);
analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
fit_start, temp);
/*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/
}
-static void analyse_donor_props(const char *fn, t_hbdata *hb, int nframes, real t,
- const output_env_t oenv)
+static FILE *open_donor_properties_file(const char *fn,
+ t_hbdata *hb,
+ const output_env_t oenv)
{
- static FILE *fp = NULL;
- const char *leg[] = { "Nbound", "Nfree" };
- int i, j, k, nbound, nb, nhtot;
+ FILE *fp = NULL;
+ const char *leg[] = { "Nbound", "Nfree" };
- if (!fn)
+ if (!fn || !hb)
{
- return;
+ return NULL;
}
- if (!fp)
+
+ fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
+ xvgr_legend(fp, asize(leg), leg, oenv);
+
+ return fp;
+}
+
+static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
+{
+ int i, j, k, nbound, nb, nhtot;
+
+ if (!fp || !hb)
{
- fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
- xvgr_legend(fp, asize(leg), leg, oenv);
+ return;
}
nbound = 0;
nhtot = 0;
for (i = 0; i < isize[grp]; i++)
{
fprintf(fp, (i%15) ? " " : "\n");
- fprintf(fp, " %4u", index[grp][i]+1);
+ fprintf(fp, " %4d", index[grp][i]+1);
}
fprintf(fp, "\n");
/*
{
for (j = 0; (j < hb->d.nhydro[i]); j++)
{
- fprintf(fp, " %4u %4u", hb->d.don[i]+1,
+ fprintf(fp, " %4d %4d", hb->d.don[i]+1,
hb->d.hydro[i][j]+1);
}
fprintf(fp, "\n");
if (hb->a.grp[i] == grp)
{
fprintf(fp, (i%15 && !first) ? " " : "\n");
- fprintf(fp, " %4u", hb->a.acc[i]+1);
+ fprintf(fp, " %4d", hb->a.acc[i]+1);
first = FALSE;
}
}
sprintf(as, "%s", mkatomname(atoms, aaa));
if (bContact)
{
- fprintf(fp, " %6u %6u\n", ddd+1, aaa+1);
+ fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
if (fplog)
{
fprintf(fplog, "%12s %12s\n", ds, as);
{
hhh = hb->d.hydro[i][m];
sprintf(hs, "%s", mkatomname(atoms, hhh));
- fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1);
+ fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
if (fplog)
{
fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
/* "It is also possible to analyse specific hydrogen bonds with",
"[TT]-sel[tt]. This index file must contain a group of atom triplets",
- "Donor Hydrogen Acceptor, in the following way:[PAR]",
+ "Donor Hydrogen Acceptor, in the following way::",
+ "",
+ "[ selected ]",
+ " 20 21 24",
+ " 25 26 29",
+ " 1 3 6",
+ "",
+ "Note that the triplets need not be on separate lines.",
+ "Each atom triplet specifies a hydrogen bond to be analyzed,",
+ "note also that no check is made for the types of atoms.[PAR]",
*/
- "[TT]",
- "[ selected ][BR]",
- " 20 21 24[BR]",
- " 25 26 29[BR]",
- " 1 3 6[BR]",
- "[tt][BR]",
- "Note that the triplets need not be on separate lines.",
- "Each atom triplet specifies a hydrogen bond to be analyzed,",
- "note also that no check is made for the types of atoms.[PAR]",
-
- "[BB]Output:[bb][BR]",
- "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]",
- "[TT]-ac[tt]: average over all autocorrelations of the existence",
- "functions (either 0 or 1) of all hydrogen bonds.[BR]",
- "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]",
- "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]",
- "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
- "where n and n+i stand for residue numbers and i ranges from 0 to 6.",
- "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
- "with helices in proteins.[BR]",
- "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
- "for selected groups, all hydrogen bonded atoms from all groups and",
- "all solvent atoms involved in insertion.[BR]",
- "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
- "frames, this also contains information on solvent insertion",
- "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
- "index file.[BR]",
- "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
- "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]",
- "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
- "compare results to Raman Spectroscopy.",
- "[PAR]",
+
+ "[BB]Output:[bb]",
+ "",
+ " * [TT]-num[tt]: number of hydrogen bonds as a function of time.",
+ " * [TT]-ac[tt]: average over all autocorrelations of the existence",
+ " functions (either 0 or 1) of all hydrogen bonds.",
+ " * [TT]-dist[tt]: distance distribution of all hydrogen bonds.",
+ " * [TT]-ang[tt]: angle distribution of all hydrogen bonds.",
+ " * [TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time",
+ " where n and n+i stand for residue numbers and i ranges from 0 to 6.",
+ " This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated",
+ " with helices in proteins.",
+ " * [TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors",
+ " for selected groups, all hydrogen bonded atoms from all groups and",
+ " all solvent atoms involved in insertion.",
+ " * [TT]-hbm[tt]: existence matrix for all hydrogen bonds over all",
+ " frames, this also contains information on solvent insertion",
+ " into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]",
+ " index file.",
+ " * [TT]-dan[tt]: write out the number of donors and acceptors analyzed for",
+ " each timeframe. This is especially useful when using [TT]-shell[tt].",
+ " * [TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to",
+ " compare results to Raman Spectroscopy.",
+ "",
"Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]",
"require an amount of memory proportional to the total numbers of donors",
"times the total number of acceptors in the selected group(s)."
{ "-temp", FALSE, etREAL, {&temp},
"Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
{ "-dump", FALSE, etINT, {&nDump},
- "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" },
+ "Dump the first N hydrogen bond ACFs in a single [REF].xvg[ref] file for debugging" },
{ "-max_hb", FALSE, etREAL, {&maxnhb},
"Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" },
{ "-merge", FALSE, etBOOL, {&bMerge},
};
t_filenm fnm[] = {
{ efTRX, "-f", NULL, ffREAD },
- { efTPX, NULL, NULL, ffREAD },
+ { efTPR, NULL, NULL, ffREAD },
{ efNDX, NULL, NULL, ffOPTRD },
/* { efNDX, "-sel", "select", ffOPTRD },*/
{ efXVG, "-num", "hbnum", ffWRITE },
int grp, nabin, nrbin, bin, resdist, ihb;
char **leg;
t_hbdata *hb, *hbptr;
- FILE *fp, *fpins = NULL, *fpnhb = NULL;
+ FILE *fp, *fpins = NULL, *fpnhb = NULL, *donor_properties = NULL;
t_gridcell ***grid;
t_ncell *icell, *jcell, *kcell;
ivec ngrid;
npargs = asize(pa);
ppa = add_acf_pargs(&npargs, pa);
- if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE, NFILE, fnm, npargs,
+ if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
{
return 0;
hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
/* get topology */
- read_tpx_top(ftp2fn(efTPX, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
+ read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
snew(grpnames, grNR);
snew(index, grNR);
/*if (bSelected)
snew(donors[gr0D], dons[gr0D].nrd);*/
+ donor_properties = open_donor_properties_file(opt2fn_null("-don", NFILE, fnm), hb, oenv);
+
if (bHBmap)
{
printf("Making hbmap structure...");
#pragma omp barrier
#pragma omp single
{
- if (hb != NULL)
- {
- analyse_donor_props(opt2fn_null("-don", NFILE, fnm), hb, k, t, oenv);
- }
+ analyse_donor_properties(donor_properties, hb, k, t);
}
#pragma omp single
free_grid(ngrid, &grid);
close_trj(status);
+
+ if (donor_properties)
+ {
+ xvgrclose(donor_properties);
+ }
+
if (fpnhb)
{
- gmx_ffclose(fpnhb);
+ xvgrclose(fpnhb);
}
/* Compute maximum possible number of different hbonds */
aver_nhb += hb->nhb[i];
aver_dist += hb->ndist[i];
}
- gmx_ffclose(fp);
+ xvgrclose(fp);
aver_nhb /= nframes;
aver_dist /= nframes;
/* Print HB distance distribution */
{
fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
}
- gmx_ffclose(fp);
+ xvgrclose(fp);
}
/* Print HB angle distribution */
{
fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
}
- gmx_ffclose(fp);
+ xvgrclose(fp);
}
/* Print HB in alpha-helix */
}
fprintf(fp, "\n");
}
- gmx_ffclose(fp);
+ xvgrclose(fp);
}
if (!bNN)
{
gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
}
}
- gemstring = strdup(gemType[hb->per->gemtype]);
+ gemstring = gmx_strdup(gemType[hb->per->gemtype]);
do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
bMerge, bContact, fit_start, temp, r2cut > 0, oenv,
gemstring, nThreads, NN, bBallistic, bGemFit);
{
t_matrix mat;
int id, ia, hh, x, y;
+ mat.flags = mat.y0 = 0;
if ((nframes > 0) && (hb->nrhb > 0))
{
if (USE_THIS_GROUP(j) )
{
sprintf(buf, "Donors %s", grpnames[j]);
- legnames[i++] = strdup(buf);
+ legnames[i++] = gmx_strdup(buf);
sprintf(buf, "Acceptors %s", grpnames[j]);
- legnames[i++] = strdup(buf);
+ legnames[i++] = gmx_strdup(buf);
}
}
if (i != nleg)
}
fprintf(fp, "\n");
}
- gmx_ffclose(fp);
+ xvgrclose(fp);
}
return 0;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff