#define max_hx 7
typedef int t_hx[max_hx];
#define NRHXTYPES max_hx
-static const char *hxtypenames[NRHXTYPES] =
-{"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"};
+static const char* hxtypenames[NRHXTYPES] = { "n-n", "n-n+1", "n-n+2", "n-n+3",
+ "n-n+4", "n-n+5", "n-n>6" };
#define MAXHH 4
static const int NOTSET = -49297;
/* -----------------------------------------*/
-enum {
- gr0, gr1, grI, grNR
+enum
+{
+ gr0,
+ gr1,
+ grI,
+ grNR
};
-enum {
- hbNo, hbDist, hbHB, hbNR, hbR2
+enum
+{
+ hbNo,
+ hbDist,
+ hbHB,
+ hbNR,
+ hbR2
};
static const unsigned char c_acceptorMask = (1 << 0);
static const unsigned char c_donorMask = (1 << 1);
static const unsigned char c_inGroupMask = (1 << 2);
-static const char *grpnames[grNR] = {"0", "1", "I" };
+static const char* grpnames[grNR] = { "0", "1", "I" };
-static gmx_bool bDebug = FALSE;
+static gmx_bool bDebug = FALSE;
#define HB_NO 0
-#define HB_YES 1<<0
-#define HB_INS 1<<1
-#define HB_YESINS (HB_YES|HB_INS)
-#define HB_NR (1<<2)
+#define HB_YES 1 << 0
+#define HB_INS 1 << 1
+#define HB_YESINS (HB_YES | HB_INS)
+#define HB_NR (1 << 2)
#define MAXHYDRO 4
-#define ISHB(h) ((h) & 2)
-#define ISDIST(h) ((h) & 1)
-#define ISDIST2(h) ((h) & 4)
-#define ISACC(h) ((h) & c_acceptorMask)
-#define ISDON(h) ((h) & c_donorMask)
-#define ISINGRP(h) ((h) & c_inGroupMask)
-
-typedef struct {
- int nr;
- int maxnr;
- int *atoms;
+#define ISHB(h) ((h)&2)
+#define ISDIST(h) ((h)&1)
+#define ISDIST2(h) ((h)&4)
+#define ISACC(h) ((h)&c_acceptorMask)
+#define ISDON(h) ((h)&c_donorMask)
+#define ISINGRP(h) ((h)&c_inGroupMask)
+
+typedef struct
+{
+ int nr;
+ int maxnr;
+ int* atoms;
} t_ncell;
-typedef struct {
+typedef struct
+{
t_ncell d[grNR];
t_ncell a[grNR];
} t_gridcell;
-typedef int t_icell[grNR];
+typedef int t_icell[grNR];
typedef int h_id[MAXHYDRO];
-typedef struct {
- int history[MAXHYDRO];
+typedef struct
+{
+ int history[MAXHYDRO];
/* Has this hbond existed ever? If so as hbDist or hbHB or both.
* Result is stored as a bitmap (1 = hbDist) || (2 = hbHB)
*/
*/
int n0; /* First frame a HB was found */
int nframes, maxframes; /* Amount of frames in this hbond */
- unsigned int **h;
- unsigned int **g;
+ unsigned int** h;
+ unsigned int** g;
/* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the
* function g(t) = [1-h(t)] H(t) where H(t) is one when the donor-
* acceptor distance is less than the user-specified distance (typically
*/
} t_hbond;
-typedef struct {
- int nra, max_nra;
- int *acc; /* Atom numbers of the acceptors */
- int *grp; /* Group index */
- int *aptr; /* Map atom number to acceptor index */
+typedef struct
+{
+ int nra, max_nra;
+ int* acc; /* Atom numbers of the acceptors */
+ int* grp; /* Group index */
+ int* aptr; /* Map atom number to acceptor index */
} t_acceptors;
-typedef struct {
- int nrd, max_nrd;
- int *don; /* Atom numbers of the donors */
- int *grp; /* Group index */
- int *dptr; /* Map atom number to donor index */
- int *nhydro; /* Number of hydrogens for each donor */
- h_id *hydro; /* The atom numbers of the hydrogens */
- h_id *nhbonds; /* The number of HBs per H at current */
+typedef struct
+{
+ int nrd, max_nrd;
+ int* don; /* Atom numbers of the donors */
+ int* grp; /* Group index */
+ int* dptr; /* Map atom number to donor index */
+ int* nhydro; /* Number of hydrogens for each donor */
+ h_id* hydro; /* The atom numbers of the hydrogens */
+ h_id* nhbonds; /* The number of HBs per H at current */
} t_donors;
-typedef struct {
- gmx_bool bHBmap, bDAnr;
- int wordlen;
+typedef struct
+{
+ gmx_bool bHBmap, bDAnr;
+ int wordlen;
/* The following arrays are nframes long */
- int nframes, max_frames, maxhydro;
- int *nhb, *ndist;
- h_id *n_bound;
- real *time;
- t_icell *danr;
- t_hx *nhx;
+ int nframes, max_frames, maxhydro;
+ int * nhb, *ndist;
+ h_id* n_bound;
+ real* time;
+ t_icell* danr;
+ t_hx* nhx;
/* These structures are initialized from the topology at start up */
- t_donors d;
- t_acceptors a;
+ t_donors d;
+ t_acceptors a;
/* This holds a matrix with all possible hydrogen bonds */
- int nrhb, nrdist;
- t_hbond ***hbmap;
+ int nrhb, nrdist;
+ t_hbond*** hbmap;
} t_hbdata;
/* Changed argument 'bMerge' into 'oneHB' below,
* - Erik Marklund May 29, 2006
*/
-static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB)
+static t_hbdata* mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB)
{
- t_hbdata *hb;
+ t_hbdata* hb;
snew(hb, 1);
- hb->wordlen = 8*sizeof(unsigned int);
+ hb->wordlen = 8 * sizeof(unsigned int);
hb->bHBmap = bHBmap;
hb->bDAnr = bDAnr;
if (oneHB)
return hb;
}
-static void mk_hbmap(t_hbdata *hb)
+static void mk_hbmap(t_hbdata* hb)
{
- int i, j;
+ int i, j;
snew(hb->hbmap, hb->d.nrd);
for (i = 0; (i < hb->d.nrd); i++)
}
}
-static void add_frames(t_hbdata *hb, int nframes)
+static void add_frames(t_hbdata* hb, int nframes)
{
if (nframes >= hb->max_frames)
{
}
#define OFFSET(frame) ((frame) / 32)
-#define MASK(frame) (1 << ((frame) % 32))
+#define MASK(frame) (1 << ((frame) % 32))
static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue)
{
return (hbexist[OFFSET(frame)] & MASK(frame)) != 0;
}
-static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb)
+static void set_hb(t_hbdata* hb, int id, int ih, int ia, int frame, int ihb)
{
- unsigned int *ghptr = nullptr;
+ unsigned int* ghptr = nullptr;
if (ihb == hbHB)
{
gmx_fatal(FARGS, "Incomprehensible iValue %d in set_hb", ihb);
}
- _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE);
+ _set_hb(ghptr, frame - hb->hbmap[id][ia]->n0, TRUE);
}
-static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb)
+static void add_ff(t_hbdata* hbd, int id, int h, int ia, int frame, int ihb)
{
- int i, j, n;
- t_hbond *hb = hbd->hbmap[id][ia];
- int maxhydro = std::min(hbd->maxhydro, hbd->d.nhydro[id]);
- int wlen = hbd->wordlen;
- int delta = 32*wlen;
+ int i, j, n;
+ t_hbond* hb = hbd->hbmap[id][ia];
+ int maxhydro = std::min(hbd->maxhydro, hbd->d.nhydro[id]);
+ int wlen = hbd->wordlen;
+ int delta = 32 * wlen;
if (!hb->h[0])
{
hb->maxframes = delta;
for (i = 0; (i < maxhydro); i++)
{
- snew(hb->h[i], hb->maxframes/wlen);
- snew(hb->g[i], hb->maxframes/wlen);
+ snew(hb->h[i], hb->maxframes / wlen);
+ snew(hb->g[i], hb->maxframes / wlen);
}
}
else
{
- hb->nframes = frame-hb->n0;
+ hb->nframes = frame - hb->n0;
/* We need a while loop here because hbonds may be returning
* after a long time.
*/
n = hb->maxframes + delta;
for (i = 0; (i < maxhydro); i++)
{
- srenew(hb->h[i], n/wlen);
- srenew(hb->g[i], n/wlen);
- for (j = hb->maxframes/wlen; (j < n/wlen); j++)
+ srenew(hb->h[i], n / wlen);
+ srenew(hb->g[i], n / wlen);
+ for (j = hb->maxframes / wlen; (j < n / wlen); j++)
{
hb->h[i][j] = 0;
hb->g[i][j] = 0;
{
set_hb(hbd, id, h, ia, frame, ihb);
}
-
}
-static void inc_nhbonds(t_donors *ddd, int d, int h)
+static void inc_nhbonds(t_donors* ddd, int d, int h)
{
int j;
int dptr = ddd->dptr[d];
}
if (j == ddd->nhydro[dptr])
{
- gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h+1, d+1);
+ gmx_fatal(FARGS, "No such hydrogen %d on donor %d\n", h + 1, d + 1);
}
}
-static int _acceptor_index(t_acceptors *a, int grp, int i,
- const char *file, int line)
+static int _acceptor_index(t_acceptors* a, int grp, int i, const char* file, int line)
{
int ai = a->aptr[i];
{
if (debug && bDebug)
{
- fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
- ai, a->grp[ai], grp, file, line);
+ fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n", ai,
+ a->grp[ai], grp, file, line);
}
return NOTSET;
}
}
#define acceptor_index(a, grp, i) _acceptor_index(a, grp, i, __FILE__, __LINE__)
-static int _donor_index(t_donors *d, int grp, int i, const char *file, int line)
+static int _donor_index(t_donors* d, int grp, int i, const char* file, int line)
{
int di = d->dptr[i];
{
if (debug && bDebug)
{
- fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n",
- di, d->grp[di], grp, file, line);
+ fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n", di,
+ d->grp[di], grp, file, line);
}
return NOTSET;
}
}
#define donor_index(d, grp, i) _donor_index(d, grp, i, __FILE__, __LINE__)
-static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd)
+static gmx_bool isInterchangable(t_hbdata* hb, int d, int a, int grpa, int grpd)
{
/* g_hbond doesn't allow overlapping groups */
if (grpa != grpd)
{
return FALSE;
}
- return
- donor_index(&hb->d, grpd, a) != NOTSET
- && acceptor_index(&hb->a, grpa, d) != NOTSET;
+ return donor_index(&hb->d, grpd, a) != NOTSET && acceptor_index(&hb->a, grpa, d) != NOTSET;
}
-static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa,
- int frame, gmx_bool bMerge, int ihb, gmx_bool bContact)
+static void
+add_hbond(t_hbdata* hb, int d, int a, int h, int grpd, int grpa, int frame, gmx_bool bMerge, int ihb, gmx_bool bContact)
{
int k, id, ia, hh;
gmx_bool daSwap = FALSE;
if ((id = hb->d.dptr[d]) == NOTSET)
{
- gmx_fatal(FARGS, "No donor atom %d", d+1);
+ gmx_fatal(FARGS, "No donor atom %d", d + 1);
}
else if (grpd != hb->d.grp[id])
{
- gmx_fatal(FARGS, "Inconsistent donor groups, %d instead of %d, atom %d",
- grpd, hb->d.grp[id], d+1);
+ gmx_fatal(FARGS, "Inconsistent donor groups, %d instead of %d, atom %d", grpd,
+ hb->d.grp[id], d + 1);
}
if ((ia = hb->a.aptr[a]) == NOTSET)
{
- gmx_fatal(FARGS, "No acceptor atom %d", a+1);
+ gmx_fatal(FARGS, "No acceptor atom %d", a + 1);
}
else if (grpa != hb->a.grp[ia])
{
- gmx_fatal(FARGS, "Inconsistent acceptor groups, %d instead of %d, atom %d",
- grpa, hb->a.grp[ia], a+1);
+ gmx_fatal(FARGS, "Inconsistent acceptor groups, %d instead of %d, atom %d", grpa,
+ hb->a.grp[ia], a + 1);
}
if (bMerge)
/* Now repeat donor/acc check. */
if ((id = hb->d.dptr[d]) == NOTSET)
{
- gmx_fatal(FARGS, "No donor atom %d", d+1);
+ gmx_fatal(FARGS, "No donor atom %d", d + 1);
}
else if (grpd != hb->d.grp[id])
{
- gmx_fatal(FARGS, "Inconsistent donor groups, %d instead of %d, atom %d",
- grpd, hb->d.grp[id], d+1);
+ gmx_fatal(FARGS, "Inconsistent donor groups, %d instead of %d, atom %d", grpd,
+ hb->d.grp[id], d + 1);
}
if ((ia = hb->a.aptr[a]) == NOTSET)
{
- gmx_fatal(FARGS, "No acceptor atom %d", a+1);
+ gmx_fatal(FARGS, "No acceptor atom %d", a + 1);
}
else if (grpa != hb->a.grp[ia])
{
- gmx_fatal(FARGS, "Inconsistent acceptor groups, %d instead of %d, atom %d",
- grpa, hb->a.grp[ia], a+1);
+ gmx_fatal(FARGS, "Inconsistent acceptor groups, %d instead of %d, atom %d", grpa,
+ hb->a.grp[ia], a + 1);
}
}
}
}
if (k == hb->d.nhydro[id])
{
- gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)",
- d+1, h+1, a+1);
+ gmx_fatal(FARGS, "Donor %d does not have hydrogen %d (a = %d)", d + 1, h + 1, a + 1);
}
}
else
}
add_ff(hb, id, k, ia, frame, ihb);
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
}
}
}
-static char *mkatomname(const t_atoms *atoms, int i)
+static char* mkatomname(const t_atoms* atoms, int i)
{
static char buf[32];
int rnr;
rnr = atoms->atom[i].resind;
- sprintf(buf, "%4s%d%-4s",
- *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
+ sprintf(buf, "%4s%d%-4s", *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]);
return buf;
}
-static void gen_datable(int *index, int isize, unsigned char *datable, int natoms)
+static void gen_datable(int* index, int isize, unsigned char* datable, int natoms)
{
/* Generates table of all atoms and sets the ingroup bit for atoms in index[] */
int i;
}
}
-static void clear_datable_grp(unsigned char *datable, int size)
+static void clear_datable_grp(unsigned char* datable, int size)
{
/* Clears group information from the table */
- int i;
+ int i;
if (size > 0)
{
for (i = 0; i < size; i++)
}
}
-static void add_acc(t_acceptors *a, int ia, int grp)
+static void add_acc(t_acceptors* a, int ia, int grp)
{
if (a->nra >= a->max_nra)
{
a->acc[a->nra++] = ia;
}
-static void search_acceptors(const t_topology *top, int isize,
- const int *index, t_acceptors *a, int grp,
- gmx_bool bNitAcc,
- gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable)
+static void search_acceptors(const t_topology* top,
+ int isize,
+ const int* index,
+ t_acceptors* a,
+ int grp,
+ gmx_bool bNitAcc,
+ gmx_bool bContact,
+ gmx_bool bDoIt,
+ unsigned char* datable)
{
int i, n;
for (i = 0; (i < isize); i++)
{
n = index[i];
- if ((bContact ||
- (((*top->atoms.atomname[n])[0] == 'O') ||
- (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
- ISINGRP(datable[n]))
+ if ((bContact
+ || (((*top->atoms.atomname[n])[0] == 'O')
+ || (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N'))))
+ && ISINGRP(datable[n]))
{
datable[n] |= c_acceptorMask;
add_acc(a, n, grp);
}
}
-static void add_h2d(int id, int ih, t_donors *ddd)
+static void add_h2d(int id, int ih, t_donors* ddd)
{
int i;
{
if (ddd->hydro[id][i] == ih)
{
- printf("Hm. This isn't the first time I found this donor (%d,%d)\n",
- ddd->don[id], ih);
+ printf("Hm. This isn't the first time I found this donor (%d,%d)\n", ddd->don[id], ih);
break;
}
}
{
if (ddd->nhydro[id] >= MAXHYDRO)
{
- gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!",
- ddd->don[id], MAXHYDRO);
+ gmx_fatal(FARGS, "Donor %d has more than %d hydrogens!", ddd->don[id], MAXHYDRO);
}
ddd->hydro[id][i] = ih;
ddd->nhydro[id]++;
}
}
-static void add_dh(t_donors *ddd, int id, int ih, int grp, const unsigned char *datable)
+static void add_dh(t_donors* ddd, int id, int ih, int grp, const unsigned char* datable)
{
int i;
if (!datable || ISDON(datable[id]))
{
- if (ddd->dptr[id] == NOTSET) /* New donor */
+ if (ddd->dptr[id] == NOTSET) /* New donor */
{
i = ddd->nrd;
ddd->dptr[id] = i;
}
add_h2d(i, ih, ddd);
}
- else
- if (datable)
+ else if (datable)
{
printf("Warning: Atom %d is not in the d/a-table!\n", id);
}
}
-static void search_donors(const t_topology *top, int isize, const int *index,
- t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt,
- unsigned char *datable)
+static void search_donors(const t_topology* top,
+ int isize,
+ const int* index,
+ t_donors* ddd,
+ int grp,
+ gmx_bool bContact,
+ gmx_bool bDoIt,
+ unsigned char* datable)
{
- int i, j;
- t_functype func_type;
- int nr1, nr2, nr3;
+ int i, j;
+ t_functype func_type;
+ int nr1, nr2, nr3;
if (!ddd->dptr)
{
{
for (func_type = 0; (func_type < F_NRE); func_type++)
{
- const t_ilist *interaction = &(top->idef.il[func_type]);
+ const t_ilist* interaction = &(top->idef.il[func_type]);
if (func_type == F_POSRES || func_type == F_FBPOSRES)
{
/* The ilist looks strange for posre. Bug in grompp?
continue;
}
for (i = 0; i < interaction->nr;
- i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
+ i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms + 1)
{
/* next function */
if (func_type != top->idef.functype[interaction->iatoms[i]])
{
- fprintf(stderr, "Error in func_type %s",
- interaction_function[func_type].longname);
+ fprintf(stderr, "Error in func_type %s", interaction_function[func_type].longname);
continue;
}
/* check out this functype */
if (func_type == F_SETTLE)
{
- nr1 = interaction->iatoms[i+1];
- nr2 = interaction->iatoms[i+2];
- nr3 = interaction->iatoms[i+3];
+ nr1 = interaction->iatoms[i + 1];
+ nr2 = interaction->iatoms[i + 2];
+ nr3 = interaction->iatoms[i + 3];
if (ISINGRP(datable[nr1]))
{
if (ISINGRP(datable[nr2]))
{
datable[nr1] |= c_donorMask;
- add_dh(ddd, nr1, nr1+1, grp, datable);
+ add_dh(ddd, nr1, nr1 + 1, grp, datable);
}
if (ISINGRP(datable[nr3]))
{
datable[nr1] |= c_donorMask;
- add_dh(ddd, nr1, nr1+2, grp, datable);
+ add_dh(ddd, nr1, nr1 + 2, grp, datable);
}
}
}
{
for (j = 0; j < 2; j++)
{
- nr1 = interaction->iatoms[i+1+j];
- nr2 = interaction->iatoms[i+2-j];
- if ((*top->atoms.atomname[nr1][0] == 'H') &&
- ((*top->atoms.atomname[nr2][0] == 'O') ||
- (*top->atoms.atomname[nr2][0] == 'N')) &&
- ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
+ nr1 = interaction->iatoms[i + 1 + j];
+ nr2 = interaction->iatoms[i + 2 - j];
+ if ((*top->atoms.atomname[nr1][0] == 'H')
+ && ((*top->atoms.atomname[nr2][0] == 'O') || (*top->atoms.atomname[nr2][0] == 'N'))
+ && ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
{
datable[nr2] |= c_donorMask;
add_dh(ddd, nr2, nr1, grp, datable);
{
interaction = &top->idef.il[func_type];
for (i = 0; i < interaction->nr;
- i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1)
+ i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms + 1)
{
/* next function */
if (func_type != top->idef.functype[interaction->iatoms[i]])
if (interaction_function[func_type].flags & IF_VSITE)
{
- nr1 = interaction->iatoms[i+1];
- if (*top->atoms.atomname[nr1][0] == 'H')
+ nr1 = interaction->iatoms[i + 1];
+ if (*top->atoms.atomname[nr1][0] == 'H')
{
- nr2 = nr1-1;
+ nr2 = nr1 - 1;
stop = FALSE;
- while (!stop && ( *top->atoms.atomname[nr2][0] == 'H'))
+ while (!stop && (*top->atoms.atomname[nr2][0] == 'H'))
{
if (nr2)
{
stop = TRUE;
}
}
- if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') ||
- ( *top->atoms.atomname[nr2][0] == 'N') ) &&
- ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
+ if (!stop
+ && ((*top->atoms.atomname[nr2][0] == 'O') || (*top->atoms.atomname[nr2][0] == 'N'))
+ && ISINGRP(datable[nr1]) && ISINGRP(datable[nr2]))
{
datable[nr2] |= c_donorMask;
add_dh(ddd, nr2, nr1, grp, datable);
}
}
-static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
+static t_gridcell*** init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid)
{
- t_gridcell ***grid;
+ t_gridcell*** grid;
int i, y, z;
if (bBox)
{
for (i = 0; i < DIM; i++)
{
- ngrid[i] = static_cast<int>(box[i][i]/(1.2*rcut));
+ ngrid[i] = static_cast<int>(box[i][i] / (1.2 * rcut));
}
}
- if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3) )
+ if (!bBox || (ngrid[XX] < 3) || (ngrid[YY] < 3) || (ngrid[ZZ] < 3))
{
for (i = 0; i < DIM; i++)
{
}
else
{
- printf("\nWill do grid-search on %dx%dx%d grid, rcut=%3.8f\n",
- ngrid[XX], ngrid[YY], ngrid[ZZ], rcut);
+ printf("\nWill do grid-search on %dx%dx%d grid, rcut=%3.8f\n", ngrid[XX], ngrid[YY],
+ ngrid[ZZ], rcut);
}
- if (((ngrid[XX]*ngrid[YY]*ngrid[ZZ]) * sizeof(grid)) > INT_MAX)
+ if (((ngrid[XX] * ngrid[YY] * ngrid[ZZ]) * sizeof(grid)) > INT_MAX)
{
- gmx_fatal(FARGS, "Failed to allocate memory for %d x %d x %d grid points, which is larger than the maximum of %zu. "
- "You are likely either using a box that is too large (box dimensions are %3.8f nm x%3.8f nm x%3.8f nm) or a cutoff (%3.8f nm) that is too small.",
- ngrid[XX], ngrid[YY], ngrid[ZZ], INT_MAX/sizeof(grid), box[XX][XX], box[YY][YY], box[ZZ][ZZ], rcut);
+ gmx_fatal(FARGS,
+ "Failed to allocate memory for %d x %d x %d grid points, which is larger than "
+ "the maximum of %zu. "
+ "You are likely either using a box that is too large (box dimensions are %3.8f "
+ "nm x%3.8f nm x%3.8f nm) or a cutoff (%3.8f nm) that is too small.",
+ ngrid[XX], ngrid[YY], ngrid[ZZ], INT_MAX / sizeof(grid), box[XX][XX], box[YY][YY],
+ box[ZZ][ZZ], rcut);
}
snew(grid, ngrid[ZZ]);
for (z = 0; z < ngrid[ZZ]; z++)
return grid;
}
-static void reset_nhbonds(t_donors *ddd)
+static void reset_nhbonds(t_donors* ddd)
{
int i, j;
static void pbc_correct_gem(rvec dx, matrix box, const rvec hbox);
static void pbc_in_gridbox(rvec dx, matrix box);
-static void build_grid(t_hbdata *hb, rvec x[], rvec xshell,
- gmx_bool bBox, matrix box, rvec hbox,
- real rcut, real rshell,
- ivec ngrid, t_gridcell ***grid)
+static void build_grid(t_hbdata* hb,
+ rvec x[],
+ rvec xshell,
+ gmx_bool bBox,
+ matrix box,
+ rvec hbox,
+ real rcut,
+ real rshell,
+ ivec ngrid,
+ t_gridcell*** grid)
{
- int i, m, gr, xi, yi, zi, nr;
- int *ad;
- ivec grididx;
- rvec invdelta, dshell;
- t_ncell *newgrid;
- gmx_bool bDoRshell, bInShell;
- real rshell2 = 0;
- int gx, gy, gz;
- int dum = -1;
+ int i, m, gr, xi, yi, zi, nr;
+ int* ad;
+ ivec grididx;
+ rvec invdelta, dshell;
+ t_ncell* newgrid;
+ gmx_bool bDoRshell, bInShell;
+ real rshell2 = 0;
+ int gx, gy, gz;
+ int dum = -1;
bDoRshell = (rshell > 0);
rshell2 = gmx::square(rshell);
bInShell = TRUE;
-#define DBB(x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__,#x, x)
+#define DBB(x) \
+ if (debug && bDebug) \
+ fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__, #x, x)
DBB(dum);
for (m = 0; m < DIM; m++)
{
- hbox[m] = box[m][m]*0.5;
+ hbox[m] = box[m][m] * 0.5;
if (bBox)
{
- invdelta[m] = ngrid[m]/box[m][m];
- if (1/invdelta[m] < rcut)
+ invdelta[m] = ngrid[m] / box[m][m];
+ if (1 / invdelta[m] < rcut)
{
- gmx_fatal(FARGS, "Your computational box has shrunk too much.\n"
+ gmx_fatal(FARGS,
+ "Your computational box has shrunk too much.\n"
"%s can not handle this situation, sorry.\n",
gmx::getProgramContext().displayName());
}
while (!bDone)
{
bDone = TRUE;
- for (m = DIM-1; m >= 0 && bInShell; m--)
+ for (m = DIM - 1; m >= 0 && bInShell; m--)
{
if (dshell[m] < -hbox[m])
{
}
if (dshell[m] >= hbox[m])
{
- bDone = FALSE;
- dshell[m] -= 2*hbox[m];
+ bDone = FALSE;
+ dshell[m] -= 2 * hbox[m];
}
}
}
- for (m = DIM-1; m >= 0 && bInShell; m--)
+ for (m = DIM - 1; m >= 0 && bInShell; m--)
{
/* if we're outside the cube, we're outside the sphere also! */
- if ( (dshell[m] > rshell) || (-dshell[m] > rshell) )
+ if ((dshell[m] > rshell) || (-dshell[m] > rshell))
{
bInShell = FALSE;
}
{
pbc_in_gridbox(x[ad[i]], box);
- for (m = DIM-1; m >= 0; m--)
- { /* determine grid index of atom */
- grididx[m] = static_cast<int>(x[ad[i]][m]*invdelta[m]);
- grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m];
+ for (m = DIM - 1; m >= 0; m--)
+ { /* determine grid index of atom */
+ grididx[m] = static_cast<int>(x[ad[i]][m] * invdelta[m]);
+ grididx[m] = (grididx[m] + ngrid[m]) % ngrid[m];
}
}
}
}
-static void count_da_grid(const ivec ngrid, t_gridcell ***grid, t_icell danr)
+static void count_da_grid(const ivec ngrid, t_gridcell*** grid, t_icell danr)
{
int gr, xi, yi, zi;
*/
static inline int grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
{
- return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
+ return ((n == 1) ? x : bTric && bEdge ? 0 : (x - 1));
}
static inline int grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
{
- return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
+ return ((n == 1) ? x : bTric && bEdge ? (n - 1) : (x + 1));
}
static inline int grid_mod(int j, int n)
{
- return (j+n) % (n);
+ return (j + n) % (n);
}
-static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
+static void dump_grid(FILE* fp, ivec ngrid, t_gridcell*** grid)
{
int gr, x, y, z, sum[grNR];
fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]);
for (z = 0; z < ngrid[ZZ]; z += 2)
{
- fprintf(fp, "Z=%d,%d\n", z, z+1);
+ fprintf(fp, "Z=%d,%d\n", z, z + 1);
for (y = 0; y < ngrid[YY]; y++)
{
for (x = 0; x < ngrid[XX]; x++)
sum[gr] += grid[z][y][x].d[gr].nr;
fprintf(fp, "%3d", grid[x][y][z].a[gr].nr);
sum[gr] += grid[z][y][x].a[gr].nr;
-
}
fprintf(fp, " | ");
- if ( (z+1) < ngrid[ZZ])
+ if ((z + 1) < ngrid[ZZ])
{
for (x = 0; x < ngrid[XX]; x++)
{
- fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr);
- sum[gr] += grid[z+1][y][x].d[gr].nr;
- fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr);
- sum[gr] += grid[z+1][y][x].a[gr].nr;
+ fprintf(fp, "%3d", grid[z + 1][y][x].d[gr].nr);
+ sum[gr] += grid[z + 1][y][x].d[gr].nr;
+ fprintf(fp, "%3d", grid[z + 1][y][x].a[gr].nr);
+ sum[gr] += grid[z + 1][y][x].a[gr].nr;
}
}
fprintf(fp, "\n");
/* New GMX record! 5 * in a row. Congratulations!
* Sorry, only four left.
*/
-static void free_grid(const ivec ngrid, t_gridcell ****grid)
+static void free_grid(const ivec ngrid, t_gridcell**** grid)
{
int y, z;
- t_gridcell ***g = *grid;
+ t_gridcell*** g = *grid;
for (z = 0; z < ngrid[ZZ]; z++)
{
while (!bDone)
{
bDone = TRUE;
- for (m = DIM-1; m >= 0; m--)
+ for (m = DIM - 1; m >= 0; m--)
{
if (dx[m] < -hbox[m])
{
while (!bDone)
{
bDone = TRUE;
- for (m = DIM-1; m >= 0; m--)
+ for (m = DIM - 1; m >= 0; m--)
{
if (dx[m] < 0)
{
* use of second cut-off.
* - Erik Marklund, June 29, 2006
*/
-static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a,
- real rcut, real r2cut, real ccut,
- rvec x[], gmx_bool bBox, matrix box, rvec hbox,
- real *d_ha, real *ang, gmx_bool bDA, int *hhh,
- gmx_bool bContact, gmx_bool bMerge)
+static int is_hbond(t_hbdata* hb,
+ int grpd,
+ int grpa,
+ int d,
+ int a,
+ real rcut,
+ real r2cut,
+ real ccut,
+ rvec x[],
+ gmx_bool bBox,
+ matrix box,
+ rvec hbox,
+ real* d_ha,
+ real* ang,
+ gmx_bool bDA,
+ int* hhh,
+ gmx_bool bContact,
+ gmx_bool bMerge)
{
int h, hh, id;
rvec r_da, r_ha, r_dh;
return hbNo;
}
- if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) ||
- (acceptor_index(&hb->a, grpa, a) == NOTSET))
+ if (((id = donor_index(&hb->d, grpd, d)) == NOTSET) || (acceptor_index(&hb->a, grpa, a) == NOTSET))
{
return hbNo;
}
- rc2 = rcut*rcut;
- r2c2 = r2cut*r2cut;
+ rc2 = rcut * rcut;
+ r2c2 = r2cut * r2cut;
rvec_sub(x[d], x[a], r_da);
/* Insert projection code here */
}
if (bBox)
{
- if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
- { /* return hbNo; */
- daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
+ if (d > a && bMerge
+ && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */
+ { /* return hbNo; */
+ daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
}
pbc_correct_gem(r_da, box, hbox);
}
for (h = 0; (h < hb->d.nhydro[id]); h++)
{
hh = hb->d.hydro[id][h];
- rha2 = rc2+1;
+ rha2 = rc2 + 1;
if (!bDA)
{
rvec_sub(x[hh], x[a], r_ha);
/* Merging is now done on the fly, so do_merge is most likely obsolete now.
* Will do some more testing before removing the function entirely.
* - Erik Marklund, MAY 10 2010 */
-static void do_merge(t_hbdata *hb, int ntmp,
- bool htmp[], bool gtmp[],
- t_hbond *hb0, t_hbond *hb1)
+static void do_merge(t_hbdata* hb, int ntmp, bool htmp[], bool gtmp[], t_hbond* hb0, t_hbond* hb1)
{
/* Here we need to make sure we're treating periodicity in
* the right way for the geminate recombination kinetics. */
- int m, mm, n00, n01, nn0, nnframes;
+ int m, mm, n00, n01, nn0, nnframes;
/* Decide where to start from when merging */
n00 = hb0->n0;
- Erik Marklund, June 1, 2006 */
for (m = 0; (m <= hb0->nframes); m++)
{
- mm = m+n00-nn0;
+ mm = m + n00 - nn0;
htmp[mm] = is_hb(hb0->h[0], m);
}
for (m = 0; (m <= hb0->nframes); m++)
{
- mm = m+n00-nn0;
+ mm = m + n00 - nn0;
gtmp[mm] = is_hb(hb0->g[0], m);
}
/* Next HB */
for (m = 0; (m <= hb1->nframes); m++)
{
- mm = m+n01-nn0;
+ mm = m + n01 - nn0;
htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m);
gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m);
}
/* Reallocate target array */
if (nnframes > hb0->maxframes)
{
- srenew(hb0->h[0], 4+nnframes/hb->wordlen);
- srenew(hb0->g[0], 4+nnframes/hb->wordlen);
+ srenew(hb0->h[0], 4 + nnframes / hb->wordlen);
+ srenew(hb0->g[0], 4 + nnframes / hb->wordlen);
}
/* Copy temp array to target array */
hb0->maxframes = nnframes;
}
-static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact)
+static void merge_hb(t_hbdata* hb, gmx_bool bTwo, gmx_bool bContact)
{
- int i, inrnew, indnew, j, ii, jj, id, ia, ntmp;
- bool *htmp, *gtmp;
- t_hbond *hb0, *hb1;
+ int i, inrnew, indnew, j, ii, jj, id, ia, ntmp;
+ bool * htmp, *gtmp;
+ t_hbond *hb0, *hb1;
inrnew = hb->nrhb;
indnew = hb->nrdist;
/* Check whether donors are also acceptors */
printf("Merging hbonds with Acceptor and Donor swapped\n");
- ntmp = 2*hb->max_frames;
+ ntmp = 2 * hb->max_frames;
snew(gtmp, ntmp);
snew(htmp, ntmp);
for (i = 0; (i < hb->d.nrd); i++)
{
- fprintf(stderr, "\r%d/%d", i+1, hb->d.nrd);
+ fprintf(stderr, "\r%d/%d", i + 1, hb->d.nrd);
fflush(stderr);
id = hb->d.don[i];
ii = hb->a.aptr[id];
{
ia = hb->a.acc[j];
jj = hb->d.dptr[ia];
- if ((id != ia) && (ii != NOTSET) && (jj != NOTSET) &&
- (!bTwo || (hb->d.grp[i] != hb->a.grp[j])))
+ if ((id != ia) && (ii != NOTSET) && (jj != NOTSET)
+ && (!bTwo || (hb->d.grp[i] != hb->a.grp[j])))
{
hb0 = hb->hbmap[i][j];
hb1 = hb->hbmap[jj][ii];
{
indnew--;
}
- else
- if (bContact)
+ else if (bContact)
{
gmx_incons("No contact history");
}
sfree(htmp);
}
-static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t)
+static void do_nhb_dist(FILE* fp, t_hbdata* hb, real t)
{
- int i, j, k, n_bound[MAXHH], nbtot;
+ int i, j, k, n_bound[MAXHH], nbtot;
/* Set array to 0 */
for (k = 0; (k < MAXHH); k++)
for (k = 0; (k < MAXHH); k++)
{
fprintf(fp, " %8d", n_bound[k]);
- nbtot += n_bound[k]*k;
+ nbtot += n_bound[k] * k;
}
fprintf(fp, " %8d\n", nbtot);
}
-static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact,
- const gmx_output_env_t *oenv)
+static void do_hblife(const char* fn, t_hbdata* hb, gmx_bool bMerge, gmx_bool bContact, const gmx_output_env_t* oenv)
{
- FILE *fp;
- const char *leg[] = { "p(t)", "t p(t)" };
- int *histo;
+ FILE* fp;
+ const char* leg[] = { "p(t)", "t p(t)" };
+ int* histo;
int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0;
int nframes = hb->nframes;
- unsigned int **h;
+ unsigned int** h;
real t, x1, dt;
double sum, integral;
- t_hbond *hbh;
+ t_hbond* hbh;
snew(h, hb->maxhydro);
- snew(histo, nframes+1);
+ snew(histo, nframes + 1);
/* Total number of hbonds analyzed here */
for (i = 0; (i < hb->d.nrd); i++)
{
for (j = 0; (j <= hbh->nframes); j++)
{
- ihb = static_cast<int>(is_hb(h[nh], j));
+ ihb = static_cast<int>(is_hb(h[nh], j));
if (debug && (ndump < 10))
{
fprintf(debug, "%5d %5d\n", j, ihb);
}
else
{
- histo[j-j0]++;
+ histo[j - j0]++;
}
ohb = ihb;
}
}
else
{
- fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()",
- oenv);
+ fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv),
+ "()", oenv);
}
xvgr_legend(fp, asize(leg), leg, oenv);
- j0 = nframes-1;
+ j0 = nframes - 1;
while ((j0 > 0) && (histo[j0] == 0))
{
j0--;
{
sum += histo[i];
}
- dt = hb->time[1]-hb->time[0];
- sum = dt*sum;
+ dt = hb->time[1] - hb->time[0];
+ sum = dt * sum;
integral = 0;
for (i = 1; (i <= j0); i++)
{
- t = hb->time[i] - hb->time[0] - 0.5*dt;
- x1 = t*histo[i]/sum;
- fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1);
+ t = hb->time[i] - hb->time[0] - 0.5 * dt;
+ x1 = t * histo[i] / sum;
+ fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i] / sum, x1);
integral += x1;
}
integral *= dt;
sfree(histo);
}
-static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump)
+static void dump_ac(t_hbdata* hb, gmx_bool oneHB, int nDump)
{
- FILE *fp;
- int i, j, k, m, nd, ihb, idist;
- int nframes = hb->nframes;
- gmx_bool bPrint;
- t_hbond *hbh;
+ FILE* fp;
+ int i, j, k, m, nd, ihb, idist;
+ int nframes = hb->nframes;
+ gmx_bool bPrint;
+ t_hbond* hbh;
if (nDump <= 0)
{
for (k = 0; (k < hb->a.nra) && (nd < nDump); k++)
{
bPrint = FALSE;
- ihb = idist = 0;
- hbh = hb->hbmap[i][k];
+ ihb = idist = 0;
+ hbh = hb->hbmap[i][k];
if (oneHB)
{
if (hbh->h[0])
{
for (m = 0; (m < hb->maxhydro) && !ihb; m++)
{
- ihb = static_cast<int>((ihb != 0) || (((hbh->h[m]) != nullptr) && is_hb(hbh->h[m], j)));
- idist = static_cast<int>((idist != 0) || (((hbh->g[m]) != nullptr) && is_hb(hbh->g[m], j)));
+ ihb = static_cast<int>((ihb != 0)
+ || (((hbh->h[m]) != nullptr) && is_hb(hbh->h[m], j)));
+ idist = static_cast<int>((idist != 0)
+ || (((hbh->g[m]) != nullptr) && is_hb(hbh->g[m], j)));
}
/* This is not correct! */
/* What isn't correct? -Erik M */
{
real kbt;
- kbt = BOLTZ*temp;
+ kbt = BOLTZ * temp;
if (tau <= 0)
{
return -666;
}
else
{
- return kbt*std::log(kbt*tau/PLANCK);
+ return kbt * std::log(kbt * tau / PLANCK);
}
}
-typedef struct {
+typedef struct
+{
int n0, n1, nparams, ndelta;
real kkk[2];
real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
} t_luzar;
-static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
- real kt[], real sigma_ct[], real sigma_nt[],
- real sigma_kt[], real *k, real *kp,
- real *sigma_k, real *sigma_kp,
- real fit_start)
+static real compute_weighted_rates(int n,
+ real t[],
+ real ct[],
+ real nt[],
+ real kt[],
+ real sigma_ct[],
+ real sigma_nt[],
+ real sigma_kt[],
+ real* k,
+ real* kp,
+ real* sigma_k,
+ real* sigma_kp,
+ real fit_start)
{
#define NK 10
- int i, j;
- t_luzar tl;
- real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
+ int i, j;
+ t_luzar tl;
+ real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2;
*sigma_k = 0;
*sigma_kp = 0;
tl.kkk[0] = *k;
tl.kkk[1] = *kp;
- chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
- *k = tl.kkk[0];
- *kp = tl.kkk[1] = *kp;
- tl.ndelta = NK;
+ chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
+ *k = tl.kkk[0];
+ *kp = tl.kkk[1] = *kp;
+ tl.ndelta = NK;
for (j = 0; (j < NK); j++)
{
kkk += tl.kkk[0];
kp2 += gmx::square(tl.kkk[1]);
tl.n0++;
}
- *sigma_k = std::sqrt(kk2/NK - gmx::square(kkk/NK));
- *sigma_kp = std::sqrt(kp2/NK - gmx::square(kkp/NK));
+ *sigma_k = std::sqrt(kk2 / NK - gmx::square(kkk / NK));
+ *sigma_kp = std::sqrt(kp2 / NK - gmx::square(kkp / NK));
return chi2;
}
-void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
- real sigma_ct[], real sigma_nt[], real sigma_kt[],
- real fit_start, real temp)
+void analyse_corr(int n,
+ real t[],
+ real ct[],
+ real nt[],
+ real kt[],
+ real sigma_ct[],
+ real sigma_nt[],
+ real sigma_kt[],
+ real fit_start,
+ real temp)
{
- int i0, i;
- real k = 1, kp = 1, kow = 1;
- real Q = 0, chi2, tau_hb, dtau, tau_rlx, e_1, sigma_k, sigma_kp, ddg;
- double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
- gmx_bool bError = (sigma_ct != nullptr) && (sigma_nt != nullptr) && (sigma_kt != nullptr);
+ int i0, i;
+ real k = 1, kp = 1, kow = 1;
+ real Q = 0, chi2, tau_hb, dtau, tau_rlx, e_1, sigma_k, sigma_kp, ddg;
+ double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0;
+ gmx_bool bError = (sigma_ct != nullptr) && (sigma_nt != nullptr) && (sigma_kt != nullptr);
- for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++)
- {
- ;
- }
- if (i0 < n-2)
+ for (i0 = 0; (i0 < n - 2) && ((t[i0] - t[0]) < fit_start); i0++) {}
+ if (i0 < n - 2)
{
for (i = i0; (i < n); i++)
{
sc2 += gmx::square(ct[i]);
sn2 += gmx::square(nt[i]);
sk2 += gmx::square(kt[i]);
- sck += ct[i]*kt[i];
- snk += nt[i]*kt[i];
- scn += ct[i]*nt[i];
+ sck += ct[i] * kt[i];
+ snk += nt[i] * kt[i];
+ scn += ct[i] * nt[i];
}
printf("Hydrogen bond thermodynamics at T = %g K\n", temp);
- tmp = (sn2*sc2-gmx::square(scn));
+ tmp = (sn2 * sc2 - gmx::square(scn));
if ((tmp > 0) && (sn2 > 0))
{
- k = (sn2*sck-scn*snk)/tmp;
- kp = (k*scn-snk)/sn2;
+ k = (sn2 * sck - scn * snk) / tmp;
+ kp = (k * scn - snk) / sn2;
if (bError)
{
chi2 = 0;
for (i = i0; (i < n); i++)
{
- chi2 += gmx::square(k*ct[i]-kp*nt[i]-kt[i]);
+ chi2 += gmx::square(k * ct[i] - kp * nt[i] - kt[i]);
}
- compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
- sigma_kt, &k, &kp,
+ compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt, sigma_kt, &k, &kp,
&sigma_k, &sigma_kp, fit_start);
Q = 0; /* quality_of_fit(chi2, 2);*/
- ddg = BOLTZ*temp*sigma_k/k;
- printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
- chi2, Q);
+ ddg = BOLTZ * temp * sigma_k / k;
+ printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n", chi2, Q);
printf("The Rate and Delta G are followed by an error estimate\n");
printf("----------------------------------------------------------\n"
"Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n");
- printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
- k, sigma_k, 1/k, calc_dg(1/k, temp), ddg);
- ddg = BOLTZ*temp*sigma_kp/kp;
- printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n",
- kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg);
+ printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n", k, sigma_k, 1 / k,
+ calc_dg(1 / k, temp), ddg);
+ ddg = BOLTZ * temp * sigma_kp / kp;
+ printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n", kp, sigma_kp, 1 / kp,
+ calc_dg(1 / kp, temp), ddg);
}
else
{
chi2 = 0;
for (i = i0; (i < n); i++)
{
- chi2 += gmx::square(k*ct[i]-kp*nt[i]-kt[i]);
+ chi2 += gmx::square(k * ct[i] - kp * nt[i] - kt[i]);
}
- printf("Fitting parameters chi^2 = %10g\nQ = %10g\n",
- chi2, Q);
+ printf("Fitting parameters chi^2 = %10g\nQ = %10g\n", chi2, Q);
printf("--------------------------------------------------\n"
"Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n");
- printf("Forward %10.3f %8.3f %10.3f %10g\n",
- k, 1/k, calc_dg(1/k, temp), chi2);
- printf("Backward %10.3f %8.3f %10.3f\n",
- kp, 1/kp, calc_dg(1/kp, temp));
+ printf("Forward %10.3f %8.3f %10.3f %10g\n", k, 1 / k, calc_dg(1 / k, temp), chi2);
+ printf("Backward %10.3f %8.3f %10.3f\n", kp, 1 / kp, calc_dg(1 / kp, temp));
}
}
if (sc2 > 0)
{
- kow = 2*sck/sc2;
- printf("One-way %10.3f %s%8.3f %10.3f\n",
- kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp));
+ kow = 2 * sck / sc2;
+ printf("One-way %10.3f %s%8.3f %10.3f\n", kow, bError ? " " : "", 1 / kow,
+ calc_dg(1 / kow, temp));
}
else
{
sc2, sn2, sk2, sck, snk, scn);
}
/* Determine integral of the correlation function */
- tau_hb = evaluate_integral(n, t, ct, nullptr, (t[n-1]-t[0])/2, &dtau);
- printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb,
- bError ? " " : "", tau_hb, calc_dg(tau_hb, temp));
+ tau_hb = evaluate_integral(n, t, ct, nullptr, (t[n - 1] - t[0]) / 2, &dtau);
+ printf("Integral %10.3f %s%8.3f %10.3f\n", 1 / tau_hb, bError ? " " : "", tau_hb,
+ calc_dg(tau_hb, temp));
e_1 = std::exp(-1.0);
- for (i = 0; (i < n-2); i++)
+ for (i = 0; (i < n - 2); i++)
{
- if ((ct[i] > e_1) && (ct[i+1] <= e_1))
+ if ((ct[i] > e_1) && (ct[i + 1] <= e_1))
{
break;
}
}
- if (i < n-2)
+ if (i < n - 2)
{
/* Determine tau_relax from linear interpolation */
- tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]);
- printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx,
- tau_rlx, bError ? " " : "",
- calc_dg(tau_rlx, temp));
+ tau_rlx = t[i] - t[0] + (e_1 - ct[i]) * (t[i + 1] - t[i]) / (ct[i + 1] - ct[i]);
+ printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1 / tau_rlx, tau_rlx,
+ bError ? " " : "", calc_dg(tau_rlx, temp));
}
}
else
int j;
/* Compute k(t) = dc(t)/dt */
- for (j = 1; (j < nn-1); j++)
+ for (j = 1; (j < nn - 1); j++)
{
- dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]);
+ dydx[j] = (y[j + 1] - y[j - 1]) / (x[j + 1] - x[j - 1]);
}
/* Extrapolate endpoints */
- dydx[0] = 2*dydx[1] - dydx[2];
- dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3];
+ dydx[0] = 2 * dydx[1] - dydx[2];
+ dydx[nn - 1] = 2 * dydx[nn - 2] - dydx[nn - 3];
}
-static void normalizeACF(real *ct, real *gt, int nhb, int len)
+static void normalizeACF(real* ct, real* gt, int nhb, int len)
{
real ct_fac, gt_fac = 0;
int i;
/* Xu and Berne use the same normalization constant */
- ct_fac = 1.0/ct[0];
+ ct_fac = 1.0 / ct[0];
if (nhb != 0)
{
- gt_fac = 1.0/nhb;
+ gt_fac = 1.0 / nhb;
}
printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
}
}
-static void do_hbac(const char *fn, t_hbdata *hb,
- int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start,
- real temp, gmx_bool R2, const gmx_output_env_t *oenv,
- int nThreads)
+static void do_hbac(const char* fn,
+ t_hbdata* hb,
+ int nDump,
+ gmx_bool bMerge,
+ gmx_bool bContact,
+ real fit_start,
+ real temp,
+ gmx_bool R2,
+ const gmx_output_env_t* oenv,
+ int nThreads)
{
- FILE *fp;
- int i, j, k, m, ihb, idist, n2, nn;
-
- const char *legLuzar[] = {
- "Ac\\sfin sys\\v{}\\z{}(t)",
- "Ac(t)",
- "Cc\\scontact,hb\\v{}\\z{}(t)",
- "-dAc\\sfs\\v{}\\z{}/dt"
- };
- gmx_bool bNorm = FALSE;
- double nhb = 0;
- real *rhbex = nullptr, *ht, *gt, *ght, *dght, *kt;
- real *ct, tail, tail2, dtail, *cct;
- const real tol = 1e-3;
- int nframes = hb->nframes;
- unsigned int **h = nullptr, **g = nullptr;
+ FILE* fp;
+ int i, j, k, m, ihb, idist, n2, nn;
+
+ const char* legLuzar[] = { "Ac\\sfin sys\\v{}\\z{}(t)", "Ac(t)", "Cc\\scontact,hb\\v{}\\z{}(t)",
+ "-dAc\\sfs\\v{}\\z{}/dt" };
+ gmx_bool bNorm = FALSE;
+ double nhb = 0;
+ real * rhbex = nullptr, *ht, *gt, *ght, *dght, *kt;
+ real * ct, tail, tail2, dtail, *cct;
+ const real tol = 1e-3;
+ int nframes = hb->nframes;
+ unsigned int **h = nullptr, **g = nullptr;
int nh, nhbonds, nhydro;
- t_hbond *hbh;
+ t_hbond* hbh;
int acType;
- int *dondata = nullptr;
+ int* dondata = nullptr;
- enum {
- AC_NONE, AC_NN, AC_GEM, AC_LUZAR
+ enum
+ {
+ AC_NONE,
+ AC_NN,
+ AC_GEM,
+ AC_LUZAR
};
const bool bOMP = GMX_OPENMP;
n2 *= 2;
}
- nn = nframes/2;
+ nn = nframes / 2;
if (acType != AC_NN || bOMP)
{
dondata[i] = -1;
}
printf("ACF calculations parallelized with OpenMP using %i threads.\n"
- "Expect close to linear scaling over this donor-loop.\n", nThreads);
+ "Expect close to linear scaling over this donor-loop.\n",
+ nThreads);
fflush(stdout);
fprintf(stderr, "Donors: [thread no]\n");
{
/* Build the ACF */
- snew(rhbex, 2*n2);
- snew(ct, 2*n2);
- snew(gt, 2*n2);
- snew(ht, 2*n2);
- snew(ght, 2*n2);
- snew(dght, 2*n2);
+ snew(rhbex, 2 * n2);
+ snew(ct, 2 * n2);
+ snew(gt, 2 * n2);
+ snew(ht, 2 * n2);
+ snew(ght, 2 * n2);
+ snew(dght, 2 * n2);
snew(kt, nn);
snew(cct, nn);
for (nh = 0; (nh < nhydro); nh++)
{
int nrint = bContact ? hb->nrdist : hb->nrhb;
- if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint))
+ if ((((nhbonds + 1) % 10) == 0) || (nhbonds + 1 == nrint))
{
- fprintf(stderr, "\rACF %d/%d", nhbonds+1, nrint);
+ fprintf(stderr, "\rACF %d/%d", nhbonds + 1, nrint);
fflush(stderr);
}
nhbonds++;
* otherwise use g(t) = 1-h(t) */
if (!R2 && bContact)
{
- gt[j] = 1-ihb;
+ gt[j] = 1 - ihb;
}
else
{
- gt[j] = idist*(1-ihb);
+ gt[j] = idist * (1 - ihb);
}
- ht[j] = rhbex[j];
- nhb += ihb;
+ ht[j] = rhbex[j];
+ nhb += ihb;
}
/* The autocorrelation function is normalized after summation only */
- low_do_autocorr(nullptr, oenv, nullptr, nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0],
- eacNormal, 1, FALSE, bNorm, FALSE, 0, -1, 0);
+ low_do_autocorr(nullptr, oenv, nullptr, nframes, 1, -1, &rhbex,
+ hb->time[1] - hb->time[0], eacNormal, 1, FALSE, bNorm, FALSE, 0,
+ -1, 0);
/* Cross correlation analysis for thermodynamics */
for (j = nframes; (j < n2); j++)
for (j = 0; (j < nn); j++)
{
- ct[j] += rhbex[j];
+ ct[j] += rhbex[j];
ght[j] += dght[j];
}
}
/* Determine tail value for statistics */
tail = 0;
tail2 = 0;
- for (j = nn/2; (j < nn); j++)
+ for (j = nn / 2; (j < nn); j++)
{
- tail += ct[j];
- tail2 += ct[j]*ct[j];
+ tail += ct[j];
+ tail2 += ct[j] * ct[j];
}
- tail /= (nn - int{nn/2});
- tail2 /= (nn - int{nn/2});
- dtail = std::sqrt(tail2-tail*tail);
+ tail /= (nn - int{ nn / 2 });
+ tail2 /= (nn - int{ nn / 2 });
+ dtail = std::sqrt(tail2 - tail * tail);
/* Check whether the ACF is long enough */
if (dtail > tol)
for (j = 0; (j < nn); j++)
{
cct[j] = ct[j];
- ct[j] = (cct[j]-tail)/(1-tail);
+ ct[j] = (cct[j] - tail) / (1 - tail);
}
/* Compute negative derivative k(t) = -dc(t)/dt */
compute_derivative(nn, hb->time, ct, kt);
for (j = 0; (j < nn); j++)
{
- fprintf(fp, "%10g %10g %10g %10g %10g\n",
- hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
+ fprintf(fp, "%10g %10g %10g %10g %10g\n", hb->time[j] - hb->time[0], ct[j], cct[j],
+ ght[j], kt[j]);
}
xvgrclose(fp);
- analyse_corr(nn, hb->time, ct, ght, kt, nullptr, nullptr, nullptr,
- fit_start, temp);
+ analyse_corr(nn, hb->time, ct, ght, kt, nullptr, nullptr, nullptr, fit_start, temp);
do_view(oenv, fn, nullptr);
sfree(rhbex);
sfree(kt);
}
-static void init_hbframe(t_hbdata *hb, int nframes, real t)
+static void init_hbframe(t_hbdata* hb, int nframes, real t)
{
int i;
- hb->time[nframes] = t;
- hb->nhb[nframes] = 0;
- hb->ndist[nframes] = 0;
+ hb->time[nframes] = t;
+ hb->nhb[nframes] = 0;
+ hb->ndist[nframes] = 0;
for (i = 0; (i < max_hx); i++)
{
hb->nhx[nframes][i] = 0;
}
}
-static FILE *open_donor_properties_file(const char *fn,
- t_hbdata *hb,
- const gmx_output_env_t *oenv)
+static FILE* open_donor_properties_file(const char* fn, t_hbdata* hb, const gmx_output_env_t* oenv)
{
- FILE *fp = nullptr;
- const char *leg[] = { "Nbound", "Nfree" };
+ FILE* fp = nullptr;
+ const char* leg[] = { "Nbound", "Nfree" };
if (!fn || !hb)
{
return fp;
}
-static void analyse_donor_properties(FILE *fp, t_hbdata *hb, int nframes, real t)
+static void analyse_donor_properties(FILE* fp, t_hbdata* hb, int nframes, real t)
{
int i, j, k, nbound, nb, nhtot;
nhtot++;
for (j = 0; (j < hb->a.nra) && (nb == 0); j++)
{
- if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] &&
- is_hb(hb->hbmap[i][j]->h[k], nframes))
+ if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] && is_hb(hb->hbmap[i][j]->h[k], nframes))
{
nb = 1;
}
nbound += nb;
}
}
- fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound);
+ fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot - nbound);
}
-static void dump_hbmap(t_hbdata *hb,
- int nfile, t_filenm fnm[], gmx_bool bTwo,
- gmx_bool bContact, const int isize[], int *index[], char *grpnames[],
- const t_atoms *atoms)
+static void dump_hbmap(t_hbdata* hb,
+ int nfile,
+ t_filenm fnm[],
+ gmx_bool bTwo,
+ gmx_bool bContact,
+ const int isize[],
+ int* index[],
+ char* grpnames[],
+ const t_atoms* atoms)
{
- FILE *fp, *fplog;
+ FILE * fp, *fplog;
int ddd, hhh, aaa, i, j, k, m, grp;
char ds[32], hs[32], as[32];
gmx_bool first;
fprintf(fp, "[ %s ]", grpnames[grp]);
for (i = 0; i < isize[grp]; i++)
{
- fprintf(fp, (i%15) ? " " : "\n");
- fprintf(fp, " %4d", index[grp][i]+1);
+ fprintf(fp, (i % 15) ? " " : "\n");
+ fprintf(fp, " %4d", index[grp][i] + 1);
}
fprintf(fp, "\n");
{
for (j = 0; (j < hb->d.nhydro[i]); j++)
{
- fprintf(fp, " %4d %4d", hb->d.don[i]+1,
- hb->d.hydro[i][j]+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, " %4d", hb->a.acc[i]+1);
+ fprintf(fp, (i % 15 && !first) ? " " : "\n");
+ fprintf(fp, " %4d", hb->a.acc[i] + 1);
first = FALSE;
}
}
}
if (bTwo)
{
- fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" :
- "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
+ fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" : "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]);
}
else
{
sprintf(as, "%s", mkatomname(atoms, aaa));
if (bContact)
{
- fprintf(fp, " %6d %6d\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, " %6d %6d %6d\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);
/* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
* It mimics add_frames() and init_frame() to some extent. */
-static void sync_hbdata(t_hbdata *p_hb, int nframes)
+static void sync_hbdata(t_hbdata* p_hb, int nframes)
{
if (nframes >= p_hb->max_frames)
{
p_hb->max_frames += 4096;
- srenew(p_hb->nhb, p_hb->max_frames);
+ srenew(p_hb->nhb, p_hb->max_frames);
srenew(p_hb->ndist, p_hb->max_frames);
srenew(p_hb->n_bound, p_hb->max_frames);
srenew(p_hb->nhx, p_hb->max_frames);
{
srenew(p_hb->danr, p_hb->max_frames);
}
- std::memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
- std::memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes));
+ std::memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames - nframes));
+ std::memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames - nframes));
p_hb->nhb[nframes] = 0;
p_hb->ndist[nframes] = 0;
-
}
p_hb->nframes = nframes;
- std::memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx); /* zero the helix count for this frame */
+ std::memset(&(p_hb->nhx[nframes]), 0, sizeof(int) * max_hx); /* zero the helix count for this frame */
}
-int gmx_hbond(int argc, char *argv[])
+int gmx_hbond(int argc, char* argv[])
{
- const char *desc[] = {
+ const char* desc[] = {
"[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are",
"determined based on cutoffs for the angle Hydrogen - Donor - Acceptor",
"(zero is extended) and the distance Donor - Acceptor",
"If you set [TT]-shell[tt], you will be asked for an additional index group",
"which should contain exactly one atom. In this case, only hydrogen",
- "bonds between atoms within the shell distance from the one atom are",
- "considered.[PAR]",
+ "bonds between atoms within the shell distance from the one atom are", "considered.[PAR]",
"With option -ac, rate constants for hydrogen bonding can be derived with the",
"model of Luzar and Chandler (Nature 379:55, 1996; J. Chem. Phys. 113:23, 2000).",
"n(t) can be defined as either all pairs that are not within contact distance r at time t",
"(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
"are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
- "See mentioned literature for more details and definitions.",
- "[PAR]",
+ "See mentioned literature for more details and definitions.", "[PAR]",
/* "It is also possible to analyse specific hydrogen bonds with",
"[TT]-sel[tt]. This index file must contain a group of atom triplets",
"note also that no check is made for the types of atoms.[PAR]",
*/
- "[BB]Output:[bb]",
- "",
- " * [TT]-num[tt]: number of hydrogen bonds as a function of time.",
+ "[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.",
" 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.",
+ " 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.",
- "",
+ " 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)."
};
- static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
- static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15;
- static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
- static int nDump = 0;
- static int nThreads = 0;
+ static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1;
+ static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15;
+ static gmx_bool bNitAcc = TRUE, bDA = TRUE, bMerge = TRUE;
+ static int nDump = 0;
+ static int nThreads = 0;
- static gmx_bool bContact = FALSE;
+ static gmx_bool bContact = FALSE;
/* options */
- t_pargs pa [] = {
- { "-a", FALSE, etREAL, {&acut},
- "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
- { "-r", FALSE, etREAL, {&rcut},
- "Cutoff radius (nm, X - Acceptor, see next option)" },
- { "-da", FALSE, etBOOL, {&bDA},
+ t_pargs pa[] = {
+ { "-a", FALSE, etREAL, { &acut }, "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" },
+ { "-r", FALSE, etREAL, { &rcut }, "Cutoff radius (nm, X - Acceptor, see next option)" },
+ { "-da",
+ FALSE,
+ etBOOL,
+ { &bDA },
"Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" },
- { "-r2", FALSE, etREAL, {&r2cut},
- "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"},
- { "-abin", FALSE, etREAL, {&abin},
- "Binwidth angle distribution (degrees)" },
- { "-rbin", FALSE, etREAL, {&rbin},
- "Binwidth distance distribution (nm)" },
- { "-nitacc", FALSE, etBOOL, {&bNitAcc},
- "Regard nitrogen atoms as acceptors" },
- { "-contact", FALSE, etBOOL, {&bContact},
+ { "-r2",
+ FALSE,
+ etREAL,
+ { &r2cut },
+ "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]" },
+ { "-abin", FALSE, etREAL, { &abin }, "Binwidth angle distribution (degrees)" },
+ { "-rbin", FALSE, etREAL, { &rbin }, "Binwidth distance distribution (nm)" },
+ { "-nitacc", FALSE, etBOOL, { &bNitAcc }, "Regard nitrogen atoms as acceptors" },
+ { "-contact",
+ FALSE,
+ etBOOL,
+ { &bContact },
"Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" },
- { "-shell", FALSE, etREAL, {&rshell},
+ { "-shell",
+ FALSE,
+ etREAL,
+ { &rshell },
"when > 0, only calculate hydrogen bonds within # nm shell around "
"one particle" },
- { "-fitstart", FALSE, etREAL, {&fit_start},
- "Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation. With [TT]-gemfit[tt] we suggest [TT]-fitstart 0[tt]" },
- { "-fitend", FALSE, etREAL, {&fit_end},
- "Time (ps) to which to stop fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation (only with [TT]-gemfit[tt])" },
- { "-temp", FALSE, etREAL, {&temp},
- "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
- { "-dump", FALSE, etINT, {&nDump},
+ { "-fitstart",
+ FALSE,
+ etREAL,
+ { &fit_start },
+ "Time (ps) from which to start fitting the correlation functions in order to obtain the "
+ "forward and backward rate constants for HB breaking and formation. With [TT]-gemfit[tt] "
+ "we suggest [TT]-fitstart 0[tt]" },
+ { "-fitend",
+ FALSE,
+ etREAL,
+ { &fit_end },
+ "Time (ps) to which to stop fitting the correlation functions in order to obtain the "
+ "forward and backward rate constants for HB breaking and formation (only with "
+ "[TT]-gemfit[tt])" },
+ { "-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 [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},
- "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
+ { "-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 },
+ "H-bonds between the same donor and acceptor, but with different hydrogen are treated as "
+ "a single H-bond. Mainly important for the ACF." },
#if GMX_OPENMP
- { "-nthreads", FALSE, etINT, {&nThreads},
- "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
+ { "-nthreads",
+ FALSE,
+ etINT,
+ { &nThreads },
+ "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means "
+ "maximum number of threads. Requires linking with OpenMP. The number of threads is "
+ "limited by the number of cores (before OpenMP v.3 ) or environment variable "
+ "OMP_THREAD_LIMIT (OpenMP v.3)" },
#endif
};
- const char *bugs[] = {
- "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being."
+ const char* bugs[] = {
+ "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and "
+ "therefore not available for the time being."
};
- t_filenm fnm[] = {
- { efTRX, "-f", nullptr, ffREAD },
- { efTPR, nullptr, nullptr, ffREAD },
- { efNDX, nullptr, nullptr, ffOPTRD },
- /* { efNDX, "-sel", "select", ffOPTRD },*/
- { efXVG, "-num", "hbnum", ffWRITE },
- { efLOG, "-g", "hbond", ffOPTWR },
- { efXVG, "-ac", "hbac", ffOPTWR },
- { efXVG, "-dist", "hbdist", ffOPTWR },
- { efXVG, "-ang", "hbang", ffOPTWR },
- { efXVG, "-hx", "hbhelix", ffOPTWR },
- { efNDX, "-hbn", "hbond", ffOPTWR },
- { efXPM, "-hbm", "hbmap", ffOPTWR },
- { efXVG, "-don", "donor", ffOPTWR },
- { efXVG, "-dan", "danum", ffOPTWR },
- { efXVG, "-life", "hblife", ffOPTWR },
- { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
+ t_filenm fnm[] = { { efTRX, "-f", nullptr, ffREAD },
+ { efTPR, nullptr, nullptr, ffREAD },
+ { efNDX, nullptr, nullptr, ffOPTRD },
+ /* { efNDX, "-sel", "select", ffOPTRD },*/
+ { efXVG, "-num", "hbnum", ffWRITE },
+ { efLOG, "-g", "hbond", ffOPTWR },
+ { efXVG, "-ac", "hbac", ffOPTWR },
+ { efXVG, "-dist", "hbdist", ffOPTWR },
+ { efXVG, "-ang", "hbang", ffOPTWR },
+ { efXVG, "-hx", "hbhelix", ffOPTWR },
+ { efNDX, "-hbn", "hbond", ffOPTWR },
+ { efXPM, "-hbm", "hbmap", ffOPTWR },
+ { efXVG, "-don", "donor", ffOPTWR },
+ { efXVG, "-dan", "danum", ffOPTWR },
+ { efXVG, "-life", "hblife", ffOPTWR },
+ { efXVG, "-nhbdist", "nhbdist", ffOPTWR }
};
#define NFILE asize(fnm)
- char hbmap [HB_NR] = { ' ', 'o', '-', '*' };
- const char *hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
- t_rgb hbrgb [HB_NR] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} };
-
- t_trxstatus *status;
- bool trrStatus = true;
- t_topology top;
- t_pargs *ppa;
- int npargs, natoms, nframes = 0, shatom;
- int *isize;
- char **grpnames;
- int **index;
- rvec *x, hbox;
- matrix box;
- real t, ccut, dist = 0.0, ang = 0.0;
- double max_nhb, aver_nhb, aver_dist;
- int h = 0, i = 0, j, k = 0, ogrp, nsel;
- int xi = 0, yi, zi, ai;
- int xj, yj, zj, aj, xjj, yjj, zjj;
- gmx_bool bSelected, bHBmap, bStop, bTwo, bBox, bTric;
- int *adist, *rdist;
- int grp, nabin, nrbin, resdist, ihb;
- char **leg;
- t_hbdata *hb;
- FILE *fp, *fpnhb = nullptr, *donor_properties = nullptr;
- t_gridcell ***grid;
- t_ncell *icell, *jcell;
- ivec ngrid;
- unsigned char *datable;
- gmx_output_env_t *oenv;
- int ii, hh, actual_nThreads;
- int threadNr = 0;
- gmx_bool bParallel;
- gmx_bool bEdge_yjj, bEdge_xjj;
-
- t_hbdata **p_hb = nullptr; /* one per thread, then merge after the frame loop */
- int **p_adist = nullptr, **p_rdist = nullptr; /* a histogram for each thread. */
-
- const bool bOMP = GMX_OPENMP;
+ char hbmap[HB_NR] = { ' ', 'o', '-', '*' };
+ const char* hbdesc[HB_NR] = { "None", "Present", "Inserted", "Present & Inserted" };
+ t_rgb hbrgb[HB_NR] = { { 1, 1, 1 }, { 1, 0, 0 }, { 0, 0, 1 }, { 1, 0, 1 } };
+
+ t_trxstatus* status;
+ bool trrStatus = true;
+ t_topology top;
+ t_pargs* ppa;
+ int npargs, natoms, nframes = 0, shatom;
+ int* isize;
+ char** grpnames;
+ int** index;
+ rvec * x, hbox;
+ matrix box;
+ real t, ccut, dist = 0.0, ang = 0.0;
+ double max_nhb, aver_nhb, aver_dist;
+ int h = 0, i = 0, j, k = 0, ogrp, nsel;
+ int xi = 0, yi, zi, ai;
+ int xj, yj, zj, aj, xjj, yjj, zjj;
+ gmx_bool bSelected, bHBmap, bStop, bTwo, bBox, bTric;
+ int * adist, *rdist;
+ int grp, nabin, nrbin, resdist, ihb;
+ char** leg;
+ t_hbdata* hb;
+ FILE * fp, *fpnhb = nullptr, *donor_properties = nullptr;
+ t_gridcell*** grid;
+ t_ncell * icell, *jcell;
+ ivec ngrid;
+ unsigned char* datable;
+ gmx_output_env_t* oenv;
+ int ii, hh, actual_nThreads;
+ int threadNr = 0;
+ gmx_bool bParallel;
+ gmx_bool bEdge_yjj, bEdge_xjj;
+
+ t_hbdata** p_hb = nullptr; /* one per thread, then merge after the frame loop */
+ int ** p_adist = nullptr, **p_rdist = nullptr; /* a histogram for each thread. */
+
+ const bool bOMP = GMX_OPENMP;
npargs = asize(pa);
ppa = add_acf_pargs(&npargs, pa);
- if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
- ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
+ if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs, ppa,
+ asize(desc), desc, asize(bugs), bugs, &oenv))
{
sfree(ppa);
return 0;
/* process input */
bSelected = FALSE;
- ccut = std::cos(acut*DEG2RAD);
+ ccut = std::cos(acut * DEG2RAD);
if (bContact)
{
}
/* Initiate main data structure! */
- bHBmap = (opt2bSet("-ac", NFILE, fnm) ||
- opt2bSet("-life", NFILE, fnm) ||
- opt2bSet("-hbn", NFILE, fnm) ||
- opt2bSet("-hbm", NFILE, fnm));
+ bHBmap = (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)
+ || opt2bSet("-hbn", NFILE, fnm) || opt2bSet("-hbm", NFILE, fnm));
if (opt2bSet("-nhbdist", NFILE, fnm))
{
- const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
- fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm),
- "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv);
+ const char* leg[MAXHH + 1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
+ fpnhb = xvgropen(opt2fn("-nhbdist", NFILE, fnm), "Number of donor-H with N HBs",
+ output_env_get_xvgr_tlabel(oenv), "N", oenv);
xvgr_legend(fpnhb, asize(leg), leg, oenv);
}
hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact);
/* get topology */
- t_inputrec irInstance;
- t_inputrec *ir = &irInstance;
+ t_inputrec irInstance;
+ t_inputrec* ir = &irInstance;
read_tpx_top(ftp2fn(efTPR, NFILE, fnm), ir, box, &natoms, nullptr, nullptr, &top);
snew(grpnames, grNR);
{
/* analyze selected hydrogen bonds */
printf("Select group with selected atoms:\n");
- get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm),
- 1, &nsel, index, grpnames);
+ get_index(&(top.atoms), opt2fn("-sel", NFILE, fnm), 1, &nsel, index, grpnames);
if (nsel % 3)
{
- gmx_fatal(FARGS, "Number of atoms in group '%s' not a multiple of 3\n"
+ gmx_fatal(FARGS,
+ "Number of atoms in group '%s' not a multiple of 3\n"
"and therefore cannot contain triplets of "
- "Donor-Hydrogen-Acceptor", grpnames[0]);
+ "Donor-Hydrogen-Acceptor",
+ grpnames[0]);
}
bTwo = FALSE;
for (i = 0; (i < nsel); i += 3)
{
- int dd = index[0][i];
- int aa = index[0][i+2];
- /* int */ hh = index[0][i+1];
- add_dh (&hb->d, dd, hh, i, datable);
+ int dd = index[0][i];
+ int aa = index[0][i + 2];
+ /* int */ hh = index[0][i + 1];
+ add_dh(&hb->d, dd, hh, i, datable);
add_acc(&hb->a, aa, i);
/* Should this be here ? */
snew(hb->d.dptr, top.atoms.nr);
snew(hb->a.aptr, top.atoms.nr);
add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact);
}
- printf("Analyzing %d selected hydrogen bonds from '%s'\n",
- isize[0], grpnames[0]);
+ printf("Analyzing %d selected hydrogen bonds from '%s'\n", isize[0], grpnames[0]);
}
else
{
/* analyze all hydrogen bonds: get group(s) */
printf("Specify 2 groups to analyze:\n");
- get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
- 2, isize, index, grpnames);
+ get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpnames);
/* check if we have two identical or two non-overlapping groups */
bTwo = isize[0] != isize[1];
}
if (bTwo)
{
- printf("Checking for overlap in atoms between %s and %s\n",
- grpnames[0], grpnames[1]);
+ printf("Checking for overlap in atoms between %s and %s\n", grpnames[0], grpnames[1]);
gen_datable(index[0], isize[0], datable, top.atoms.nr);
{
if (ISINGRP(datable[index[1][i]]))
{
- gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'",
- grpnames[0], grpnames[1]);
+ gmx_fatal(FARGS, "Partial overlap between groups '%s' and '%s'", grpnames[0],
+ grpnames[1]);
}
}
}
{
printf("Calculating %s "
"between %s (%d atoms) and %s (%d atoms)\n",
- bContact ? "contacts" : "hydrogen bonds",
- grpnames[0], isize[0], grpnames[1], isize[1]);
+ bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0], grpnames[1],
+ isize[1]);
}
else
{
snew(datable, top.atoms.nr);
for (i = 0; (i < grNR); i++)
{
- if ( ((i == gr0) && !bSelected ) ||
- ((i == gr1) && bTwo ))
+ if (((i == gr0) && !bSelected) || ((i == gr1) && bTwo))
{
gen_datable(index[i], isize[i], datable, top.atoms.nr);
if (bContact)
{
- search_acceptors(&top, isize[i], index[i], &hb->a, i,
- bNitAcc, TRUE, (bTwo && (i == gr0)) || !bTwo, datable);
- search_donors (&top, isize[i], index[i], &hb->d, i,
- TRUE, (bTwo && (i == gr1)) || !bTwo, datable);
+ search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, TRUE,
+ (bTwo && (i == gr0)) || !bTwo, datable);
+ search_donors(&top, isize[i], index[i], &hb->d, i, TRUE,
+ (bTwo && (i == gr1)) || !bTwo, datable);
}
else
{
search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE, TRUE, datable);
- search_donors (&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
+ search_donors(&top, isize[i], index[i], &hb->d, i, FALSE, TRUE, datable);
}
if (bTwo)
{
shatom = 0;
if (rshell > 0)
{
- int shisz;
- int *shidx;
- char *shgrpnm;
+ int shisz;
+ int* shidx;
+ char* shgrpnm;
/* get index group with atom for shell */
do
{
printf("Select atom for shell (1 atom):\n");
- get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
- 1, &shisz, &shidx, &shgrpnm);
+ get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &shisz, &shidx, &shgrpnm);
if (shisz != 1)
{
printf("group contains %d atoms, should be 1 (one)\n", shisz);
}
- }
- while (shisz != 1);
+ } while (shisz != 1);
shatom = shidx[0];
printf("Will calculate hydrogen bonds within a shell "
- "of %g nm around atom %i\n", rshell, shatom+1);
+ "of %g nm around atom %i\n",
+ rshell, shatom + 1);
}
/* Analyze trajectory */
natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
if (natoms > top.atoms.nr)
{
- gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
- top.atoms.nr, natoms);
+ gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)", top.atoms.nr, natoms);
}
bBox = (ir->ePBC != epbcNONE);
grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid);
- nabin = static_cast<int>(acut/abin);
- nrbin = static_cast<int>(rcut/rbin);
- snew(adist, nabin+1);
- snew(rdist, nrbin+1);
+ nabin = static_cast<int>(acut / abin);
+ nrbin = static_cast<int>(rcut / rbin);
+ snew(adist, nabin + 1);
+ snew(rdist, nrbin + 1);
#if !GMX_OPENMP
-#define __ADIST adist
-#define __RDIST rdist
-#define __HBDATA hb
-#else /* GMX_OPENMP ================================================== \
- * Set up the OpenMP stuff, |
- * like the number of threads and such |
- * Also start the parallel loop. |
+# define __ADIST adist
+# define __RDIST rdist
+# define __HBDATA hb
+#else /* GMX_OPENMP ================================================== \
+ * Set up the OpenMP stuff, | \
+ * like the number of threads and such | \
+ * Also start the parallel loop. | \
*/
-#define __ADIST p_adist[threadNr]
-#define __RDIST p_rdist[threadNr]
-#define __HBDATA p_hb[threadNr]
+# define __ADIST p_adist[threadNr]
+# define __RDIST p_rdist[threadNr]
+# define __HBDATA p_hb[threadNr]
#endif
if (bOMP)
{
actual_nThreads = 1;
}
- snew(p_hb, actual_nThreads);
+ snew(p_hb, actual_nThreads);
snew(p_adist, actual_nThreads);
snew(p_rdist, actual_nThreads);
for (i = 0; i < actual_nThreads; i++)
{
snew(p_hb[i], 1);
- snew(p_adist[i], nabin+1);
- snew(p_rdist[i], nrbin+1);
+ snew(p_adist[i], nabin + 1);
+ snew(p_rdist[i], nrbin + 1);
p_hb[i]->max_frames = 0;
p_hb[i]->nhb = nullptr;
p_hb[i]->time = nullptr;
p_hb[i]->nhx = nullptr;
- p_hb[i]->bHBmap = hb->bHBmap;
- p_hb[i]->bDAnr = hb->bDAnr;
- p_hb[i]->wordlen = hb->wordlen;
- p_hb[i]->nframes = hb->nframes;
- p_hb[i]->maxhydro = hb->maxhydro;
- p_hb[i]->danr = hb->danr;
- p_hb[i]->d = hb->d;
- p_hb[i]->a = hb->a;
- p_hb[i]->hbmap = hb->hbmap;
- p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
+ p_hb[i]->bHBmap = hb->bHBmap;
+ p_hb[i]->bDAnr = hb->bDAnr;
+ p_hb[i]->wordlen = hb->wordlen;
+ p_hb[i]->nframes = hb->nframes;
+ p_hb[i]->maxhydro = hb->maxhydro;
+ p_hb[i]->danr = hb->danr;
+ p_hb[i]->d = hb->d;
+ p_hb[i]->a = hb->a;
+ p_hb[i]->hbmap = hb->hbmap;
+ p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
p_hb[i]->nrhb = 0;
p_hb[i]->nrdist = 0;
/* Make a thread pool here,
* instead of forking anew at every frame. */
-#pragma omp parallel \
- firstprivate(i) \
- private(j, h, ii, hh, \
- xi, yi, zi, xj, yj, zj, threadNr, \
- dist, ang, icell, jcell, \
- grp, ogrp, ai, aj, xjj, yjj, zjj, \
- ihb, resdist, \
- k, bTric, \
- bEdge_xjj, bEdge_yjj) \
- default(shared)
- { /* Start of parallel region */
+#pragma omp parallel firstprivate(i) private( \
+ j, h, ii, hh, xi, yi, zi, xj, yj, zj, threadNr, dist, ang, icell, jcell, grp, ogrp, ai, \
+ aj, xjj, yjj, zjj, ihb, resdist, k, bTric, bEdge_xjj, bEdge_yjj) default(shared)
+ { /* Start of parallel region */
if (bOMP)
{
threadNr = gmx_omp_get_thread_num();
{
sync_hbdata(p_hb[threadNr], nframes);
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
#pragma omp single
{
count_da_grid(ngrid, grid, hb->danr[nframes]);
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
} /* omp single */
if (bOMP)
/* int ii; */
for (ii = 0; (ii < nsel); ii++)
{
- int dd = index[0][i];
- int aa = index[0][i+2];
- /* int */ hh = index[0][i+1];
- ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
- hbox, &dist, &ang, bDA, &h, bContact, bMerge);
+ int dd = index[0][i];
+ int aa = index[0][i + 2];
+ /* int */ hh = index[0][i + 1];
+ ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box,
+ hbox, &dist, &ang, bDA, &h, bContact, bMerge);
if (ihb)
{
}
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
} /* omp single */
} /* if (bSelected) */
else
if (bTwo)
{
- ogrp = 1-grp;
+ ogrp = 1 - grp;
}
else
{
*/
for (ai = 0; (ai < icell->nr); ai++)
{
- i = icell->atoms[ai];
+ i = icell->atoms[ai];
/* loop over all adjacent gridcells (xj,yj,zj) */
for (zjj = grid_loop_begin(ngrid[ZZ], zi, bTric, FALSE);
- zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE);
- zjj++)
+ zjj <= grid_loop_end(ngrid[ZZ], zi, bTric, FALSE); zjj++)
{
zj = grid_mod(zjj, ngrid[ZZ]);
bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
yjj++)
{
yj = grid_mod(yjj, ngrid[YY]);
- bEdge_xjj =
- (yj == 0) || (yj == ngrid[YY] - 1) ||
- (zj == 0) || (zj == ngrid[ZZ] - 1);
+ bEdge_xjj = (yj == 0) || (yj == ngrid[YY] - 1)
+ || (zj == 0) || (zj == ngrid[ZZ] - 1);
for (xjj = grid_loop_begin(ngrid[XX], xi, bTric, bEdge_xjj);
xjj <= grid_loop_end(ngrid[XX], xi, bTric, bEdge_xjj);
xjj++)
{
xj = grid_mod(xjj, ngrid[XX]);
jcell = &(grid[zj][yj][xj].a[ogrp]);
- /* loop over acceptor atoms from other group (ogrp)
- * in this adjacent gridcell (jcell)
+ /* loop over acceptor atoms from other group
+ * (ogrp) in this adjacent gridcell (jcell)
*/
for (aj = 0; (aj < jcell->nr); aj++)
{
j = jcell->atoms[aj];
/* check if this once was a h-bond */
- ihb = is_hbond(__HBDATA, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box,
- hbox, &dist, &ang, bDA, &h, bContact, bMerge);
+ ihb = is_hbond(__HBDATA, grp, ogrp, i, j,
+ rcut, r2cut, ccut, x, bBox,
+ box, hbox, &dist, &ang, bDA,
+ &h, bContact, bMerge);
if (ihb)
{
/* add to index if not already there */
/* Add a hbond */
- add_hbond(__HBDATA, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact);
+ add_hbond(__HBDATA, i, j, h, grp, ogrp,
+ nframes, bMerge, ihb, bContact);
/* make angle and distance distributions */
if (ihb == hbHB && !bContact)
{
if (dist > rcut)
{
- gmx_fatal(FARGS, "distance is higher than what is allowed for an hbond: %f", dist);
+ gmx_fatal(
+ FARGS,
+ "distance is higher "
+ "than what is allowed "
+ "for an hbond: %f",
+ dist);
}
ang *= RAD2DEG;
- __ADIST[static_cast<int>( ang/abin)]++;
- __RDIST[static_cast<int>(dist/rbin)]++;
+ __ADIST[static_cast<int>(ang / abin)]++;
+ __RDIST[static_cast<int>(dist / rbin)]++;
if (!bTwo)
{
if (donor_index(&hb->d, grp, i) == NOTSET)
{
- gmx_fatal(FARGS, "Invalid donor %d", i);
+ gmx_fatal(
+ FARGS,
+ "Invalid donor %d", i);
}
- if (acceptor_index(&hb->a, ogrp, j) == NOTSET)
+ if (acceptor_index(&hb->a, ogrp, j)
+ == NOTSET)
{
- gmx_fatal(FARGS, "Invalid acceptor %d", j);
+ gmx_fatal(FARGS,
+ "Invalid "
+ "acceptor %d",
+ j);
}
- resdist = std::abs(top.atoms.atom[i].resind-top.atoms.atom[j].resind);
+ resdist = std::abs(
+ top.atoms.atom[i].resind
+ - top.atoms.atom[j].resind);
if (resdist >= max_hx)
{
- resdist = max_hx-1;
+ resdist = max_hx - 1;
}
__HBDATA->nhx[nframes][resdist]++;
}
}
-
}
} /* for aj */
} /* for xjj */
} /* for xi,yi,zi */
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
} /* if (bSelected) {...} else */
/* Sum up histograms and counts from p_hb[] into hb */
if (bOMP)
{
- hb->nhb[k] += p_hb[threadNr]->nhb[k];
+ hb->nhb[k] += p_hb[threadNr]->nhb[k];
hb->ndist[k] += p_hb[threadNr]->ndist[k];
for (j = 0; j < max_hx; j++)
{
- hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
+ hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
}
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
/* Here are a handful of single constructs
{
analyse_donor_properties(donor_properties, hb, k, t);
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
#pragma omp single
do_nhb_dist(fpnhb, hb, t);
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
#pragma omp single
trrStatus = (read_next_x(oenv, status, &t, x, box));
nframes++;
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
#pragma omp barrier
- }
- while (trrStatus);
+ } while (trrStatus);
if (bOMP)
{
#pragma omp critical
{
- hb->nrhb += p_hb[threadNr]->nrhb;
+ hb->nrhb += p_hb[threadNr]->nrhb;
hb->nrdist += p_hb[threadNr]->nrdist;
}
adist[i] += p_adist[j][i];
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
#pragma omp for
for (i = 0; i <= nrbin; i++)
rdist[i] += p_rdist[j][i];
}
}
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
sfree(p_adist[threadNr]);
if (nframes < 2 && (opt2bSet("-ac", NFILE, fnm) || opt2bSet("-life", NFILE, fnm)))
{
- gmx_fatal(FARGS, "Cannot calculate autocorrelation of life times with less than two frames");
+ gmx_fatal(FARGS,
+ "Cannot calculate autocorrelation of life times with less than two frames");
}
free_grid(ngrid, &grid);
}
else
{
- max_nhb = 0.5*(hb->d.nrd*hb->a.nra);
+ max_nhb = 0.5 * (hb->d.nrd * hb->a.nra);
}
if (bHBmap)
{
printf("Found %d different %s in trajectory\n"
"Found %d different atom-pairs within %s distance\n",
- hb->nrhb, bContact ? "contacts" : "hydrogen bonds",
- hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
+ hb->nrhb, bContact ? "contacts" : "hydrogen bonds", hb->nrdist,
+ (r2cut > 0) ? "second cut-off" : "hydrogen bonding");
if (bMerge)
{
/* Print out number of hbonds and distances */
aver_nhb = 0;
aver_dist = 0;
- fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" :
- "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv);
+ fp = xvgropen(opt2fn("-num", NFILE, fnm), bContact ? "Contacts" : "Hydrogen Bonds",
+ output_env_get_xvgr_tlabel(oenv), "Number", oenv);
snew(leg, 2);
snew(leg[0], STRLEN);
snew(leg[1], STRLEN);
for (i = 0; (i < nframes); i++)
{
fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]);
- aver_nhb += hb->nhb[i];
+ aver_nhb += hb->nhb[i];
aver_dist += hb->ndist[i];
}
xvgrclose(fp);
- aver_nhb /= nframes;
+ aver_nhb /= nframes;
aver_dist /= nframes;
/* Print HB distance distribution */
if (opt2bSet("-dist", NFILE, fnm))
sum += rdist[i];
}
- fp = xvgropen(opt2fn("-dist", NFILE, fnm),
- "Hydrogen Bond Distribution",
- bDA ?
- "Donor - Acceptor Distance (nm)" :
- "Hydrogen - Acceptor Distance (nm)", "", oenv);
+ fp = xvgropen(opt2fn("-dist", NFILE, fnm), "Hydrogen Bond Distribution",
+ bDA ? "Donor - Acceptor Distance (nm)" : "Hydrogen - Acceptor Distance (nm)",
+ "", oenv);
for (i = 0; i < nrbin; i++)
{
- fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*sum));
+ fprintf(fp, "%10g %10g\n", (i + 0.5) * rbin, rdist[i] / (rbin * sum));
}
xvgrclose(fp);
}
sum += adist[i];
}
- fp = xvgropen(opt2fn("-ang", NFILE, fnm),
- "Hydrogen Bond Distribution",
+ fp = xvgropen(opt2fn("-ang", NFILE, fnm), "Hydrogen Bond Distribution",
"Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv);
for (i = 0; i < nabin; i++)
{
- fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*sum));
+ fprintf(fp, "%10g %10g\n", (i + 0.5) * abin, adist[i] / (abin * sum));
}
xvgrclose(fp);
}
/* Print HB in alpha-helix */
if (opt2bSet("-hx", NFILE, fnm))
{
- fp = xvgropen(opt2fn("-hx", NFILE, fnm),
- "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
+ fp = xvgropen(opt2fn("-hx", NFILE, fnm), "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv),
+ "Count", oenv);
xvgr_legend(fp, NRHXTYPES, hxtypenames, oenv);
for (i = 0; i < nframes; i++)
{
}
printf("Average number of %s per timeframe %.3f out of %g possible\n",
- bContact ? "contacts" : "hbonds",
- bContact ? aver_dist : aver_nhb, max_nhb);
+ bContact ? "contacts" : "hbonds", bContact ? aver_dist : aver_nhb, max_nhb);
/* Do Autocorrelation etc. */
if (hb->bHBmap)
}
if (opt2bSet("-ac", NFILE, fnm))
{
- do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
- bMerge, bContact, fit_start, temp, r2cut > 0, oenv,
- nThreads);
+ do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump, bMerge, bContact, fit_start, temp,
+ r2cut > 0, oenv, nThreads);
}
if (opt2bSet("-life", NFILE, fnm))
{
mat.ny = hb->nrhb;
mat.matrix.resize(mat.nx, mat.ny);
- for (auto &value : mat.matrix.toArrayRef())
+ for (auto& value : mat.matrix.toArrayRef())
{
value = 0;
}
{
int nn0 = hb->hbmap[id][ia]->n0;
range_check(y, 0, mat.ny);
- mat.matrix(x+nn0, y) = static_cast<t_matelmt>(is_hb(hb->hbmap[id][ia]->h[hh], x));
+ mat.matrix(x + nn0, y) = static_cast<t_matelmt>(
+ is_hb(hb->hbmap[id][ia]->h[hh], x));
}
y++;
}
std::copy(hb->time, hb->time + mat.nx, mat.axis_x.begin());
mat.axis_y.resize(mat.ny);
std::iota(mat.axis_y.begin(), mat.axis_y.end(), 0);
- mat.title = (bContact ? "Contact Existence Map" :
- "Hydrogen Bond Existence Map");
- mat.legend = bContact ? "Contacts" : "Hydrogen Bonds";
- mat.label_x = output_env_get_xvgr_tlabel(oenv);
- mat.label_y = bContact ? "Contact Index" : "Hydrogen Bond Index";
+ mat.title = (bContact ? "Contact Existence Map" : "Hydrogen Bond Existence Map");
+ mat.legend = bContact ? "Contacts" : "Hydrogen Bonds";
+ mat.label_x = output_env_get_xvgr_tlabel(oenv);
+ mat.label_y = bContact ? "Contact Index" : "Hydrogen Bond Index";
mat.bDiscrete = true;
mat.map.resize(2);
- for (auto &m : mat.map)
+ for (auto& m : mat.map)
{
m.code.c1 = hbmap[i];
m.desc = hbdesc[i];
}
else
{
- fprintf(stderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n");
+ fprintf(stderr,
+ "No hydrogen bonds/contacts found. No hydrogen bond map will be "
+ "printed.\n");
}
}
}
if (hb->bDAnr)
{
int i, j, nleg;
- char **legnames;
+ char** legnames;
char buf[STRLEN];
-#define USE_THIS_GROUP(j) ( ((j) == gr0) || (bTwo && ((j) == gr1)) )
+#define USE_THIS_GROUP(j) (((j) == gr0) || (bTwo && ((j) == gr1)))
- fp = xvgropen(opt2fn("-dan", NFILE, fnm),
- "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv);
- nleg = (bTwo ? 2 : 1)*2;
+ fp = xvgropen(opt2fn("-dan", NFILE, fnm), "Donors and Acceptors",
+ output_env_get_xvgr_tlabel(oenv), "Count", oenv);
+ nleg = (bTwo ? 2 : 1) * 2;
snew(legnames, nleg);
i = 0;
for (j = 0; j < grNR; j++)
{
- if (USE_THIS_GROUP(j) )
+ if (USE_THIS_GROUP(j))
{
sprintf(buf, "Donors %s", grpnames[j]);
legnames[i++] = gmx_strdup(buf);
fprintf(fp, "%10g", hb->time[i]);
for (j = 0; (j < grNR); j++)
{
- if (USE_THIS_GROUP(j) )
+ if (USE_THIS_GROUP(j))
{
fprintf(fp, " %6d", hb->danr[i][j]);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff