| File: | gromacs/gmxana/gmx_hbond.c |
| Location: | line 2686, column 5 |
| Description: | Value stored to 'anhb' is never read |
| 1 | /* |
| 2 | * This file is part of the GROMACS molecular simulation package. |
| 3 | * |
| 4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
| 5 | * Copyright (c) 2001-2008, The GROMACS development team. |
| 6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
| 7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
| 8 | * and including many others, as listed in the AUTHORS file in the |
| 9 | * top-level source directory and at http://www.gromacs.org. |
| 10 | * |
| 11 | * GROMACS is free software; you can redistribute it and/or |
| 12 | * modify it under the terms of the GNU Lesser General Public License |
| 13 | * as published by the Free Software Foundation; either version 2.1 |
| 14 | * of the License, or (at your option) any later version. |
| 15 | * |
| 16 | * GROMACS is distributed in the hope that it will be useful, |
| 17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 19 | * Lesser General Public License for more details. |
| 20 | * |
| 21 | * You should have received a copy of the GNU Lesser General Public |
| 22 | * License along with GROMACS; if not, see |
| 23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
| 24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
| 25 | * |
| 26 | * If you want to redistribute modifications to GROMACS, please |
| 27 | * consider that scientific software is very special. Version |
| 28 | * control is crucial - bugs must be traceable. We will be happy to |
| 29 | * consider code for inclusion in the official distribution, but |
| 30 | * derived work must not be called official GROMACS. Details are found |
| 31 | * in the README & COPYING files - if they are missing, get the |
| 32 | * official version at http://www.gromacs.org. |
| 33 | * |
| 34 | * To help us fund GROMACS development, we humbly ask that you cite |
| 35 | * the research papers on the package. Check out http://www.gromacs.org. |
| 36 | */ |
| 37 | #ifdef HAVE_CONFIG_H1 |
| 38 | #include <config.h> |
| 39 | #endif |
| 40 | #include <math.h> |
| 41 | |
| 42 | /*#define HAVE_NN_LOOPS*/ |
| 43 | |
| 44 | #include "gromacs/commandline/pargs.h" |
| 45 | #include "copyrite.h" |
| 46 | #include "txtdump.h" |
| 47 | #include "physics.h" |
| 48 | #include "macros.h" |
| 49 | #include "gromacs/utility/fatalerror.h" |
| 50 | #include "index.h" |
| 51 | #include "gromacs/utility/smalloc.h" |
| 52 | #include "gromacs/math/vec.h" |
| 53 | #include "gromacs/fileio/xvgr.h" |
| 54 | #include "viewit.h" |
| 55 | #include "gstat.h" |
| 56 | #include "gromacs/utility/cstringutil.h" |
| 57 | #include "pbc.h" |
| 58 | #include "correl.h" |
| 59 | #include "gmx_ana.h" |
| 60 | #include "geminate.h" |
| 61 | |
| 62 | #include "gromacs/utility/futil.h" |
| 63 | #include "gromacs/fileio/matio.h" |
| 64 | #include "gromacs/fileio/tpxio.h" |
| 65 | #include "gromacs/fileio/trxio.h" |
| 66 | #include "gromacs/utility/gmxomp.h" |
| 67 | |
| 68 | typedef short int t_E; |
| 69 | typedef int t_EEst; |
| 70 | #define max_hx7 7 |
| 71 | typedef int t_hx[max_hx7]; |
| 72 | #define NRHXTYPES7 max_hx7 |
| 73 | const char *hxtypenames[NRHXTYPES7] = |
| 74 | {"n-n", "n-n+1", "n-n+2", "n-n+3", "n-n+4", "n-n+5", "n-n>6"}; |
| 75 | #define MAXHH4 4 |
| 76 | |
| 77 | #ifdef GMX_OPENMP |
| 78 | #define MASTER_THREAD_ONLY(threadNr)((threadNr) == (threadNr)) ((threadNr) == 0) |
| 79 | #else |
| 80 | #define MASTER_THREAD_ONLY(threadNr)((threadNr) == (threadNr)) ((threadNr) == (threadNr)) |
| 81 | #endif |
| 82 | |
| 83 | /* -----------------------------------------*/ |
| 84 | |
| 85 | enum { |
| 86 | gr0, gr1, grI, grNR |
| 87 | }; |
| 88 | enum { |
| 89 | hbNo, hbDist, hbHB, hbNR, hbR2 |
| 90 | }; |
| 91 | enum { |
| 92 | noDA, ACC, DON, DA, INGROUP |
| 93 | }; |
| 94 | enum { |
| 95 | NN_NULL, NN_NONE, NN_BINARY, NN_1_over_r3, NN_dipole, NN_NR |
| 96 | }; |
| 97 | |
| 98 | static const char *grpnames[grNR] = {"0", "1", "I" }; |
| 99 | |
| 100 | static gmx_bool bDebug = FALSE0; |
| 101 | |
| 102 | #define HB_NO0 0 |
| 103 | #define HB_YES1<<0 1<<0 |
| 104 | #define HB_INS1<<1 1<<1 |
| 105 | #define HB_YESINS1<<0|1<<1 HB_YES1<<0|HB_INS1<<1 |
| 106 | #define HB_NR(1<<2) (1<<2) |
| 107 | #define MAXHYDRO4 4 |
| 108 | |
| 109 | #define ISHB(h)(((h) & 2) == 2) (((h) & 2) == 2) |
| 110 | #define ISDIST(h)(((h) & 1) == 1) (((h) & 1) == 1) |
| 111 | #define ISDIST2(h)(((h) & 4) == 4) (((h) & 4) == 4) |
| 112 | #define ISACC(h)(((h) & 1) == 1) (((h) & 1) == 1) |
| 113 | #define ISDON(h)(((h) & 2) == 2) (((h) & 2) == 2) |
| 114 | #define ISINGRP(h)(((h) & 4) == 4) (((h) & 4) == 4) |
| 115 | |
| 116 | typedef struct { |
| 117 | int nr; |
| 118 | int maxnr; |
| 119 | atom_id *atoms; |
| 120 | } t_ncell; |
| 121 | |
| 122 | typedef struct { |
| 123 | t_ncell d[grNR]; |
| 124 | t_ncell a[grNR]; |
| 125 | } t_gridcell; |
| 126 | |
| 127 | typedef int t_icell[grNR]; |
| 128 | typedef atom_id h_id[MAXHYDRO4]; |
| 129 | |
| 130 | typedef struct { |
| 131 | int history[MAXHYDRO4]; |
| 132 | /* Has this hbond existed ever? If so as hbDist or hbHB or both. |
| 133 | * Result is stored as a bitmap (1 = hbDist) || (2 = hbHB) |
| 134 | */ |
| 135 | /* Bitmask array which tells whether a hbond is present |
| 136 | * at a given time. Either of these may be NULL |
| 137 | */ |
| 138 | int n0; /* First frame a HB was found */ |
| 139 | int nframes, maxframes; /* Amount of frames in this hbond */ |
| 140 | unsigned int **h; |
| 141 | unsigned int **g; |
| 142 | /* See Xu and Berne, JPCB 105 (2001), p. 11929. We define the |
| 143 | * function g(t) = [1-h(t)] H(t) where H(t) is one when the donor- |
| 144 | * acceptor distance is less than the user-specified distance (typically |
| 145 | * 0.35 nm). |
| 146 | */ |
| 147 | } t_hbond; |
| 148 | |
| 149 | typedef struct { |
| 150 | int nra, max_nra; |
| 151 | atom_id *acc; /* Atom numbers of the acceptors */ |
| 152 | int *grp; /* Group index */ |
| 153 | int *aptr; /* Map atom number to acceptor index */ |
| 154 | } t_acceptors; |
| 155 | |
| 156 | typedef struct { |
| 157 | int nrd, max_nrd; |
| 158 | int *don; /* Atom numbers of the donors */ |
| 159 | int *grp; /* Group index */ |
| 160 | int *dptr; /* Map atom number to donor index */ |
| 161 | int *nhydro; /* Number of hydrogens for each donor */ |
| 162 | h_id *hydro; /* The atom numbers of the hydrogens */ |
| 163 | h_id *nhbonds; /* The number of HBs per H at current */ |
| 164 | } t_donors; |
| 165 | |
| 166 | /* Tune this to match memory requirements. It should be a signed integer type, e.g. signed char.*/ |
| 167 | #define PSTYPEint int |
| 168 | |
| 169 | typedef struct { |
| 170 | int len; /* The length of frame and p. */ |
| 171 | int *frame; /* The frames at which transitio*/ |
| 172 | PSTYPEint *p; |
| 173 | } t_pShift; |
| 174 | |
| 175 | typedef struct { |
| 176 | /* Periodicity history. Used for the reversible geminate recombination. */ |
| 177 | t_pShift **pHist; /* The periodicity of every hbond in t_hbdata->hbmap: |
| 178 | * pHist[d][a]. We can safely assume that the same |
| 179 | * periodic shift holds for all hydrogens of a da-pair. |
| 180 | * |
| 181 | * Nowadays it only stores TRANSITIONS, and not the shift at every frame. |
| 182 | * That saves a LOT of memory, an hopefully kills a mysterious bug where |
| 183 | * pHist gets contaminated. */ |
| 184 | |
| 185 | PSTYPEint nper; /* The length of p2i */ |
| 186 | ivec *p2i; /* Maps integer to periodic shift for a pair.*/ |
| 187 | matrix P; /* Projection matrix to find the box shifts. */ |
| 188 | int gemtype; /* enumerated type */ |
| 189 | } t_gemPeriod; |
| 190 | |
| 191 | typedef struct { |
| 192 | int nframes; |
| 193 | int *Etot; /* Total energy for each frame */ |
| 194 | t_E ****E; /* Energy estimate for [d][a][h][frame-n0] */ |
| 195 | } t_hbEmap; |
| 196 | |
| 197 | typedef struct { |
| 198 | gmx_bool bHBmap, bDAnr, bGem; |
| 199 | int wordlen; |
| 200 | /* The following arrays are nframes long */ |
| 201 | int nframes, max_frames, maxhydro; |
| 202 | int *nhb, *ndist; |
| 203 | h_id *n_bound; |
| 204 | real *time; |
| 205 | t_icell *danr; |
| 206 | t_hx *nhx; |
| 207 | /* These structures are initialized from the topology at start up */ |
| 208 | t_donors d; |
| 209 | t_acceptors a; |
| 210 | /* This holds a matrix with all possible hydrogen bonds */ |
| 211 | int nrhb, nrdist; |
| 212 | t_hbond ***hbmap; |
| 213 | #ifdef HAVE_NN_LOOPS |
| 214 | t_hbEmap hbE; |
| 215 | #endif |
| 216 | /* For parallelization reasons this will have to be a pointer. |
| 217 | * Otherwise discrepancies may arise between the periodicity data |
| 218 | * seen by different threads. */ |
| 219 | t_gemPeriod *per; |
| 220 | } t_hbdata; |
| 221 | |
| 222 | static void clearPshift(t_pShift *pShift) |
| 223 | { |
| 224 | if (pShift->len > 0) |
| 225 | { |
| 226 | sfree(pShift->p)save_free("pShift->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 226, (pShift->p)); |
| 227 | sfree(pShift->frame)save_free("pShift->frame", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 227, (pShift->frame)); |
| 228 | pShift->len = 0; |
| 229 | } |
| 230 | } |
| 231 | |
| 232 | static void calcBoxProjection(matrix B, matrix P) |
| 233 | { |
| 234 | const int vp[] = {XX0, YY1, ZZ2}; |
| 235 | int i, j; |
| 236 | int m, n; |
| 237 | matrix M, N, U; |
| 238 | |
| 239 | for (i = 0; i < 3; i++) |
| 240 | { |
| 241 | m = vp[i]; |
| 242 | for (j = 0; j < 3; j++) |
| 243 | { |
| 244 | n = vp[j]; |
| 245 | U[m][n] = i == j ? 1 : 0; |
| 246 | } |
| 247 | } |
| 248 | m_inv(B, M); |
| 249 | for (i = 0; i < 3; i++) |
| 250 | { |
| 251 | m = vp[i]; |
| 252 | mvmul(M, U[m], P[m]); |
| 253 | } |
| 254 | transpose(P, N); |
| 255 | } |
| 256 | |
| 257 | static void calcBoxDistance(matrix P, rvec d, ivec ibd) |
| 258 | { |
| 259 | /* returns integer distance in box coordinates. |
| 260 | * P is the projection matrix from cartesian coordinates |
| 261 | * obtained with calcBoxProjection(). */ |
| 262 | int i; |
| 263 | rvec bd; |
| 264 | mvmul(P, d, bd); |
| 265 | /* extend it by 0.5 in all directions since (int) rounds toward 0.*/ |
| 266 | for (i = 0; i < 3; i++) |
| 267 | { |
| 268 | bd[i] = bd[i] + (bd[i] < 0 ? -0.5 : 0.5); |
| 269 | } |
| 270 | ibd[XX0] = (int)bd[XX0]; |
| 271 | ibd[YY1] = (int)bd[YY1]; |
| 272 | ibd[ZZ2] = (int)bd[ZZ2]; |
| 273 | } |
| 274 | |
| 275 | /* Changed argument 'bMerge' into 'oneHB' below, |
| 276 | * since -contact should cause maxhydro to be 1, |
| 277 | * not just -merge. |
| 278 | * - Erik Marklund May 29, 2006 |
| 279 | */ |
| 280 | |
| 281 | static PSTYPEint periodicIndex(ivec r, t_gemPeriod *per, gmx_bool daSwap) |
| 282 | { |
| 283 | /* Try to merge hbonds on the fly. That means that if the |
| 284 | * acceptor and donor are mergable, then: |
| 285 | * 1) store the hb-info so that acceptor id > donor id, |
| 286 | * 2) add the periodic shift in pairs, so that [-x,-y,-z] is |
| 287 | * stored in per.p2i[] whenever acceptor id < donor id. |
| 288 | * Note that [0,0,0] should already be the first element of per.p2i |
| 289 | * by the time this function is called. */ |
| 290 | |
| 291 | /* daSwap is TRUE if the donor and acceptor were swapped. |
| 292 | * If so, then the negative vector should be used. */ |
| 293 | PSTYPEint i; |
| 294 | |
| 295 | if (per->p2i == NULL((void*)0) || per->nper == 0) |
| 296 | { |
| 297 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 297, "'per' not initialized properly."); |
| 298 | } |
| 299 | for (i = 0; i < per->nper; i++) |
| 300 | { |
| 301 | if (r[XX0] == per->p2i[i][XX0] && |
| 302 | r[YY1] == per->p2i[i][YY1] && |
| 303 | r[ZZ2] == per->p2i[i][ZZ2]) |
| 304 | { |
| 305 | return i; |
| 306 | } |
| 307 | } |
| 308 | /* Not found apparently. Add it to the list! */ |
| 309 | /* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */ |
| 310 | |
| 311 | #pragma omp critical |
| 312 | { |
| 313 | if (!per->p2i) |
| 314 | { |
| 315 | fprintf(stderrstderr, "p2i not initialized. This shouldn't happen!\n"); |
| 316 | snew(per->p2i, 1)(per->p2i) = save_calloc("per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 316, (1), sizeof(*(per->p2i))); |
| 317 | } |
| 318 | else |
| 319 | { |
| 320 | srenew(per->p2i, per->nper+2)(per->p2i) = save_realloc("per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 320, (per->p2i), (per->nper+2), sizeof(*(per->p2i) )); |
| 321 | } |
| 322 | copy_ivec(r, per->p2i[per->nper]); |
| 323 | (per->nper)++; |
| 324 | |
| 325 | /* Add the mirror too. It's rather likely that it'll be needed. */ |
| 326 | per->p2i[per->nper][XX0] = -r[XX0]; |
| 327 | per->p2i[per->nper][YY1] = -r[YY1]; |
| 328 | per->p2i[per->nper][ZZ2] = -r[ZZ2]; |
| 329 | (per->nper)++; |
| 330 | } /* omp critical */ |
| 331 | return per->nper - 1 - (daSwap ? 0 : 1); |
| 332 | } |
| 333 | |
| 334 | static t_hbdata *mk_hbdata(gmx_bool bHBmap, gmx_bool bDAnr, gmx_bool oneHB, gmx_bool bGem, int gemmode) |
| 335 | { |
| 336 | t_hbdata *hb; |
| 337 | |
| 338 | snew(hb, 1)(hb) = save_calloc("hb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 338, (1), sizeof(*(hb))); |
| 339 | hb->wordlen = 8*sizeof(unsigned int); |
| 340 | hb->bHBmap = bHBmap; |
| 341 | hb->bDAnr = bDAnr; |
| 342 | hb->bGem = bGem; |
| 343 | if (oneHB) |
| 344 | { |
| 345 | hb->maxhydro = 1; |
| 346 | } |
| 347 | else |
| 348 | { |
| 349 | hb->maxhydro = MAXHYDRO4; |
| 350 | } |
| 351 | snew(hb->per, 1)(hb->per) = save_calloc("hb->per", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 351, (1), sizeof(*(hb->per))); |
| 352 | hb->per->gemtype = bGem ? gemmode : 0; |
| 353 | |
| 354 | return hb; |
| 355 | } |
| 356 | |
| 357 | static void mk_hbmap(t_hbdata *hb) |
| 358 | { |
| 359 | int i, j; |
| 360 | |
| 361 | snew(hb->hbmap, hb->d.nrd)(hb->hbmap) = save_calloc("hb->hbmap", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 361, (hb->d.nrd), sizeof(*(hb->hbmap))); |
| 362 | for (i = 0; (i < hb->d.nrd); i++) |
| 363 | { |
| 364 | snew(hb->hbmap[i], hb->a.nra)(hb->hbmap[i]) = save_calloc("hb->hbmap[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 364, (hb->a.nra), sizeof(*(hb->hbmap[i]))); |
| 365 | if (hb->hbmap[i] == NULL((void*)0)) |
| 366 | { |
| 367 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 367, "Could not allocate enough memory for hbmap"); |
| 368 | } |
| 369 | for (j = 0; (j > hb->a.nra); j++) |
| 370 | { |
| 371 | hb->hbmap[i][j] = NULL((void*)0); |
| 372 | } |
| 373 | } |
| 374 | } |
| 375 | |
| 376 | /* Consider redoing pHist so that is only stores transitions between |
| 377 | * periodicities and not the periodicity for all frames. This eats heaps of memory. */ |
| 378 | static void mk_per(t_hbdata *hb) |
| 379 | { |
| 380 | int i, j; |
| 381 | if (hb->bGem) |
| 382 | { |
| 383 | snew(hb->per->pHist, hb->d.nrd)(hb->per->pHist) = save_calloc("hb->per->pHist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 383, (hb->d.nrd), sizeof(*(hb->per->pHist))); |
| 384 | for (i = 0; i < hb->d.nrd; i++) |
| 385 | { |
| 386 | snew(hb->per->pHist[i], hb->a.nra)(hb->per->pHist[i]) = save_calloc("hb->per->pHist[i]" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 386, (hb->a.nra), sizeof(*(hb->per->pHist[i]))); |
| 387 | if (hb->per->pHist[i] == NULL((void*)0)) |
| 388 | { |
| 389 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 389, "Could not allocate enough memory for per->pHist"); |
| 390 | } |
| 391 | for (j = 0; j < hb->a.nra; j++) |
| 392 | { |
| 393 | clearPshift(&(hb->per->pHist[i][j])); |
| 394 | } |
| 395 | } |
| 396 | /* add the [0,0,0] shift to element 0 of p2i. */ |
| 397 | snew(hb->per->p2i, 1)(hb->per->p2i) = save_calloc("hb->per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 397, (1), sizeof(*(hb->per->p2i))); |
| 398 | clear_ivec(hb->per->p2i[0]); |
| 399 | hb->per->nper = 1; |
| 400 | } |
| 401 | } |
| 402 | |
| 403 | #ifdef HAVE_NN_LOOPS |
| 404 | static void mk_hbEmap (t_hbdata *hb, int n0) |
| 405 | { |
| 406 | int i, j, k; |
| 407 | hb->hbE.E = NULL((void*)0); |
| 408 | hb->hbE.nframes = 0; |
| 409 | snew(hb->hbE.E, hb->d.nrd)(hb->hbE.E) = save_calloc("hb->hbE.E", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 409, (hb->d.nrd), sizeof(*(hb->hbE.E))); |
| 410 | for (i = 0; i < hb->d.nrd; i++) |
| 411 | { |
| 412 | snew(hb->hbE.E[i], hb->a.nra)(hb->hbE.E[i]) = save_calloc("hb->hbE.E[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 412, (hb->a.nra), sizeof(*(hb->hbE.E[i]))); |
| 413 | for (j = 0; j < hb->a.nra; j++) |
| 414 | { |
| 415 | snew(hb->hbE.E[i][j], MAXHYDRO)(hb->hbE.E[i][j]) = save_calloc("hb->hbE.E[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 415, (4), sizeof(*(hb->hbE.E[i][j]))); |
| 416 | for (k = 0; k < MAXHYDRO4; k++) |
| 417 | { |
| 418 | hb->hbE.E[i][j][k] = NULL((void*)0); |
| 419 | } |
| 420 | } |
| 421 | } |
| 422 | hb->hbE.Etot = NULL((void*)0); |
| 423 | } |
| 424 | |
| 425 | static void free_hbEmap (t_hbdata *hb) |
| 426 | { |
| 427 | int i, j, k; |
| 428 | for (i = 0; i < hb->d.nrd; i++) |
| 429 | { |
| 430 | for (j = 0; j < hb->a.nra; j++) |
| 431 | { |
| 432 | for (k = 0; k < MAXHYDRO4; k++) |
| 433 | { |
| 434 | sfree(hb->hbE.E[i][j][k])save_free("hb->hbE.E[i][j][k]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 434, (hb->hbE.E[i][j][k])); |
| 435 | } |
| 436 | sfree(hb->hbE.E[i][j])save_free("hb->hbE.E[i][j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 436, (hb->hbE.E[i][j])); |
| 437 | } |
| 438 | sfree(hb->hbE.E[i])save_free("hb->hbE.E[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 438, (hb->hbE.E[i])); |
| 439 | } |
| 440 | sfree(hb->hbE.E)save_free("hb->hbE.E", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 440, (hb->hbE.E)); |
| 441 | sfree(hb->hbE.Etot)save_free("hb->hbE.Etot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 441, (hb->hbE.Etot)); |
| 442 | } |
| 443 | |
| 444 | static void addFramesNN(t_hbdata *hb, int frame) |
| 445 | { |
| 446 | |
| 447 | #define DELTAFRAMES_HBE 10 |
| 448 | |
| 449 | int d, a, h, nframes; |
| 450 | |
| 451 | if (frame >= hb->hbE.nframes) |
| 452 | { |
| 453 | nframes = hb->hbE.nframes + DELTAFRAMES_HBE; |
| 454 | srenew(hb->hbE.Etot, nframes)(hb->hbE.Etot) = save_realloc("hb->hbE.Etot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 454, (hb->hbE.Etot), (nframes), sizeof(*(hb->hbE.Etot ))); |
| 455 | |
| 456 | for (d = 0; d < hb->d.nrd; d++) |
| 457 | { |
| 458 | for (a = 0; a < hb->a.nra; a++) |
| 459 | { |
| 460 | for (h = 0; h < hb->d.nhydro[d]; h++) |
| 461 | { |
| 462 | srenew(hb->hbE.E[d][a][h], nframes)(hb->hbE.E[d][a][h]) = save_realloc("hb->hbE.E[d][a][h]" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 462, (hb->hbE.E[d][a][h]), (nframes), sizeof(*(hb->hbE .E[d][a][h]))); |
| 463 | } |
| 464 | } |
| 465 | } |
| 466 | |
| 467 | hb->hbE.nframes += DELTAFRAMES_HBE; |
| 468 | } |
| 469 | } |
| 470 | |
| 471 | static t_E calcHbEnergy(int d, int a, int h, rvec x[], t_EEst EEst, |
| 472 | matrix box, rvec hbox, t_donors *donors) |
| 473 | { |
| 474 | /* d - donor atom |
| 475 | * a - acceptor atom |
| 476 | * h - hydrogen |
| 477 | * alpha - angle between dipoles |
| 478 | * x[] - atomic positions |
| 479 | * EEst - the type of energy estimate (see enum in hbplugin.h) |
| 480 | * box - the box vectors \ |
| 481 | * hbox - half box lengths _These two are only needed for the pbc correction |
| 482 | */ |
| 483 | |
| 484 | t_E E; |
| 485 | rvec dist; |
| 486 | rvec dipole[2], xmol[3], xmean[2]; |
| 487 | int i; |
| 488 | real r, realE; |
| 489 | |
| 490 | if (d == a) |
| 491 | { |
| 492 | /* Self-interaction */ |
| 493 | return NONSENSE_E; |
| 494 | } |
| 495 | |
| 496 | switch (EEst) |
| 497 | { |
| 498 | case NN_BINARY: |
| 499 | /* This is a simple binary existence function that sets E=1 whenever |
| 500 | * the distance between the oxygens is equal too or less than 0.35 nm. |
| 501 | */ |
| 502 | rvec_sub(x[d], x[a], dist); |
| 503 | pbc_correct_gem(dist, box, hbox); |
| 504 | if (norm(dist) <= 0.35) |
| 505 | { |
| 506 | E = 1; |
| 507 | } |
| 508 | else |
| 509 | { |
| 510 | E = 0; |
| 511 | } |
| 512 | break; |
| 513 | |
| 514 | case NN_1_over_r3: |
| 515 | /* Negative potential energy of a dipole. |
| 516 | * E = -cos(alpha) * 1/r^3 */ |
| 517 | |
| 518 | copy_rvec(x[d], xmol[0]); /* donor */ |
| 519 | copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ |
| 520 | copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ |
| 521 | |
| 522 | svmul(15.9994*(1/1.008), xmol[0], xmean[0]); |
| 523 | rvec_inc(xmean[0], xmol[1]); |
| 524 | rvec_inc(xmean[0], xmol[2]); |
| 525 | for (i = 0; i < 3; i++) |
| 526 | { |
| 527 | xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
| 528 | } |
| 529 | |
| 530 | /* Assumes that all acceptors are also donors. */ |
| 531 | copy_rvec(x[a], xmol[0]); /* acceptor */ |
| 532 | copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ |
| 533 | copy_rvec(x[donors->hydro[donors->dptr[a]][1]], xmol[2]); /* hydrogen */ |
| 534 | |
| 535 | |
| 536 | svmul(15.9994*(1/1.008), xmol[0], xmean[1]); |
| 537 | rvec_inc(xmean[1], xmol[1]); |
| 538 | rvec_inc(xmean[1], xmol[2]); |
| 539 | for (i = 0; i < 3; i++) |
| 540 | { |
| 541 | xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
| 542 | } |
| 543 | |
| 544 | rvec_sub(xmean[0], xmean[1], dist); |
| 545 | pbc_correct_gem(dist, box, hbox); |
| 546 | r = norm(dist); |
| 547 | |
| 548 | realE = pow(r, -3.0); |
| 549 | E = (t_E)(SCALEFACTOR_E * realE); |
| 550 | break; |
| 551 | |
| 552 | case NN_dipole: |
| 553 | /* Negative potential energy of a (unpolarizable) dipole. |
| 554 | * E = -cos(alpha) * 1/r^3 */ |
| 555 | clear_rvec(dipole[1]); |
| 556 | clear_rvec(dipole[0]); |
| 557 | |
| 558 | copy_rvec(x[d], xmol[0]); /* donor */ |
| 559 | copy_rvec(x[donors->hydro[donors->dptr[d]][0]], xmol[1]); /* hydrogen */ |
| 560 | copy_rvec(x[donors->hydro[donors->dptr[d]][1]], xmol[2]); /* hydrogen */ |
| 561 | |
| 562 | rvec_inc(dipole[0], xmol[1]); |
| 563 | rvec_inc(dipole[0], xmol[2]); |
| 564 | for (i = 0; i < 3; i++) |
| 565 | { |
| 566 | dipole[0][i] *= 0.5; |
| 567 | } |
| 568 | rvec_dec(dipole[0], xmol[0]); |
| 569 | |
| 570 | svmul(15.9994*(1/1.008), xmol[0], xmean[0]); |
| 571 | rvec_inc(xmean[0], xmol[1]); |
| 572 | rvec_inc(xmean[0], xmol[2]); |
| 573 | for (i = 0; i < 3; i++) |
| 574 | { |
| 575 | xmean[0][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
| 576 | } |
| 577 | |
| 578 | /* Assumes that all acceptors are also donors. */ |
| 579 | copy_rvec(x[a], xmol[0]); /* acceptor */ |
| 580 | copy_rvec(x[donors->hydro[donors->dptr[a]][0]], xmol[1]); /* hydrogen */ |
| 581 | copy_rvec(x[donors->hydro[donors->dptr[a]][2]], xmol[2]); /* hydrogen */ |
| 582 | |
| 583 | |
| 584 | rvec_inc(dipole[1], xmol[1]); |
| 585 | rvec_inc(dipole[1], xmol[2]); |
| 586 | for (i = 0; i < 3; i++) |
| 587 | { |
| 588 | dipole[1][i] *= 0.5; |
| 589 | } |
| 590 | rvec_dec(dipole[1], xmol[0]); |
| 591 | |
| 592 | svmul(15.9994*(1/1.008), xmol[0], xmean[1]); |
| 593 | rvec_inc(xmean[1], xmol[1]); |
| 594 | rvec_inc(xmean[1], xmol[2]); |
| 595 | for (i = 0; i < 3; i++) |
| 596 | { |
| 597 | xmean[1][i] /= (15.9994 + 1.008 + 1.008)/1.008; |
| 598 | } |
| 599 | |
| 600 | rvec_sub(xmean[0], xmean[1], dist); |
| 601 | pbc_correct_gem(dist, box, hbox); |
| 602 | r = norm(dist); |
| 603 | |
| 604 | double cosalpha = cos_angle(dipole[0], dipole[1]); |
| 605 | realE = cosalpha * pow(r, -3.0); |
| 606 | E = (t_E)(SCALEFACTOR_E * realE); |
| 607 | break; |
| 608 | |
| 609 | default: |
| 610 | printf("Can't do that type of energy estimate: %i\n.", EEst); |
| 611 | E = NONSENSE_E; |
| 612 | } |
| 613 | |
| 614 | return E; |
| 615 | } |
| 616 | |
| 617 | static void storeHbEnergy(t_hbdata *hb, int d, int a, int h, t_E E, int frame) |
| 618 | { |
| 619 | /* hb - hbond data structure |
| 620 | d - donor |
| 621 | a - acceptor |
| 622 | h - hydrogen |
| 623 | E - estimate of the energy |
| 624 | frame - the current frame. |
| 625 | */ |
| 626 | |
| 627 | /* Store the estimated energy */ |
| 628 | if (E == NONSENSE_E) |
| 629 | { |
| 630 | E = 0; |
| 631 | } |
| 632 | |
| 633 | hb->hbE.E[d][a][h][frame] = E; |
| 634 | |
| 635 | #pragma omp critical |
| 636 | { |
| 637 | hb->hbE.Etot[frame] += E; |
| 638 | } |
| 639 | } |
| 640 | #endif /* HAVE_NN_LOOPS */ |
| 641 | |
| 642 | |
| 643 | /* Finds -v[] in the periodicity index */ |
| 644 | static int findMirror(PSTYPEint p, ivec v[], PSTYPEint nper) |
| 645 | { |
| 646 | PSTYPEint i; |
| 647 | ivec u; |
| 648 | for (i = 0; i < nper; i++) |
| 649 | { |
| 650 | if (v[i][XX0] == -(v[p][XX0]) && |
| 651 | v[i][YY1] == -(v[p][YY1]) && |
| 652 | v[i][ZZ2] == -(v[p][ZZ2])) |
| 653 | { |
| 654 | return (int)i; |
| 655 | } |
| 656 | } |
| 657 | printf("Couldn't find mirror of [%i, %i, %i], index \n", |
| 658 | v[p][XX0], |
| 659 | v[p][YY1], |
| 660 | v[p][ZZ2]); |
| 661 | return -1; |
| 662 | } |
| 663 | |
| 664 | |
| 665 | static void add_frames(t_hbdata *hb, int nframes) |
| 666 | { |
| 667 | int i, j, k, l; |
| 668 | |
| 669 | if (nframes >= hb->max_frames) |
| 670 | { |
| 671 | hb->max_frames += 4096; |
| 672 | srenew(hb->time, hb->max_frames)(hb->time) = save_realloc("hb->time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 672, (hb->time), (hb->max_frames), sizeof(*(hb->time ))); |
| 673 | srenew(hb->nhb, hb->max_frames)(hb->nhb) = save_realloc("hb->nhb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 673, (hb->nhb), (hb->max_frames), sizeof(*(hb->nhb ))); |
| 674 | srenew(hb->ndist, hb->max_frames)(hb->ndist) = save_realloc("hb->ndist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 674, (hb->ndist), (hb->max_frames), sizeof(*(hb->ndist ))); |
| 675 | srenew(hb->n_bound, hb->max_frames)(hb->n_bound) = save_realloc("hb->n_bound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 675, (hb->n_bound), (hb->max_frames), sizeof(*(hb-> n_bound))); |
| 676 | srenew(hb->nhx, hb->max_frames)(hb->nhx) = save_realloc("hb->nhx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 676, (hb->nhx), (hb->max_frames), sizeof(*(hb->nhx ))); |
| 677 | if (hb->bDAnr) |
| 678 | { |
| 679 | srenew(hb->danr, hb->max_frames)(hb->danr) = save_realloc("hb->danr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 679, (hb->danr), (hb->max_frames), sizeof(*(hb->danr ))); |
| 680 | } |
| 681 | } |
| 682 | hb->nframes = nframes; |
| 683 | } |
| 684 | |
| 685 | #define OFFSET(frame)(frame / 32) (frame / 32) |
| 686 | #define MASK(frame)(1 << (frame % 32)) (1 << (frame % 32)) |
| 687 | |
| 688 | static void _set_hb(unsigned int hbexist[], unsigned int frame, gmx_bool bValue) |
| 689 | { |
| 690 | if (bValue) |
| 691 | { |
| 692 | hbexist[OFFSET(frame)(frame / 32)] |= MASK(frame)(1 << (frame % 32)); |
| 693 | } |
| 694 | else |
| 695 | { |
| 696 | hbexist[OFFSET(frame)(frame / 32)] &= ~MASK(frame)(1 << (frame % 32)); |
| 697 | } |
| 698 | } |
| 699 | |
| 700 | static gmx_bool is_hb(unsigned int hbexist[], int frame) |
| 701 | { |
| 702 | return ((hbexist[OFFSET(frame)(frame / 32)] & MASK(frame)(1 << (frame % 32))) != 0) ? 1 : 0; |
| 703 | } |
| 704 | |
| 705 | static void set_hb(t_hbdata *hb, int id, int ih, int ia, int frame, int ihb) |
| 706 | { |
| 707 | unsigned int *ghptr = NULL((void*)0); |
| 708 | |
| 709 | if (ihb == hbHB) |
| 710 | { |
| 711 | ghptr = hb->hbmap[id][ia]->h[ih]; |
| 712 | } |
| 713 | else if (ihb == hbDist) |
| 714 | { |
| 715 | ghptr = hb->hbmap[id][ia]->g[ih]; |
| 716 | } |
| 717 | else |
| 718 | { |
| 719 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 719, "Incomprehensible iValue %d in set_hb", ihb); |
| 720 | } |
| 721 | |
| 722 | _set_hb(ghptr, frame-hb->hbmap[id][ia]->n0, TRUE1); |
| 723 | } |
| 724 | |
| 725 | static void addPshift(t_pShift *pHist, PSTYPEint p, int frame) |
| 726 | { |
| 727 | if (pHist->len == 0) |
| 728 | { |
| 729 | snew(pHist->frame, 1)(pHist->frame) = save_calloc("pHist->frame", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 729, (1), sizeof(*(pHist->frame))); |
| 730 | snew(pHist->p, 1)(pHist->p) = save_calloc("pHist->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 730, (1), sizeof(*(pHist->p))); |
| 731 | pHist->len = 1; |
| 732 | pHist->frame[0] = frame; |
| 733 | pHist->p[0] = p; |
| 734 | return; |
| 735 | } |
| 736 | else |
| 737 | if (pHist->p[pHist->len-1] != p) |
| 738 | { |
| 739 | pHist->len++; |
| 740 | srenew(pHist->frame, pHist->len)(pHist->frame) = save_realloc("pHist->frame", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 740, (pHist->frame), (pHist->len), sizeof(*(pHist-> frame))); |
| 741 | srenew(pHist->p, pHist->len)(pHist->p) = save_realloc("pHist->p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 741, (pHist->p), (pHist->len), sizeof(*(pHist->p)) ); |
| 742 | pHist->frame[pHist->len-1] = frame; |
| 743 | pHist->p[pHist->len-1] = p; |
| 744 | } /* Otherwise, there is no transition. */ |
| 745 | return; |
| 746 | } |
| 747 | |
| 748 | static PSTYPEint getPshift(t_pShift pHist, int frame) |
| 749 | { |
| 750 | int f, i; |
| 751 | |
| 752 | if (pHist.len == 0 |
| 753 | || (pHist.len > 0 && pHist.frame[0] > frame)) |
| 754 | { |
| 755 | return -1; |
| 756 | } |
| 757 | |
| 758 | for (i = 0; i < pHist.len; i++) |
| 759 | { |
| 760 | f = pHist.frame[i]; |
| 761 | if (f == frame) |
| 762 | { |
| 763 | return pHist.p[i]; |
| 764 | } |
| 765 | if (f > frame) |
| 766 | { |
| 767 | return pHist.p[i-1]; |
| 768 | } |
| 769 | } |
| 770 | |
| 771 | /* It seems that frame is after the last periodic transition. Return the last periodicity. */ |
| 772 | return pHist.p[pHist.len-1]; |
| 773 | } |
| 774 | |
| 775 | static void add_ff(t_hbdata *hbd, int id, int h, int ia, int frame, int ihb, PSTYPEint p) |
| 776 | { |
| 777 | int i, j, n; |
| 778 | t_hbond *hb = hbd->hbmap[id][ia]; |
| 779 | int maxhydro = min(hbd->maxhydro, hbd->d.nhydro[id])(((hbd->maxhydro) < (hbd->d.nhydro[id])) ? (hbd-> maxhydro) : (hbd->d.nhydro[id]) ); |
| 780 | int wlen = hbd->wordlen; |
| 781 | int delta = 32*wlen; |
| 782 | gmx_bool bGem = hbd->bGem; |
| 783 | |
| 784 | if (!hb->h[0]) |
| 785 | { |
| 786 | hb->n0 = frame; |
| 787 | hb->maxframes = delta; |
| 788 | for (i = 0; (i < maxhydro); i++) |
| 789 | { |
| 790 | snew(hb->h[i], hb->maxframes/wlen)(hb->h[i]) = save_calloc("hb->h[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 790, (hb->maxframes/wlen), sizeof(*(hb->h[i]))); |
| 791 | snew(hb->g[i], hb->maxframes/wlen)(hb->g[i]) = save_calloc("hb->g[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 791, (hb->maxframes/wlen), sizeof(*(hb->g[i]))); |
| 792 | } |
| 793 | } |
| 794 | else |
| 795 | { |
| 796 | hb->nframes = frame-hb->n0; |
| 797 | /* We need a while loop here because hbonds may be returning |
| 798 | * after a long time. |
| 799 | */ |
| 800 | while (hb->nframes >= hb->maxframes) |
| 801 | { |
| 802 | n = hb->maxframes + delta; |
| 803 | for (i = 0; (i < maxhydro); i++) |
| 804 | { |
| 805 | srenew(hb->h[i], n/wlen)(hb->h[i]) = save_realloc("hb->h[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 805, (hb->h[i]), (n/wlen), sizeof(*(hb->h[i]))); |
| 806 | srenew(hb->g[i], n/wlen)(hb->g[i]) = save_realloc("hb->g[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 806, (hb->g[i]), (n/wlen), sizeof(*(hb->g[i]))); |
| 807 | for (j = hb->maxframes/wlen; (j < n/wlen); j++) |
| 808 | { |
| 809 | hb->h[i][j] = 0; |
| 810 | hb->g[i][j] = 0; |
| 811 | } |
| 812 | } |
| 813 | |
| 814 | hb->maxframes = n; |
| 815 | } |
| 816 | } |
| 817 | if (frame >= 0) |
| 818 | { |
| 819 | set_hb(hbd, id, h, ia, frame, ihb); |
| 820 | if (bGem) |
| 821 | { |
| 822 | if (p >= hbd->per->nper) |
| 823 | { |
| 824 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 824, "invalid shift: p=%u, nper=%u", p, hbd->per->nper); |
| 825 | } |
| 826 | else |
| 827 | { |
| 828 | addPshift(&(hbd->per->pHist[id][ia]), p, frame); |
| 829 | } |
| 830 | |
| 831 | } |
| 832 | } |
| 833 | |
| 834 | } |
| 835 | |
| 836 | static void inc_nhbonds(t_donors *ddd, int d, int h) |
| 837 | { |
| 838 | int j; |
| 839 | int dptr = ddd->dptr[d]; |
| 840 | |
| 841 | for (j = 0; (j < ddd->nhydro[dptr]); j++) |
| 842 | { |
| 843 | if (ddd->hydro[dptr][j] == h) |
| 844 | { |
| 845 | ddd->nhbonds[dptr][j]++; |
| 846 | break; |
| 847 | } |
| 848 | } |
| 849 | if (j == ddd->nhydro[dptr]) |
| 850 | { |
| 851 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 851, "No such hydrogen %d on donor %d\n", h+1, d+1); |
| 852 | } |
| 853 | } |
| 854 | |
| 855 | static int _acceptor_index(t_acceptors *a, int grp, atom_id i, |
| 856 | const char *file, int line) |
| 857 | { |
| 858 | int ai = a->aptr[i]; |
| 859 | |
| 860 | if (a->grp[ai] != grp) |
| 861 | { |
| 862 | if (debug && bDebug) |
| 863 | { |
| 864 | fprintf(debug, "Acc. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n", |
| 865 | ai, a->grp[ai], grp, file, line); |
| 866 | } |
| 867 | return NOTSET-12345; |
| 868 | } |
| 869 | else |
| 870 | { |
| 871 | return ai; |
| 872 | } |
| 873 | } |
| 874 | #define acceptor_index(a, grp, i)_acceptor_index(a, grp, i, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 874) _acceptor_index(a, grp, i, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c", __LINE__874) |
| 875 | |
| 876 | static int _donor_index(t_donors *d, int grp, atom_id i, const char *file, int line) |
| 877 | { |
| 878 | int di = d->dptr[i]; |
| 879 | |
| 880 | if (di == NOTSET-12345) |
| 881 | { |
| 882 | return NOTSET-12345; |
| 883 | } |
| 884 | |
| 885 | if (d->grp[di] != grp) |
| 886 | { |
| 887 | if (debug && bDebug) |
| 888 | { |
| 889 | fprintf(debug, "Don. group inconsist.. grp[%d] = %d, grp = %d (%s, %d)\n", |
| 890 | di, d->grp[di], grp, file, line); |
| 891 | } |
| 892 | return NOTSET-12345; |
| 893 | } |
| 894 | else |
| 895 | { |
| 896 | return di; |
| 897 | } |
| 898 | } |
| 899 | #define donor_index(d, grp, i)_donor_index(d, grp, i, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 899) _donor_index(d, grp, i, __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c", __LINE__899) |
| 900 | |
| 901 | static gmx_bool isInterchangable(t_hbdata *hb, int d, int a, int grpa, int grpd) |
| 902 | { |
| 903 | /* g_hbond doesn't allow overlapping groups */ |
| 904 | if (grpa != grpd) |
| 905 | { |
| 906 | return FALSE0; |
| 907 | } |
| 908 | return |
| 909 | donor_index(&hb->d, grpd, a)_donor_index(&hb->d, grpd, a, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 909) != NOTSET-12345 |
| 910 | && acceptor_index(&hb->a, grpa, d)_acceptor_index(&hb->a, grpa, d, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 910) != NOTSET-12345; |
| 911 | } |
| 912 | |
| 913 | |
| 914 | static void add_hbond(t_hbdata *hb, int d, int a, int h, int grpd, int grpa, |
| 915 | int frame, gmx_bool bMerge, int ihb, gmx_bool bContact, PSTYPEint p) |
| 916 | { |
| 917 | int k, id, ia, hh; |
| 918 | gmx_bool daSwap = FALSE0; |
| 919 | |
| 920 | if ((id = hb->d.dptr[d]) == NOTSET-12345) |
| 921 | { |
| 922 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 922, "No donor atom %d", d+1); |
| 923 | } |
| 924 | else if (grpd != hb->d.grp[id]) |
| 925 | { |
| 926 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 926, "Inconsistent donor groups, %d iso %d, atom %d", |
| 927 | grpd, hb->d.grp[id], d+1); |
| 928 | } |
| 929 | if ((ia = hb->a.aptr[a]) == NOTSET-12345) |
| 930 | { |
| 931 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 931, "No acceptor atom %d", a+1); |
| 932 | } |
| 933 | else if (grpa != hb->a.grp[ia]) |
| 934 | { |
| 935 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 935, "Inconsistent acceptor groups, %d iso %d, atom %d", |
| 936 | grpa, hb->a.grp[ia], a+1); |
| 937 | } |
| 938 | |
| 939 | if (bMerge) |
| 940 | { |
| 941 | |
| 942 | if (isInterchangable(hb, d, a, grpd, grpa) && d > a) |
| 943 | /* Then swap identity so that the id of d is lower then that of a. |
| 944 | * |
| 945 | * This should really be redundant by now, as is_hbond() now ought to return |
| 946 | * hbNo in the cases where this conditional is TRUE. */ |
| 947 | { |
| 948 | daSwap = TRUE1; |
| 949 | k = d; |
| 950 | d = a; |
| 951 | a = k; |
| 952 | |
| 953 | /* Now repeat donor/acc check. */ |
| 954 | if ((id = hb->d.dptr[d]) == NOTSET-12345) |
| 955 | { |
| 956 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 956, "No donor atom %d", d+1); |
| 957 | } |
| 958 | else if (grpd != hb->d.grp[id]) |
| 959 | { |
| 960 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 960, "Inconsistent donor groups, %d iso %d, atom %d", |
| 961 | grpd, hb->d.grp[id], d+1); |
| 962 | } |
| 963 | if ((ia = hb->a.aptr[a]) == NOTSET-12345) |
| 964 | { |
| 965 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 965, "No acceptor atom %d", a+1); |
| 966 | } |
| 967 | else if (grpa != hb->a.grp[ia]) |
| 968 | { |
| 969 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 969, "Inconsistent acceptor groups, %d iso %d, atom %d", |
| 970 | grpa, hb->a.grp[ia], a+1); |
| 971 | } |
| 972 | } |
| 973 | } |
| 974 | |
| 975 | if (hb->hbmap) |
| 976 | { |
| 977 | /* Loop over hydrogens to find which hydrogen is in this particular HB */ |
| 978 | if ((ihb == hbHB) && !bMerge && !bContact) |
| 979 | { |
| 980 | for (k = 0; (k < hb->d.nhydro[id]); k++) |
| 981 | { |
| 982 | if (hb->d.hydro[id][k] == h) |
| 983 | { |
| 984 | break; |
| 985 | } |
| 986 | } |
| 987 | if (k == hb->d.nhydro[id]) |
| 988 | { |
| 989 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 989, "Donor %d does not have hydrogen %d (a = %d)", |
| 990 | d+1, h+1, a+1); |
| 991 | } |
| 992 | } |
| 993 | else |
| 994 | { |
| 995 | k = 0; |
| 996 | } |
| 997 | |
| 998 | if (hb->bHBmap) |
| 999 | { |
| 1000 | |
| 1001 | #pragma omp critical |
| 1002 | { |
| 1003 | if (hb->hbmap[id][ia] == NULL((void*)0)) |
| 1004 | { |
| 1005 | snew(hb->hbmap[id][ia], 1)(hb->hbmap[id][ia]) = save_calloc("hb->hbmap[id][ia]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1005, (1), sizeof(*(hb->hbmap[id][ia]))); |
| 1006 | snew(hb->hbmap[id][ia]->h, hb->maxhydro)(hb->hbmap[id][ia]->h) = save_calloc("hb->hbmap[id][ia]->h" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1006, (hb->maxhydro), sizeof(*(hb->hbmap[id][ia]-> h))); |
| 1007 | snew(hb->hbmap[id][ia]->g, hb->maxhydro)(hb->hbmap[id][ia]->g) = save_calloc("hb->hbmap[id][ia]->g" , "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1007, (hb->maxhydro), sizeof(*(hb->hbmap[id][ia]-> g))); |
| 1008 | } |
| 1009 | add_ff(hb, id, k, ia, frame, ihb, p); |
| 1010 | } |
| 1011 | } |
| 1012 | |
| 1013 | /* Strange construction with frame >=0 is a relic from old code |
| 1014 | * for selected hbond analysis. It may be necessary again if that |
| 1015 | * is made to work again. |
| 1016 | */ |
| 1017 | if (frame >= 0) |
| 1018 | { |
| 1019 | hh = hb->hbmap[id][ia]->history[k]; |
| 1020 | if (ihb == hbHB) |
| 1021 | { |
| 1022 | hb->nhb[frame]++; |
| 1023 | if (!(ISHB(hh)(((hh) & 2) == 2))) |
| 1024 | { |
| 1025 | hb->hbmap[id][ia]->history[k] = hh | 2; |
| 1026 | hb->nrhb++; |
| 1027 | } |
| 1028 | } |
| 1029 | else |
| 1030 | { |
| 1031 | if (ihb == hbDist) |
| 1032 | { |
| 1033 | hb->ndist[frame]++; |
| 1034 | if (!(ISDIST(hh)(((hh) & 1) == 1))) |
| 1035 | { |
| 1036 | hb->hbmap[id][ia]->history[k] = hh | 1; |
| 1037 | hb->nrdist++; |
| 1038 | } |
| 1039 | } |
| 1040 | } |
| 1041 | } |
| 1042 | } |
| 1043 | else |
| 1044 | { |
| 1045 | if (frame >= 0) |
| 1046 | { |
| 1047 | if (ihb == hbHB) |
| 1048 | { |
| 1049 | hb->nhb[frame]++; |
| 1050 | } |
| 1051 | else |
| 1052 | { |
| 1053 | if (ihb == hbDist) |
| 1054 | { |
| 1055 | hb->ndist[frame]++; |
| 1056 | } |
| 1057 | } |
| 1058 | } |
| 1059 | } |
| 1060 | if (bMerge && daSwap) |
| 1061 | { |
| 1062 | h = hb->d.hydro[id][0]; |
| 1063 | } |
| 1064 | /* Increment number if HBonds per H */ |
| 1065 | if (ihb == hbHB && !bContact) |
| 1066 | { |
| 1067 | inc_nhbonds(&(hb->d), d, h); |
| 1068 | } |
| 1069 | } |
| 1070 | |
| 1071 | static char *mkatomname(t_atoms *atoms, int i) |
| 1072 | { |
| 1073 | static char buf[32]; |
| 1074 | int rnr; |
| 1075 | |
| 1076 | rnr = atoms->atom[i].resind; |
| 1077 | sprintf(buf, "%4s%d%-4s", |
| 1078 | *atoms->resinfo[rnr].name, atoms->resinfo[rnr].nr, *atoms->atomname[i]); |
| 1079 | |
| 1080 | return buf; |
| 1081 | } |
| 1082 | |
| 1083 | static void gen_datable(atom_id *index, int isize, unsigned char *datable, int natoms) |
| 1084 | { |
| 1085 | /* Generates table of all atoms and sets the ingroup bit for atoms in index[] */ |
| 1086 | int i; |
| 1087 | |
| 1088 | for (i = 0; i < isize; i++) |
| 1089 | { |
| 1090 | if (index[i] >= natoms) |
| 1091 | { |
| 1092 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1092, "Atom has index %d larger than number of atoms %d.", index[i], natoms); |
| 1093 | } |
| 1094 | datable[index[i]] |= INGROUP; |
| 1095 | } |
| 1096 | } |
| 1097 | |
| 1098 | static void clear_datable_grp(unsigned char *datable, int size) |
| 1099 | { |
| 1100 | /* Clears group information from the table */ |
| 1101 | int i; |
| 1102 | const char mask = !(char)INGROUP; |
| 1103 | if (size > 0) |
| 1104 | { |
| 1105 | for (i = 0; i < size; i++) |
| 1106 | { |
| 1107 | datable[i] &= mask; |
| 1108 | } |
| 1109 | } |
| 1110 | } |
| 1111 | |
| 1112 | static void add_acc(t_acceptors *a, int ia, int grp) |
| 1113 | { |
| 1114 | if (a->nra >= a->max_nra) |
| 1115 | { |
| 1116 | a->max_nra += 16; |
| 1117 | srenew(a->acc, a->max_nra)(a->acc) = save_realloc("a->acc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1117, (a->acc), (a->max_nra), sizeof(*(a->acc))); |
| 1118 | srenew(a->grp, a->max_nra)(a->grp) = save_realloc("a->grp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1118, (a->grp), (a->max_nra), sizeof(*(a->grp))); |
| 1119 | } |
| 1120 | a->grp[a->nra] = grp; |
| 1121 | a->acc[a->nra++] = ia; |
| 1122 | } |
| 1123 | |
| 1124 | static void search_acceptors(t_topology *top, int isize, |
| 1125 | atom_id *index, t_acceptors *a, int grp, |
| 1126 | gmx_bool bNitAcc, |
| 1127 | gmx_bool bContact, gmx_bool bDoIt, unsigned char *datable) |
| 1128 | { |
| 1129 | int i, n; |
| 1130 | |
| 1131 | if (bDoIt) |
| 1132 | { |
| 1133 | for (i = 0; (i < isize); i++) |
| 1134 | { |
| 1135 | n = index[i]; |
| 1136 | if ((bContact || |
| 1137 | (((*top->atoms.atomname[n])[0] == 'O') || |
| 1138 | (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) && |
| 1139 | ISINGRP(datable[n])(((datable[n]) & 4) == 4)) |
| 1140 | { |
| 1141 | datable[n] |= ACC; /* set the atom's acceptor flag in datable. */ |
| 1142 | add_acc(a, n, grp); |
| 1143 | } |
| 1144 | } |
| 1145 | } |
| 1146 | snew(a->aptr, top->atoms.nr)(a->aptr) = save_calloc("a->aptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1146, (top->atoms.nr), sizeof(*(a->aptr))); |
| 1147 | for (i = 0; (i < top->atoms.nr); i++) |
| 1148 | { |
| 1149 | a->aptr[i] = NOTSET-12345; |
| 1150 | } |
| 1151 | for (i = 0; (i < a->nra); i++) |
| 1152 | { |
| 1153 | a->aptr[a->acc[i]] = i; |
| 1154 | } |
| 1155 | } |
| 1156 | |
| 1157 | static void add_h2d(int id, int ih, t_donors *ddd) |
| 1158 | { |
| 1159 | int i; |
| 1160 | |
| 1161 | for (i = 0; (i < ddd->nhydro[id]); i++) |
| 1162 | { |
| 1163 | if (ddd->hydro[id][i] == ih) |
| 1164 | { |
| 1165 | printf("Hm. This isn't the first time I found this donor (%d,%d)\n", |
| 1166 | ddd->don[id], ih); |
| 1167 | break; |
| 1168 | } |
| 1169 | } |
| 1170 | if (i == ddd->nhydro[id]) |
| 1171 | { |
| 1172 | if (ddd->nhydro[id] >= MAXHYDRO4) |
| 1173 | { |
| 1174 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1174, "Donor %d has more than %d hydrogens!", |
| 1175 | ddd->don[id], MAXHYDRO4); |
| 1176 | } |
| 1177 | ddd->hydro[id][i] = ih; |
| 1178 | ddd->nhydro[id]++; |
| 1179 | } |
| 1180 | } |
| 1181 | |
| 1182 | static void add_dh(t_donors *ddd, int id, int ih, int grp, unsigned char *datable) |
| 1183 | { |
| 1184 | int i; |
| 1185 | |
| 1186 | if (ISDON(datable[id])(((datable[id]) & 2) == 2) || !datable) |
| 1187 | { |
| 1188 | if (ddd->dptr[id] == NOTSET-12345) /* New donor */ |
| 1189 | { |
| 1190 | i = ddd->nrd; |
| 1191 | ddd->dptr[id] = i; |
| 1192 | } |
| 1193 | else |
| 1194 | { |
| 1195 | i = ddd->dptr[id]; |
| 1196 | } |
| 1197 | |
| 1198 | if (i == ddd->nrd) |
| 1199 | { |
| 1200 | if (ddd->nrd >= ddd->max_nrd) |
| 1201 | { |
| 1202 | ddd->max_nrd += 128; |
| 1203 | srenew(ddd->don, ddd->max_nrd)(ddd->don) = save_realloc("ddd->don", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1203, (ddd->don), (ddd->max_nrd), sizeof(*(ddd->don ))); |
| 1204 | srenew(ddd->nhydro, ddd->max_nrd)(ddd->nhydro) = save_realloc("ddd->nhydro", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1204, (ddd->nhydro), (ddd->max_nrd), sizeof(*(ddd-> nhydro))); |
| 1205 | srenew(ddd->hydro, ddd->max_nrd)(ddd->hydro) = save_realloc("ddd->hydro", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1205, (ddd->hydro), (ddd->max_nrd), sizeof(*(ddd-> hydro))); |
| 1206 | srenew(ddd->nhbonds, ddd->max_nrd)(ddd->nhbonds) = save_realloc("ddd->nhbonds", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1206, (ddd->nhbonds), (ddd->max_nrd), sizeof(*(ddd-> nhbonds))); |
| 1207 | srenew(ddd->grp, ddd->max_nrd)(ddd->grp) = save_realloc("ddd->grp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1207, (ddd->grp), (ddd->max_nrd), sizeof(*(ddd->grp ))); |
| 1208 | } |
| 1209 | ddd->don[ddd->nrd] = id; |
| 1210 | ddd->nhydro[ddd->nrd] = 0; |
| 1211 | ddd->grp[ddd->nrd] = grp; |
| 1212 | ddd->nrd++; |
| 1213 | } |
| 1214 | else |
| 1215 | { |
| 1216 | ddd->don[i] = id; |
| 1217 | } |
| 1218 | add_h2d(i, ih, ddd); |
| 1219 | } |
| 1220 | else |
| 1221 | if (datable) |
| 1222 | { |
| 1223 | printf("Warning: Atom %d is not in the d/a-table!\n", id); |
| 1224 | } |
| 1225 | } |
| 1226 | |
| 1227 | static void search_donors(t_topology *top, int isize, atom_id *index, |
| 1228 | t_donors *ddd, int grp, gmx_bool bContact, gmx_bool bDoIt, |
| 1229 | unsigned char *datable) |
| 1230 | { |
| 1231 | int i, j, nra, n; |
| 1232 | t_functype func_type; |
| 1233 | t_ilist *interaction; |
| 1234 | atom_id nr1, nr2, nr3; |
| 1235 | gmx_bool stop; |
| 1236 | |
| 1237 | if (!ddd->dptr) |
| 1238 | { |
| 1239 | snew(ddd->dptr, top->atoms.nr)(ddd->dptr) = save_calloc("ddd->dptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1239, (top->atoms.nr), sizeof(*(ddd->dptr))); |
| 1240 | for (i = 0; (i < top->atoms.nr); i++) |
| 1241 | { |
| 1242 | ddd->dptr[i] = NOTSET-12345; |
| 1243 | } |
| 1244 | } |
| 1245 | |
| 1246 | if (bContact) |
| 1247 | { |
| 1248 | if (bDoIt) |
| 1249 | { |
| 1250 | for (i = 0; (i < isize); i++) |
| 1251 | { |
| 1252 | datable[index[i]] |= DON; |
| 1253 | add_dh(ddd, index[i], -1, grp, datable); |
| 1254 | } |
| 1255 | } |
| 1256 | } |
| 1257 | else |
| 1258 | { |
| 1259 | for (func_type = 0; (func_type < F_NRE); func_type++) |
| 1260 | { |
| 1261 | interaction = &(top->idef.il[func_type]); |
| 1262 | if (func_type == F_POSRES || func_type == F_FBPOSRES) |
| 1263 | { |
| 1264 | /* The ilist looks strange for posre. Bug in grompp? |
| 1265 | * We don't need posre interactions for hbonds anyway.*/ |
| 1266 | continue; |
| 1267 | } |
| 1268 | for (i = 0; i < interaction->nr; |
| 1269 | i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) |
| 1270 | { |
| 1271 | /* next function */ |
| 1272 | if (func_type != top->idef.functype[interaction->iatoms[i]]) |
| 1273 | { |
| 1274 | fprintf(stderrstderr, "Error in func_type %s", |
| 1275 | interaction_function[func_type].longname); |
| 1276 | continue; |
| 1277 | } |
| 1278 | |
| 1279 | /* check out this functype */ |
| 1280 | if (func_type == F_SETTLE) |
| 1281 | { |
| 1282 | nr1 = interaction->iatoms[i+1]; |
| 1283 | nr2 = interaction->iatoms[i+2]; |
| 1284 | nr3 = interaction->iatoms[i+3]; |
| 1285 | |
| 1286 | if (ISINGRP(datable[nr1])(((datable[nr1]) & 4) == 4)) |
| 1287 | { |
| 1288 | if (ISINGRP(datable[nr2])(((datable[nr2]) & 4) == 4)) |
| 1289 | { |
| 1290 | datable[nr1] |= DON; |
| 1291 | add_dh(ddd, nr1, nr1+1, grp, datable); |
| 1292 | } |
| 1293 | if (ISINGRP(datable[nr3])(((datable[nr3]) & 4) == 4)) |
| 1294 | { |
| 1295 | datable[nr1] |= DON; |
| 1296 | add_dh(ddd, nr1, nr1+2, grp, datable); |
| 1297 | } |
| 1298 | } |
| 1299 | } |
| 1300 | else if (IS_CHEMBOND(func_type)(interaction_function[(func_type)].nratoms == 2 && (interaction_function [(func_type)].flags & 1<<3))) |
| 1301 | { |
| 1302 | for (j = 0; j < 2; j++) |
| 1303 | { |
| 1304 | nr1 = interaction->iatoms[i+1+j]; |
| 1305 | nr2 = interaction->iatoms[i+2-j]; |
| 1306 | if ((*top->atoms.atomname[nr1][0] == 'H') && |
| 1307 | ((*top->atoms.atomname[nr2][0] == 'O') || |
| 1308 | (*top->atoms.atomname[nr2][0] == 'N')) && |
| 1309 | ISINGRP(datable[nr1])(((datable[nr1]) & 4) == 4) && ISINGRP(datable[nr2])(((datable[nr2]) & 4) == 4)) |
| 1310 | { |
| 1311 | datable[nr2] |= DON; |
| 1312 | add_dh(ddd, nr2, nr1, grp, datable); |
| 1313 | } |
| 1314 | } |
| 1315 | } |
| 1316 | } |
| 1317 | } |
| 1318 | #ifdef SAFEVSITES |
| 1319 | for (func_type = 0; func_type < F_NRE; func_type++) |
| 1320 | { |
| 1321 | interaction = &top->idef.il[func_type]; |
| 1322 | for (i = 0; i < interaction->nr; |
| 1323 | i += interaction_function[top->idef.functype[interaction->iatoms[i]]].nratoms+1) |
| 1324 | { |
| 1325 | /* next function */ |
| 1326 | if (func_type != top->idef.functype[interaction->iatoms[i]]) |
| 1327 | { |
| 1328 | gmx_incons("function type in search_donors")_gmx_error("incons", "function type in search_donors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1328); |
| 1329 | } |
| 1330 | |
| 1331 | if (interaction_function[func_type].flags & IF_VSITE1<<1) |
| 1332 | { |
| 1333 | nr1 = interaction->iatoms[i+1]; |
| 1334 | if (*top->atoms.atomname[nr1][0] == 'H') |
| 1335 | { |
| 1336 | nr2 = nr1-1; |
| 1337 | stop = FALSE0; |
| 1338 | while (!stop && ( *top->atoms.atomname[nr2][0] == 'H')) |
| 1339 | { |
| 1340 | if (nr2) |
| 1341 | { |
| 1342 | nr2--; |
| 1343 | } |
| 1344 | else |
| 1345 | { |
| 1346 | stop = TRUE1; |
| 1347 | } |
| 1348 | } |
| 1349 | if (!stop && ( ( *top->atoms.atomname[nr2][0] == 'O') || |
| 1350 | ( *top->atoms.atomname[nr2][0] == 'N') ) && |
| 1351 | ISINGRP(datable[nr1])(((datable[nr1]) & 4) == 4) && ISINGRP(datable[nr2])(((datable[nr2]) & 4) == 4)) |
| 1352 | { |
| 1353 | datable[nr2] |= DON; |
| 1354 | add_dh(ddd, nr2, nr1, grp, datable); |
| 1355 | } |
| 1356 | } |
| 1357 | } |
| 1358 | } |
| 1359 | } |
| 1360 | #endif |
| 1361 | } |
| 1362 | } |
| 1363 | |
| 1364 | static t_gridcell ***init_grid(gmx_bool bBox, rvec box[], real rcut, ivec ngrid) |
| 1365 | { |
| 1366 | t_gridcell ***grid; |
| 1367 | int i, y, z; |
| 1368 | |
| 1369 | if (bBox) |
| 1370 | { |
| 1371 | for (i = 0; i < DIM3; i++) |
| 1372 | { |
| 1373 | ngrid[i] = (box[i][i]/(1.2*rcut)); |
| 1374 | } |
| 1375 | } |
| 1376 | |
| 1377 | if (!bBox || (ngrid[XX0] < 3) || (ngrid[YY1] < 3) || (ngrid[ZZ2] < 3) ) |
| 1378 | { |
| 1379 | for (i = 0; i < DIM3; i++) |
| 1380 | { |
| 1381 | ngrid[i] = 1; |
| 1382 | } |
| 1383 | } |
| 1384 | else |
| 1385 | { |
| 1386 | printf("\nWill do grid-seach on %dx%dx%d grid, rcut=%g\n", |
| 1387 | ngrid[XX0], ngrid[YY1], ngrid[ZZ2], rcut); |
| 1388 | } |
| 1389 | snew(grid, ngrid[ZZ])(grid) = save_calloc("grid", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1389, (ngrid[2]), sizeof(*(grid))); |
| 1390 | for (z = 0; z < ngrid[ZZ2]; z++) |
| 1391 | { |
| 1392 | snew((grid)[z], ngrid[YY])((grid)[z]) = save_calloc("(grid)[z]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1392, (ngrid[1]), sizeof(*((grid)[z]))); |
| 1393 | for (y = 0; y < ngrid[YY1]; y++) |
| 1394 | { |
| 1395 | snew((grid)[z][y], ngrid[XX])((grid)[z][y]) = save_calloc("(grid)[z][y]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1395, (ngrid[0]), sizeof(*((grid)[z][y]))); |
| 1396 | } |
| 1397 | } |
| 1398 | return grid; |
| 1399 | } |
| 1400 | |
| 1401 | static void reset_nhbonds(t_donors *ddd) |
| 1402 | { |
| 1403 | int i, j; |
| 1404 | |
| 1405 | for (i = 0; (i < ddd->nrd); i++) |
| 1406 | { |
| 1407 | for (j = 0; (j < MAXHH4); j++) |
| 1408 | { |
| 1409 | ddd->nhbonds[i][j] = 0; |
| 1410 | } |
| 1411 | } |
| 1412 | } |
| 1413 | |
| 1414 | void pbc_correct_gem(rvec dx, matrix box, rvec hbox); |
| 1415 | |
| 1416 | static void build_grid(t_hbdata *hb, rvec x[], rvec xshell, |
| 1417 | gmx_bool bBox, matrix box, rvec hbox, |
| 1418 | real rcut, real rshell, |
| 1419 | ivec ngrid, t_gridcell ***grid) |
| 1420 | { |
| 1421 | int i, m, gr, xi, yi, zi, nr; |
| 1422 | atom_id *ad; |
| 1423 | ivec grididx; |
| 1424 | rvec invdelta, dshell, xtemp = {0, 0, 0}; |
| 1425 | t_ncell *newgrid; |
| 1426 | gmx_bool bDoRshell, bInShell, bAcc; |
| 1427 | real rshell2 = 0; |
| 1428 | int gx, gy, gz; |
| 1429 | int dum = -1; |
| 1430 | |
| 1431 | bDoRshell = (rshell > 0); |
| 1432 | rshell2 = sqr(rshell); |
| 1433 | bInShell = TRUE1; |
| 1434 | |
| 1435 | #define DBB(x)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1435,"x", x) if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n", __LINE__1435,#x, x) |
| 1436 | DBB(dum)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1436,"dum", dum); |
| 1437 | for (m = 0; m < DIM3; m++) |
| 1438 | { |
| 1439 | hbox[m] = box[m][m]*0.5; |
| 1440 | if (bBox) |
| 1441 | { |
| 1442 | invdelta[m] = ngrid[m]/box[m][m]; |
| 1443 | if (1/invdelta[m] < rcut) |
| 1444 | { |
| 1445 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1445, "Your computational box has shrunk too much.\n" |
| 1446 | "%s can not handle this situation, sorry.\n", |
| 1447 | ShortProgram()); |
| 1448 | } |
| 1449 | } |
| 1450 | else |
| 1451 | { |
| 1452 | invdelta[m] = 0; |
| 1453 | } |
| 1454 | } |
| 1455 | grididx[XX0] = 0; |
| 1456 | grididx[YY1] = 0; |
| 1457 | grididx[ZZ2] = 0; |
| 1458 | DBB(dum)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1458,"dum", dum); |
| 1459 | /* resetting atom counts */ |
| 1460 | for (gr = 0; (gr < grNR); gr++) |
| 1461 | { |
| 1462 | for (zi = 0; zi < ngrid[ZZ2]; zi++) |
| 1463 | { |
| 1464 | for (yi = 0; yi < ngrid[YY1]; yi++) |
| 1465 | { |
| 1466 | for (xi = 0; xi < ngrid[XX0]; xi++) |
| 1467 | { |
| 1468 | grid[zi][yi][xi].d[gr].nr = 0; |
| 1469 | grid[zi][yi][xi].a[gr].nr = 0; |
| 1470 | } |
| 1471 | } |
| 1472 | } |
| 1473 | DBB(dum)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1473,"dum", dum); |
| 1474 | |
| 1475 | /* put atoms in grid cells */ |
| 1476 | for (bAcc = FALSE0; (bAcc <= TRUE1); bAcc++) |
| 1477 | { |
| 1478 | if (bAcc) |
| 1479 | { |
| 1480 | nr = hb->a.nra; |
| 1481 | ad = hb->a.acc; |
| 1482 | } |
| 1483 | else |
| 1484 | { |
| 1485 | nr = hb->d.nrd; |
| 1486 | ad = hb->d.don; |
| 1487 | } |
| 1488 | DBB(bAcc)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1488,"bAcc", bAcc); |
| 1489 | for (i = 0; (i < nr); i++) |
| 1490 | { |
| 1491 | /* check if we are inside the shell */ |
| 1492 | /* if bDoRshell=FALSE then bInShell=TRUE always */ |
| 1493 | DBB(i)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1493,"i", i); |
| 1494 | if (bDoRshell) |
| 1495 | { |
| 1496 | bInShell = TRUE1; |
| 1497 | rvec_sub(x[ad[i]], xshell, dshell); |
| 1498 | if (bBox) |
| 1499 | { |
| 1500 | if (FALSE0 && !hb->bGem) |
| 1501 | { |
| 1502 | for (m = DIM3-1; m >= 0 && bInShell; m--) |
| 1503 | { |
| 1504 | if (dshell[m] < -hbox[m]) |
| 1505 | { |
| 1506 | rvec_inc(dshell, box[m]); |
| 1507 | } |
| 1508 | else if (dshell[m] >= hbox[m]) |
| 1509 | { |
| 1510 | dshell[m] -= 2*hbox[m]; |
| 1511 | } |
| 1512 | /* if we're outside the cube, we're outside the sphere also! */ |
| 1513 | if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) |
| 1514 | { |
| 1515 | bInShell = FALSE0; |
| 1516 | } |
| 1517 | } |
| 1518 | } |
| 1519 | else |
| 1520 | { |
| 1521 | gmx_bool bDone = FALSE0; |
| 1522 | while (!bDone) |
| 1523 | { |
| 1524 | bDone = TRUE1; |
| 1525 | for (m = DIM3-1; m >= 0 && bInShell; m--) |
| 1526 | { |
| 1527 | if (dshell[m] < -hbox[m]) |
| 1528 | { |
| 1529 | bDone = FALSE0; |
| 1530 | rvec_inc(dshell, box[m]); |
| 1531 | } |
| 1532 | if (dshell[m] >= hbox[m]) |
| 1533 | { |
| 1534 | bDone = FALSE0; |
| 1535 | dshell[m] -= 2*hbox[m]; |
| 1536 | } |
| 1537 | } |
| 1538 | } |
| 1539 | for (m = DIM3-1; m >= 0 && bInShell; m--) |
| 1540 | { |
| 1541 | /* if we're outside the cube, we're outside the sphere also! */ |
| 1542 | if ( (dshell[m] > rshell) || (-dshell[m] > rshell) ) |
| 1543 | { |
| 1544 | bInShell = FALSE0; |
| 1545 | } |
| 1546 | } |
| 1547 | } |
| 1548 | } |
| 1549 | /* if we're inside the cube, check if we're inside the sphere */ |
| 1550 | if (bInShell) |
| 1551 | { |
| 1552 | bInShell = norm2(dshell) < rshell2; |
| 1553 | } |
| 1554 | } |
| 1555 | DBB(i)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1555,"i", i); |
| 1556 | if (bInShell) |
| 1557 | { |
| 1558 | if (bBox) |
| 1559 | { |
| 1560 | if (hb->bGem) |
| 1561 | { |
| 1562 | copy_rvec(x[ad[i]], xtemp); |
| 1563 | } |
| 1564 | pbc_correct_gem(x[ad[i]], box, hbox); |
| 1565 | } |
| 1566 | for (m = DIM3-1; m >= 0; m--) |
| 1567 | { |
| 1568 | if (TRUE1 || !hb->bGem) |
| 1569 | { |
| 1570 | /* put atom in the box */ |
| 1571 | while (x[ad[i]][m] < 0) |
| 1572 | { |
| 1573 | rvec_inc(x[ad[i]], box[m]); |
| 1574 | } |
| 1575 | while (x[ad[i]][m] >= box[m][m]) |
| 1576 | { |
| 1577 | rvec_dec(x[ad[i]], box[m]); |
| 1578 | } |
| 1579 | } |
| 1580 | /* determine grid index of atom */ |
| 1581 | grididx[m] = x[ad[i]][m]*invdelta[m]; |
| 1582 | grididx[m] = (grididx[m]+ngrid[m]) % ngrid[m]; |
| 1583 | } |
| 1584 | if (hb->bGem) |
| 1585 | { |
| 1586 | copy_rvec(xtemp, x[ad[i]]); /* copy back */ |
| 1587 | } |
| 1588 | gx = grididx[XX0]; |
| 1589 | gy = grididx[YY1]; |
| 1590 | gz = grididx[ZZ2]; |
| 1591 | range_check(gx, 0, ngrid[XX])_range_check(gx, 0, ngrid[0], ((void*)0),"gx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1591); |
| 1592 | range_check(gy, 0, ngrid[YY])_range_check(gy, 0, ngrid[1], ((void*)0),"gy", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1592); |
| 1593 | range_check(gz, 0, ngrid[ZZ])_range_check(gz, 0, ngrid[2], ((void*)0),"gz", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1593); |
| 1594 | DBB(gx)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1594,"gx", gx); |
| 1595 | DBB(gy)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1595,"gy", gy); |
| 1596 | DBB(gz)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1596,"gz", gz); |
| 1597 | /* add atom to grid cell */ |
| 1598 | if (bAcc) |
| 1599 | { |
| 1600 | newgrid = &(grid[gz][gy][gx].a[gr]); |
| 1601 | } |
| 1602 | else |
| 1603 | { |
| 1604 | newgrid = &(grid[gz][gy][gx].d[gr]); |
| 1605 | } |
| 1606 | if (newgrid->nr >= newgrid->maxnr) |
| 1607 | { |
| 1608 | newgrid->maxnr += 10; |
| 1609 | DBB(newgrid->maxnr)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1609,"newgrid->maxnr", newgrid->maxnr); |
| 1610 | srenew(newgrid->atoms, newgrid->maxnr)(newgrid->atoms) = save_realloc("newgrid->atoms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1610, (newgrid->atoms), (newgrid->maxnr), sizeof(*(newgrid ->atoms))); |
| 1611 | } |
| 1612 | DBB(newgrid->nr)if (debug && bDebug) fprintf(debug, "build_grid, line %d, %s = %d\n" , 1612,"newgrid->nr", newgrid->nr); |
| 1613 | newgrid->atoms[newgrid->nr] = ad[i]; |
| 1614 | newgrid->nr++; |
| 1615 | } |
| 1616 | } |
| 1617 | } |
| 1618 | } |
| 1619 | } |
| 1620 | |
| 1621 | static void count_da_grid(ivec ngrid, t_gridcell ***grid, t_icell danr) |
| 1622 | { |
| 1623 | int gr, xi, yi, zi; |
| 1624 | |
| 1625 | for (gr = 0; (gr < grNR); gr++) |
| 1626 | { |
| 1627 | danr[gr] = 0; |
| 1628 | for (zi = 0; zi < ngrid[ZZ2]; zi++) |
| 1629 | { |
| 1630 | for (yi = 0; yi < ngrid[YY1]; yi++) |
| 1631 | { |
| 1632 | for (xi = 0; xi < ngrid[XX0]; xi++) |
| 1633 | { |
| 1634 | danr[gr] += grid[zi][yi][xi].d[gr].nr; |
| 1635 | } |
| 1636 | } |
| 1637 | } |
| 1638 | } |
| 1639 | } |
| 1640 | |
| 1641 | /* The grid loop. |
| 1642 | * Without a box, the grid is 1x1x1, so all loops are 1 long. |
| 1643 | * With a rectangular box (bTric==FALSE) all loops are 3 long. |
| 1644 | * With a triclinic box all loops are 3 long, except when a cell is |
| 1645 | * located next to one of the box edges which is not parallel to the |
| 1646 | * x/y-plane, in that case all cells in a line or layer are searched. |
| 1647 | * This could be implemented slightly more efficient, but the code |
| 1648 | * would get much more complicated. |
| 1649 | */ |
| 1650 | static gmx_inlineinline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge) |
| 1651 | { |
| 1652 | return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1)); |
| 1653 | } |
| 1654 | static gmx_inlineinline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge) |
| 1655 | { |
| 1656 | return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1)); |
| 1657 | } |
| 1658 | static gmx_inlineinline int grid_mod(int j, int n) |
| 1659 | { |
| 1660 | return (j+n) % (n); |
| 1661 | } |
| 1662 | |
| 1663 | static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid) |
| 1664 | { |
| 1665 | int gr, x, y, z, sum[grNR]; |
| 1666 | |
| 1667 | fprintf(fp, "grid %dx%dx%d\n", ngrid[XX0], ngrid[YY1], ngrid[ZZ2]); |
| 1668 | for (gr = 0; gr < grNR; gr++) |
| 1669 | { |
| 1670 | sum[gr] = 0; |
| 1671 | fprintf(fp, "GROUP %d (%s)\n", gr, grpnames[gr]); |
| 1672 | for (z = 0; z < ngrid[ZZ2]; z += 2) |
| 1673 | { |
| 1674 | fprintf(fp, "Z=%d,%d\n", z, z+1); |
| 1675 | for (y = 0; y < ngrid[YY1]; y++) |
| 1676 | { |
| 1677 | for (x = 0; x < ngrid[XX0]; x++) |
| 1678 | { |
| 1679 | fprintf(fp, "%3d", grid[x][y][z].d[gr].nr); |
| 1680 | sum[gr] += grid[z][y][x].d[gr].nr; |
| 1681 | fprintf(fp, "%3d", grid[x][y][z].a[gr].nr); |
| 1682 | sum[gr] += grid[z][y][x].a[gr].nr; |
| 1683 | |
| 1684 | } |
| 1685 | fprintf(fp, " | "); |
| 1686 | if ( (z+1) < ngrid[ZZ2]) |
| 1687 | { |
| 1688 | for (x = 0; x < ngrid[XX0]; x++) |
| 1689 | { |
| 1690 | fprintf(fp, "%3d", grid[z+1][y][x].d[gr].nr); |
| 1691 | sum[gr] += grid[z+1][y][x].d[gr].nr; |
| 1692 | fprintf(fp, "%3d", grid[z+1][y][x].a[gr].nr); |
| 1693 | sum[gr] += grid[z+1][y][x].a[gr].nr; |
| 1694 | } |
| 1695 | } |
| 1696 | fprintf(fp, "\n"); |
| 1697 | } |
| 1698 | } |
| 1699 | } |
| 1700 | fprintf(fp, "TOTALS:"); |
| 1701 | for (gr = 0; gr < grNR; gr++) |
| 1702 | { |
| 1703 | fprintf(fp, " %d=%d", gr, sum[gr]); |
| 1704 | } |
| 1705 | fprintf(fp, "\n"); |
| 1706 | } |
| 1707 | |
| 1708 | /* New GMX record! 5 * in a row. Congratulations! |
| 1709 | * Sorry, only four left. |
| 1710 | */ |
| 1711 | static void free_grid(ivec ngrid, t_gridcell ****grid) |
| 1712 | { |
| 1713 | int y, z; |
| 1714 | t_gridcell ***g = *grid; |
| 1715 | |
| 1716 | for (z = 0; z < ngrid[ZZ2]; z++) |
| 1717 | { |
| 1718 | for (y = 0; y < ngrid[YY1]; y++) |
| 1719 | { |
| 1720 | sfree(g[z][y])save_free("g[z][y]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1720, (g[z][y])); |
| 1721 | } |
| 1722 | sfree(g[z])save_free("g[z]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1722, (g[z])); |
| 1723 | } |
| 1724 | sfree(g)save_free("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1724, (g)); |
| 1725 | g = NULL((void*)0); |
| 1726 | } |
| 1727 | |
| 1728 | void pbc_correct_gem(rvec dx, matrix box, rvec hbox) |
| 1729 | { |
| 1730 | int m; |
| 1731 | gmx_bool bDone = FALSE0; |
| 1732 | while (!bDone) |
| 1733 | { |
| 1734 | bDone = TRUE1; |
| 1735 | for (m = DIM3-1; m >= 0; m--) |
| 1736 | { |
| 1737 | if (dx[m] < -hbox[m]) |
| 1738 | { |
| 1739 | bDone = FALSE0; |
| 1740 | rvec_inc(dx, box[m]); |
| 1741 | } |
| 1742 | if (dx[m] >= hbox[m]) |
| 1743 | { |
| 1744 | bDone = FALSE0; |
| 1745 | rvec_dec(dx, box[m]); |
| 1746 | } |
| 1747 | } |
| 1748 | } |
| 1749 | } |
| 1750 | |
| 1751 | /* Added argument r2cut, changed contact and implemented |
| 1752 | * use of second cut-off. |
| 1753 | * - Erik Marklund, June 29, 2006 |
| 1754 | */ |
| 1755 | static int is_hbond(t_hbdata *hb, int grpd, int grpa, int d, int a, |
| 1756 | real rcut, real r2cut, real ccut, |
| 1757 | rvec x[], gmx_bool bBox, matrix box, rvec hbox, |
| 1758 | real *d_ha, real *ang, gmx_bool bDA, int *hhh, |
| 1759 | gmx_bool bContact, gmx_bool bMerge, PSTYPEint *p) |
| 1760 | { |
| 1761 | int h, hh, id, ja, ihb; |
| 1762 | rvec r_da, r_ha, r_dh, r = {0, 0, 0}; |
| 1763 | ivec ri; |
| 1764 | real rc2, r2c2, rda2, rha2, ca; |
| 1765 | gmx_bool HAinrange = FALSE0; /* If !bDA. Needed for returning hbDist in a correct way. */ |
| 1766 | gmx_bool daSwap = FALSE0; |
| 1767 | |
| 1768 | if (d == a) |
| 1769 | { |
| 1770 | return hbNo; |
| 1771 | } |
| 1772 | |
| 1773 | if (((id = donor_index(&hb->d, grpd, d)_donor_index(&hb->d, grpd, d, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1773)) == NOTSET-12345) || |
| 1774 | ((ja = acceptor_index(&hb->a, grpa, a)_acceptor_index(&hb->a, grpa, a, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1774)) == NOTSET-12345)) |
| 1775 | { |
| 1776 | return hbNo; |
| 1777 | } |
| 1778 | |
| 1779 | rc2 = rcut*rcut; |
| 1780 | r2c2 = r2cut*r2cut; |
| 1781 | |
| 1782 | rvec_sub(x[d], x[a], r_da); |
| 1783 | /* Insert projection code here */ |
| 1784 | |
| 1785 | if (bMerge && d > a && isInterchangable(hb, d, a, grpd, grpa)) |
| 1786 | { |
| 1787 | /* Then this hbond/contact will be found again, or it has already been found. */ |
| 1788 | /*return hbNo;*/ |
| 1789 | } |
| 1790 | if (bBox) |
| 1791 | { |
| 1792 | if (d > a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) /* acceptor is also a donor and vice versa? */ |
| 1793 | { /* return hbNo; */ |
| 1794 | daSwap = TRUE1; /* If so, then their history should be filed with donor and acceptor swapped. */ |
| 1795 | } |
| 1796 | if (hb->bGem) |
| 1797 | { |
| 1798 | copy_rvec(r_da, r); /* Save this for later */ |
| 1799 | pbc_correct_gem(r_da, box, hbox); |
| 1800 | } |
| 1801 | else |
| 1802 | { |
| 1803 | pbc_correct_gem(r_da, box, hbox); |
| 1804 | } |
| 1805 | } |
| 1806 | rda2 = iprod(r_da, r_da); |
| 1807 | |
| 1808 | if (bContact) |
| 1809 | { |
| 1810 | if (daSwap && grpa == grpd) |
| 1811 | { |
| 1812 | return hbNo; |
| 1813 | } |
| 1814 | if (rda2 <= rc2) |
| 1815 | { |
| 1816 | if (hb->bGem) |
| 1817 | { |
| 1818 | calcBoxDistance(hb->per->P, r, ri); |
| 1819 | *p = periodicIndex(ri, hb->per, daSwap); /* find (or add) periodicity index. */ |
| 1820 | } |
| 1821 | return hbHB; |
| 1822 | } |
| 1823 | else if (rda2 < r2c2) |
| 1824 | { |
| 1825 | return hbDist; |
| 1826 | } |
| 1827 | else |
| 1828 | { |
| 1829 | return hbNo; |
| 1830 | } |
| 1831 | } |
| 1832 | *hhh = NOTSET-12345; |
| 1833 | |
| 1834 | if (bDA && (rda2 > rc2)) |
| 1835 | { |
| 1836 | return hbNo; |
| 1837 | } |
| 1838 | |
| 1839 | for (h = 0; (h < hb->d.nhydro[id]); h++) |
| 1840 | { |
| 1841 | hh = hb->d.hydro[id][h]; |
| 1842 | rha2 = rc2+1; |
| 1843 | if (!bDA) |
| 1844 | { |
| 1845 | rvec_sub(x[hh], x[a], r_ha); |
| 1846 | if (bBox) |
| 1847 | { |
| 1848 | pbc_correct_gem(r_ha, box, hbox); |
| 1849 | } |
| 1850 | rha2 = iprod(r_ha, r_ha); |
| 1851 | } |
| 1852 | |
| 1853 | if (hb->bGem) |
| 1854 | { |
| 1855 | calcBoxDistance(hb->per->P, r, ri); |
| 1856 | *p = periodicIndex(ri, hb->per, daSwap); /* find periodicity index. */ |
| 1857 | } |
| 1858 | |
| 1859 | if (bDA || (!bDA && (rha2 <= rc2))) |
| 1860 | { |
| 1861 | rvec_sub(x[d], x[hh], r_dh); |
| 1862 | if (bBox) |
| 1863 | { |
| 1864 | pbc_correct_gem(r_dh, box, hbox); |
| 1865 | } |
| 1866 | |
| 1867 | if (!bDA) |
| 1868 | { |
| 1869 | HAinrange = TRUE1; |
| 1870 | } |
| 1871 | ca = cos_angle(r_dh, r_da); |
| 1872 | /* if angle is smaller, cos is larger */ |
| 1873 | if (ca >= ccut) |
| 1874 | { |
| 1875 | *hhh = hh; |
| 1876 | *d_ha = sqrt(bDA ? rda2 : rha2); |
| 1877 | *ang = acos(ca); |
| 1878 | return hbHB; |
| 1879 | } |
| 1880 | } |
| 1881 | } |
| 1882 | if (bDA || (!bDA && HAinrange)) |
| 1883 | { |
| 1884 | return hbDist; |
| 1885 | } |
| 1886 | else |
| 1887 | { |
| 1888 | return hbNo; |
| 1889 | } |
| 1890 | } |
| 1891 | |
| 1892 | /* Fixed previously undiscovered bug in the merge |
| 1893 | code, where the last frame of each hbond disappears. |
| 1894 | - Erik Marklund, June 1, 2006 */ |
| 1895 | /* Added the following arguments: |
| 1896 | * ptmp[] - temporary periodicity hisory |
| 1897 | * a1 - identity of first acceptor/donor |
| 1898 | * a2 - identity of second acceptor/donor |
| 1899 | * - Erik Marklund, FEB 20 2010 */ |
| 1900 | |
| 1901 | /* Merging is now done on the fly, so do_merge is most likely obsolete now. |
| 1902 | * Will do some more testing before removing the function entirely. |
| 1903 | * - Erik Marklund, MAY 10 2010 */ |
| 1904 | static void do_merge(t_hbdata *hb, int ntmp, |
| 1905 | unsigned int htmp[], unsigned int gtmp[], PSTYPEint ptmp[], |
| 1906 | t_hbond *hb0, t_hbond *hb1, int a1, int a2) |
| 1907 | { |
| 1908 | /* Here we need to make sure we're treating periodicity in |
| 1909 | * the right way for the geminate recombination kinetics. */ |
| 1910 | |
| 1911 | int m, mm, n00, n01, nn0, nnframes; |
| 1912 | PSTYPEint pm; |
| 1913 | t_pShift *pShift; |
| 1914 | |
| 1915 | /* Decide where to start from when merging */ |
| 1916 | n00 = hb0->n0; |
| 1917 | n01 = hb1->n0; |
| 1918 | nn0 = min(n00, n01)(((n00) < (n01)) ? (n00) : (n01) ); |
| 1919 | nnframes = max(n00 + hb0->nframes, n01 + hb1->nframes)(((n00 + hb0->nframes) > (n01 + hb1->nframes)) ? (n00 + hb0->nframes) : (n01 + hb1->nframes) ) - nn0; |
| 1920 | /* Initiate tmp arrays */ |
| 1921 | for (m = 0; (m < ntmp); m++) |
| 1922 | { |
| 1923 | htmp[m] = 0; |
| 1924 | gtmp[m] = 0; |
| 1925 | ptmp[m] = 0; |
| 1926 | } |
| 1927 | /* Fill tmp arrays with values due to first HB */ |
| 1928 | /* Once again '<' had to be replaced with '<=' |
| 1929 | to catch the last frame in which the hbond |
| 1930 | appears. |
| 1931 | - Erik Marklund, June 1, 2006 */ |
| 1932 | for (m = 0; (m <= hb0->nframes); m++) |
| 1933 | { |
| 1934 | mm = m+n00-nn0; |
| 1935 | htmp[mm] = is_hb(hb0->h[0], m); |
| 1936 | if (hb->bGem) |
| 1937 | { |
| 1938 | pm = getPshift(hb->per->pHist[a1][a2], m+hb0->n0); |
| 1939 | if (pm > hb->per->nper) |
| 1940 | { |
| 1941 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1941, "Illegal shift!"); |
| 1942 | } |
| 1943 | else |
| 1944 | { |
| 1945 | ptmp[mm] = pm; /*hb->per->pHist[a1][a2][m];*/ |
| 1946 | } |
| 1947 | } |
| 1948 | } |
| 1949 | /* If we're doing geminate recompbination we usually don't need the distances. |
| 1950 | * Let's save some memory and time. */ |
| 1951 | if (TRUE1 || !hb->bGem || hb->per->gemtype == gemAD) |
| 1952 | { |
| 1953 | for (m = 0; (m <= hb0->nframes); m++) |
| 1954 | { |
| 1955 | mm = m+n00-nn0; |
| 1956 | gtmp[mm] = is_hb(hb0->g[0], m); |
| 1957 | } |
| 1958 | } |
| 1959 | /* Next HB */ |
| 1960 | for (m = 0; (m <= hb1->nframes); m++) |
| 1961 | { |
| 1962 | mm = m+n01-nn0; |
| 1963 | htmp[mm] = htmp[mm] || is_hb(hb1->h[0], m); |
| 1964 | gtmp[mm] = gtmp[mm] || is_hb(hb1->g[0], m); |
| 1965 | if (hb->bGem /* && ptmp[mm] != 0 */) |
| 1966 | { |
| 1967 | |
| 1968 | /* If this hbond has been seen before with donor and acceptor swapped, |
| 1969 | * then we need to find the mirrored (*-1) periodicity vector to truely |
| 1970 | * merge the hbond history. */ |
| 1971 | pm = findMirror(getPshift(hb->per->pHist[a2][a1], m+hb1->n0), hb->per->p2i, hb->per->nper); |
| 1972 | /* Store index of mirror */ |
| 1973 | if (pm > hb->per->nper) |
| 1974 | { |
| 1975 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1975, "Illegal shift!"); |
| 1976 | } |
| 1977 | ptmp[mm] = pm; |
| 1978 | } |
| 1979 | } |
| 1980 | /* Reallocate target array */ |
| 1981 | if (nnframes > hb0->maxframes) |
| 1982 | { |
| 1983 | srenew(hb0->h[0], 4+nnframes/hb->wordlen)(hb0->h[0]) = save_realloc("hb0->h[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1983, (hb0->h[0]), (4+nnframes/hb->wordlen), sizeof(* (hb0->h[0]))); |
| 1984 | srenew(hb0->g[0], 4+nnframes/hb->wordlen)(hb0->g[0]) = save_realloc("hb0->g[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 1984, (hb0->g[0]), (4+nnframes/hb->wordlen), sizeof(* (hb0->g[0]))); |
| 1985 | } |
| 1986 | if (NULL((void*)0) != hb->per->pHist) |
| 1987 | { |
| 1988 | clearPshift(&(hb->per->pHist[a1][a2])); |
| 1989 | } |
| 1990 | |
| 1991 | /* Copy temp array to target array */ |
| 1992 | for (m = 0; (m <= nnframes); m++) |
| 1993 | { |
| 1994 | _set_hb(hb0->h[0], m, htmp[m]); |
| 1995 | _set_hb(hb0->g[0], m, gtmp[m]); |
| 1996 | if (hb->bGem) |
| 1997 | { |
| 1998 | addPshift(&(hb->per->pHist[a1][a2]), ptmp[m], m+nn0); |
| 1999 | } |
| 2000 | } |
| 2001 | |
| 2002 | /* Set scalar variables */ |
| 2003 | hb0->n0 = nn0; |
| 2004 | hb0->maxframes = nnframes; |
| 2005 | } |
| 2006 | |
| 2007 | /* Added argument bContact for nicer output. |
| 2008 | * Erik Marklund, June 29, 2006 |
| 2009 | */ |
| 2010 | static void merge_hb(t_hbdata *hb, gmx_bool bTwo, gmx_bool bContact) |
| 2011 | { |
| 2012 | int i, inrnew, indnew, j, ii, jj, m, id, ia, grp, ogrp, ntmp; |
| 2013 | unsigned int *htmp, *gtmp; |
| 2014 | PSTYPEint *ptmp; |
| 2015 | t_hbond *hb0, *hb1; |
| 2016 | |
| 2017 | inrnew = hb->nrhb; |
| 2018 | indnew = hb->nrdist; |
| 2019 | |
| 2020 | /* Check whether donors are also acceptors */ |
| 2021 | printf("Merging hbonds with Acceptor and Donor swapped\n"); |
| 2022 | |
| 2023 | ntmp = 2*hb->max_frames; |
| 2024 | snew(gtmp, ntmp)(gtmp) = save_calloc("gtmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2024, (ntmp), sizeof(*(gtmp))); |
| 2025 | snew(htmp, ntmp)(htmp) = save_calloc("htmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2025, (ntmp), sizeof(*(htmp))); |
| 2026 | snew(ptmp, ntmp)(ptmp) = save_calloc("ptmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2026, (ntmp), sizeof(*(ptmp))); |
| 2027 | for (i = 0; (i < hb->d.nrd); i++) |
| 2028 | { |
| 2029 | fprintf(stderrstderr, "\r%d/%d", i+1, hb->d.nrd); |
| 2030 | id = hb->d.don[i]; |
| 2031 | ii = hb->a.aptr[id]; |
| 2032 | for (j = 0; (j < hb->a.nra); j++) |
| 2033 | { |
| 2034 | ia = hb->a.acc[j]; |
| 2035 | jj = hb->d.dptr[ia]; |
| 2036 | if ((id != ia) && (ii != NOTSET-12345) && (jj != NOTSET-12345) && |
| 2037 | (!bTwo || (bTwo && (hb->d.grp[i] != hb->a.grp[j])))) |
| 2038 | { |
| 2039 | hb0 = hb->hbmap[i][j]; |
| 2040 | hb1 = hb->hbmap[jj][ii]; |
| 2041 | if (hb0 && hb1 && ISHB(hb0->history[0])(((hb0->history[0]) & 2) == 2) && ISHB(hb1->history[0])(((hb1->history[0]) & 2) == 2)) |
| 2042 | { |
| 2043 | do_merge(hb, ntmp, htmp, gtmp, ptmp, hb0, hb1, i, j); |
| 2044 | if (ISHB(hb1->history[0])(((hb1->history[0]) & 2) == 2)) |
| 2045 | { |
| 2046 | inrnew--; |
| 2047 | } |
| 2048 | else if (ISDIST(hb1->history[0])(((hb1->history[0]) & 1) == 1)) |
| 2049 | { |
| 2050 | indnew--; |
| 2051 | } |
| 2052 | else |
| 2053 | if (bContact) |
| 2054 | { |
| 2055 | gmx_incons("No contact history")_gmx_error("incons", "No contact history", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2055); |
| 2056 | } |
| 2057 | else |
| 2058 | { |
| 2059 | gmx_incons("Neither hydrogen bond nor distance")_gmx_error("incons", "Neither hydrogen bond nor distance", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2059); |
| 2060 | } |
| 2061 | sfree(hb1->h[0])save_free("hb1->h[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2061, (hb1->h[0])); |
| 2062 | sfree(hb1->g[0])save_free("hb1->g[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2062, (hb1->g[0])); |
| 2063 | if (hb->bGem) |
| 2064 | { |
| 2065 | clearPshift(&(hb->per->pHist[jj][ii])); |
| 2066 | } |
| 2067 | hb1->h[0] = NULL((void*)0); |
| 2068 | hb1->g[0] = NULL((void*)0); |
| 2069 | hb1->history[0] = hbNo; |
| 2070 | } |
| 2071 | } |
| 2072 | } |
| 2073 | } |
| 2074 | fprintf(stderrstderr, "\n"); |
| 2075 | printf("- Reduced number of hbonds from %d to %d\n", hb->nrhb, inrnew); |
| 2076 | printf("- Reduced number of distances from %d to %d\n", hb->nrdist, indnew); |
| 2077 | hb->nrhb = inrnew; |
| 2078 | hb->nrdist = indnew; |
| 2079 | sfree(gtmp)save_free("gtmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2079, (gtmp)); |
| 2080 | sfree(htmp)save_free("htmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2080, (htmp)); |
| 2081 | sfree(ptmp)save_free("ptmp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2081, (ptmp)); |
| 2082 | } |
| 2083 | |
| 2084 | static void do_nhb_dist(FILE *fp, t_hbdata *hb, real t) |
| 2085 | { |
| 2086 | int i, j, k, n_bound[MAXHH4], nbtot; |
| 2087 | h_id nhb; |
| 2088 | |
| 2089 | |
| 2090 | /* Set array to 0 */ |
| 2091 | for (k = 0; (k < MAXHH4); k++) |
| 2092 | { |
| 2093 | n_bound[k] = 0; |
| 2094 | } |
| 2095 | /* Loop over possible donors */ |
| 2096 | for (i = 0; (i < hb->d.nrd); i++) |
| 2097 | { |
| 2098 | for (j = 0; (j < hb->d.nhydro[i]); j++) |
| 2099 | { |
| 2100 | n_bound[hb->d.nhbonds[i][j]]++; |
| 2101 | } |
| 2102 | } |
| 2103 | fprintf(fp, "%12.5e", t); |
| 2104 | nbtot = 0; |
| 2105 | for (k = 0; (k < MAXHH4); k++) |
| 2106 | { |
| 2107 | fprintf(fp, " %8d", n_bound[k]); |
| 2108 | nbtot += n_bound[k]*k; |
| 2109 | } |
| 2110 | fprintf(fp, " %8d\n", nbtot); |
| 2111 | } |
| 2112 | |
| 2113 | /* Added argument bContact in do_hblife(...). Also |
| 2114 | * added support for -contact in function body. |
| 2115 | * - Erik Marklund, May 31, 2006 */ |
| 2116 | /* Changed the contact code slightly. |
| 2117 | * - Erik Marklund, June 29, 2006 |
| 2118 | */ |
| 2119 | static void do_hblife(const char *fn, t_hbdata *hb, gmx_bool bMerge, gmx_bool bContact, |
| 2120 | const output_env_t oenv) |
| 2121 | { |
| 2122 | FILE *fp; |
| 2123 | const char *leg[] = { "p(t)", "t p(t)" }; |
| 2124 | int *histo; |
| 2125 | int i, j, j0, k, m, nh, ihb, ohb, nhydro, ndump = 0; |
| 2126 | int nframes = hb->nframes; |
| 2127 | unsigned int **h; |
| 2128 | real t, x1, dt; |
| 2129 | double sum, integral; |
| 2130 | t_hbond *hbh; |
| 2131 | |
| 2132 | snew(h, hb->maxhydro)(h) = save_calloc("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2132, (hb->maxhydro), sizeof(*(h))); |
| 2133 | snew(histo, nframes+1)(histo) = save_calloc("histo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2133, (nframes+1), sizeof(*(histo))); |
| 2134 | /* Total number of hbonds analyzed here */ |
| 2135 | for (i = 0; (i < hb->d.nrd); i++) |
| 2136 | { |
| 2137 | for (k = 0; (k < hb->a.nra); k++) |
| 2138 | { |
| 2139 | hbh = hb->hbmap[i][k]; |
| 2140 | if (hbh) |
| 2141 | { |
| 2142 | if (bMerge) |
| 2143 | { |
| 2144 | if (hbh->h[0]) |
| 2145 | { |
| 2146 | h[0] = hbh->h[0]; |
| 2147 | nhydro = 1; |
| 2148 | } |
| 2149 | else |
| 2150 | { |
| 2151 | nhydro = 0; |
| 2152 | } |
| 2153 | } |
| 2154 | else |
| 2155 | { |
| 2156 | nhydro = 0; |
| 2157 | for (m = 0; (m < hb->maxhydro); m++) |
| 2158 | { |
| 2159 | if (hbh->h[m]) |
| 2160 | { |
| 2161 | h[nhydro++] = bContact ? hbh->g[m] : hbh->h[m]; |
| 2162 | } |
| 2163 | } |
| 2164 | } |
| 2165 | for (nh = 0; (nh < nhydro); nh++) |
| 2166 | { |
| 2167 | ohb = 0; |
| 2168 | j0 = 0; |
| 2169 | |
| 2170 | /* Changed '<' into '<=' below, just like I |
| 2171 | did in the hbm-output-loop in the main code. |
| 2172 | - Erik Marklund, May 31, 2006 |
| 2173 | */ |
| 2174 | for (j = 0; (j <= hbh->nframes); j++) |
| 2175 | { |
| 2176 | ihb = is_hb(h[nh], j); |
| 2177 | if (debug && (ndump < 10)) |
| 2178 | { |
| 2179 | fprintf(debug, "%5d %5d\n", j, ihb); |
| 2180 | } |
| 2181 | if (ihb != ohb) |
| 2182 | { |
| 2183 | if (ihb) |
| 2184 | { |
| 2185 | j0 = j; |
| 2186 | } |
| 2187 | else |
| 2188 | { |
| 2189 | histo[j-j0]++; |
| 2190 | } |
| 2191 | ohb = ihb; |
| 2192 | } |
| 2193 | } |
| 2194 | ndump++; |
| 2195 | } |
| 2196 | } |
| 2197 | } |
| 2198 | } |
| 2199 | fprintf(stderrstderr, "\n"); |
| 2200 | if (bContact) |
| 2201 | { |
| 2202 | fp = xvgropen(fn, "Uninterrupted contact lifetime", output_env_get_xvgr_tlabel(oenv), "()", oenv); |
| 2203 | } |
| 2204 | else |
| 2205 | { |
| 2206 | fp = xvgropen(fn, "Uninterrupted hydrogen bond lifetime", output_env_get_xvgr_tlabel(oenv), "()", |
| 2207 | oenv); |
| 2208 | } |
| 2209 | |
| 2210 | xvgr_legend(fp, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); |
| 2211 | j0 = nframes-1; |
| 2212 | while ((j0 > 0) && (histo[j0] == 0)) |
| 2213 | { |
| 2214 | j0--; |
| 2215 | } |
| 2216 | sum = 0; |
| 2217 | for (i = 0; (i <= j0); i++) |
| 2218 | { |
| 2219 | sum += histo[i]; |
| 2220 | } |
| 2221 | dt = hb->time[1]-hb->time[0]; |
| 2222 | sum = dt*sum; |
| 2223 | integral = 0; |
| 2224 | for (i = 1; (i <= j0); i++) |
| 2225 | { |
| 2226 | t = hb->time[i] - hb->time[0] - 0.5*dt; |
| 2227 | x1 = t*histo[i]/sum; |
| 2228 | fprintf(fp, "%8.3f %10.3e %10.3e\n", t, histo[i]/sum, x1); |
| 2229 | integral += x1; |
| 2230 | } |
| 2231 | integral *= dt; |
| 2232 | gmx_ffclose(fp); |
| 2233 | printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral); |
| 2234 | printf("Note that the lifetime obtained in this manner is close to useless\n"); |
| 2235 | printf("Use the -ac option instead and check the Forward lifetime\n"); |
| 2236 | please_cite(stdoutstdout, "Spoel2006b"); |
| 2237 | sfree(h)save_free("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2237, (h)); |
| 2238 | sfree(histo)save_free("histo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2238, (histo)); |
| 2239 | } |
| 2240 | |
| 2241 | /* Changed argument bMerge into oneHB to handle contacts properly. |
| 2242 | * - Erik Marklund, June 29, 2006 |
| 2243 | */ |
| 2244 | static void dump_ac(t_hbdata *hb, gmx_bool oneHB, int nDump) |
| 2245 | { |
| 2246 | FILE *fp; |
| 2247 | int i, j, k, m, nd, ihb, idist; |
| 2248 | int nframes = hb->nframes; |
| 2249 | gmx_bool bPrint; |
| 2250 | t_hbond *hbh; |
| 2251 | |
| 2252 | if (nDump <= 0) |
| 2253 | { |
| 2254 | return; |
| 2255 | } |
| 2256 | fp = gmx_ffopen("debug-ac.xvg", "w"); |
| 2257 | for (j = 0; (j < nframes); j++) |
| 2258 | { |
| 2259 | fprintf(fp, "%10.3f", hb->time[j]); |
| 2260 | for (i = nd = 0; (i < hb->d.nrd) && (nd < nDump); i++) |
| 2261 | { |
| 2262 | for (k = 0; (k < hb->a.nra) && (nd < nDump); k++) |
| 2263 | { |
| 2264 | bPrint = FALSE0; |
| 2265 | ihb = idist = 0; |
| 2266 | hbh = hb->hbmap[i][k]; |
| 2267 | if (oneHB) |
| 2268 | { |
| 2269 | if (hbh->h[0]) |
| 2270 | { |
| 2271 | ihb = is_hb(hbh->h[0], j); |
| 2272 | idist = is_hb(hbh->g[0], j); |
| 2273 | bPrint = TRUE1; |
| 2274 | } |
| 2275 | } |
| 2276 | else |
| 2277 | { |
| 2278 | for (m = 0; (m < hb->maxhydro) && !ihb; m++) |
| 2279 | { |
| 2280 | ihb = ihb || ((hbh->h[m]) && is_hb(hbh->h[m], j)); |
| 2281 | idist = idist || ((hbh->g[m]) && is_hb(hbh->g[m], j)); |
| 2282 | } |
| 2283 | /* This is not correct! */ |
| 2284 | /* What isn't correct? -Erik M */ |
| 2285 | bPrint = TRUE1; |
| 2286 | } |
| 2287 | if (bPrint) |
| 2288 | { |
| 2289 | fprintf(fp, " %1d-%1d", ihb, idist); |
| 2290 | nd++; |
| 2291 | } |
| 2292 | } |
| 2293 | } |
| 2294 | fprintf(fp, "\n"); |
| 2295 | } |
| 2296 | gmx_ffclose(fp); |
| 2297 | } |
| 2298 | |
| 2299 | static real calc_dg(real tau, real temp) |
| 2300 | { |
| 2301 | real kbt; |
| 2302 | |
| 2303 | kbt = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp; |
| 2304 | if (tau <= 0) |
| 2305 | { |
| 2306 | return -666; |
| 2307 | } |
| 2308 | else |
| 2309 | { |
| 2310 | return kbt*log(kbt*tau/PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))); |
| 2311 | } |
| 2312 | } |
| 2313 | |
| 2314 | typedef struct { |
| 2315 | int n0, n1, nparams, ndelta; |
| 2316 | real kkk[2]; |
| 2317 | real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt; |
| 2318 | } t_luzar; |
| 2319 | |
| 2320 | static real compute_weighted_rates(int n, real t[], real ct[], real nt[], |
| 2321 | real kt[], real sigma_ct[], real sigma_nt[], |
| 2322 | real sigma_kt[], real *k, real *kp, |
| 2323 | real *sigma_k, real *sigma_kp, |
| 2324 | real fit_start) |
| 2325 | { |
| 2326 | #define NK10 10 |
| 2327 | int i, j; |
| 2328 | t_luzar tl; |
| 2329 | real kkk = 0, kkp = 0, kk2 = 0, kp2 = 0, chi2; |
| 2330 | |
| 2331 | *sigma_k = 0; |
| 2332 | *sigma_kp = 0; |
| 2333 | |
| 2334 | for (i = 0; (i < n); i++) |
| 2335 | { |
| 2336 | if (t[i] >= fit_start) |
| 2337 | { |
| 2338 | break; |
| 2339 | } |
| 2340 | } |
| 2341 | tl.n0 = i; |
| 2342 | tl.n1 = n; |
| 2343 | tl.nparams = 2; |
| 2344 | tl.ndelta = 1; |
| 2345 | tl.t = t; |
| 2346 | tl.ct = ct; |
| 2347 | tl.nt = nt; |
| 2348 | tl.kt = kt; |
| 2349 | tl.sigma_ct = sigma_ct; |
| 2350 | tl.sigma_nt = sigma_nt; |
| 2351 | tl.sigma_kt = sigma_kt; |
| 2352 | tl.kkk[0] = *k; |
| 2353 | tl.kkk[1] = *kp; |
| 2354 | |
| 2355 | chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */ |
| 2356 | *k = tl.kkk[0]; |
| 2357 | *kp = tl.kkk[1] = *kp; |
| 2358 | tl.ndelta = NK10; |
| 2359 | for (j = 0; (j < NK10); j++) |
| 2360 | { |
| 2361 | /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */ |
| 2362 | kkk += tl.kkk[0]; |
| 2363 | kkp += tl.kkk[1]; |
| 2364 | kk2 += sqr(tl.kkk[0]); |
| 2365 | kp2 += sqr(tl.kkk[1]); |
| 2366 | tl.n0++; |
| 2367 | } |
| 2368 | *sigma_k = sqrt(kk2/NK10 - sqr(kkk/NK10)); |
| 2369 | *sigma_kp = sqrt(kp2/NK10 - sqr(kkp/NK10)); |
| 2370 | |
| 2371 | return chi2; |
| 2372 | } |
| 2373 | |
| 2374 | static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start, |
| 2375 | const output_env_t oenv) |
| 2376 | { |
| 2377 | FILE *fp; |
| 2378 | real e_1, fitparm[4]; |
| 2379 | int i; |
| 2380 | |
| 2381 | e_1 = exp(-1); |
| 2382 | for (i = 0; (i < n); i++) |
| 2383 | { |
| 2384 | if (c[i] < e_1) |
| 2385 | { |
| 2386 | break; |
| 2387 | } |
| 2388 | } |
| 2389 | if (i < n) |
| 2390 | { |
| 2391 | fitparm[0] = t[i]; |
| 2392 | } |
| 2393 | else |
| 2394 | { |
| 2395 | fitparm[0] = 10; |
| 2396 | } |
| 2397 | fitparm[1] = 0.95; |
| 2398 | do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(), effnEXP2, fitparm, 0); |
| 2399 | } |
| 2400 | |
| 2401 | void analyse_corr(int n, real t[], real ct[], real nt[], real kt[], |
| 2402 | real sigma_ct[], real sigma_nt[], real sigma_kt[], |
| 2403 | real fit_start, real temp, real smooth_tail_start, |
| 2404 | const output_env_t oenv) |
| 2405 | { |
| 2406 | int i0, i; |
| 2407 | real k = 1, kp = 1, kow = 1; |
| 2408 | real Q = 0, chi22, chi2, dg, dgp, tau_hb, dtau, tau_rlx, e_1, dt, sigma_k, sigma_kp, ddg; |
| 2409 | double tmp, sn2 = 0, sc2 = 0, sk2 = 0, scn = 0, sck = 0, snk = 0; |
| 2410 | gmx_bool bError = (sigma_ct != NULL((void*)0)) && (sigma_nt != NULL((void*)0)) && (sigma_kt != NULL((void*)0)); |
| 2411 | |
| 2412 | if (smooth_tail_start >= 0) |
| 2413 | { |
| 2414 | smooth_tail(n, t, ct, sigma_ct, smooth_tail_start, oenv); |
| 2415 | smooth_tail(n, t, nt, sigma_nt, smooth_tail_start, oenv); |
| 2416 | smooth_tail(n, t, kt, sigma_kt, smooth_tail_start, oenv); |
| 2417 | } |
| 2418 | for (i0 = 0; (i0 < n-2) && ((t[i0]-t[0]) < fit_start); i0++) |
| 2419 | { |
| 2420 | ; |
| 2421 | } |
| 2422 | if (i0 < n-2) |
| 2423 | { |
| 2424 | for (i = i0; (i < n); i++) |
| 2425 | { |
| 2426 | sc2 += sqr(ct[i]); |
| 2427 | sn2 += sqr(nt[i]); |
| 2428 | sk2 += sqr(kt[i]); |
| 2429 | sck += ct[i]*kt[i]; |
| 2430 | snk += nt[i]*kt[i]; |
| 2431 | scn += ct[i]*nt[i]; |
| 2432 | } |
| 2433 | printf("Hydrogen bond thermodynamics at T = %g K\n", temp); |
| 2434 | tmp = (sn2*sc2-sqr(scn)); |
| 2435 | if ((tmp > 0) && (sn2 > 0)) |
| 2436 | { |
| 2437 | k = (sn2*sck-scn*snk)/tmp; |
| 2438 | kp = (k*scn-snk)/sn2; |
| 2439 | if (bError) |
| 2440 | { |
| 2441 | chi2 = 0; |
| 2442 | for (i = i0; (i < n); i++) |
| 2443 | { |
| 2444 | chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]); |
| 2445 | } |
| 2446 | chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt, |
| 2447 | sigma_kt, &k, &kp, |
| 2448 | &sigma_k, &sigma_kp, fit_start); |
| 2449 | Q = 0; /* quality_of_fit(chi2, 2);*/ |
| 2450 | ddg = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp*sigma_k/k; |
| 2451 | printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n", |
| 2452 | chi2, Q); |
| 2453 | printf("The Rate and Delta G are followed by an error estimate\n"); |
| 2454 | printf("----------------------------------------------------------\n" |
| 2455 | "Type Rate (1/ps) Sigma Time (ps) DG (kJ/mol) Sigma\n"); |
| 2456 | printf("Forward %10.3f %6.2f %8.3f %10.3f %6.2f\n", |
| 2457 | k, sigma_k, 1/k, calc_dg(1/k, temp), ddg); |
| 2458 | ddg = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*temp*sigma_kp/kp; |
| 2459 | printf("Backward %10.3f %6.2f %8.3f %10.3f %6.2f\n", |
| 2460 | kp, sigma_kp, 1/kp, calc_dg(1/kp, temp), ddg); |
| 2461 | } |
| 2462 | else |
| 2463 | { |
| 2464 | chi2 = 0; |
| 2465 | for (i = i0; (i < n); i++) |
| 2466 | { |
| 2467 | chi2 += sqr(k*ct[i]-kp*nt[i]-kt[i]); |
| 2468 | } |
| 2469 | printf("Fitting parameters chi^2 = %10g\nQ = %10g\n", |
| 2470 | chi2, Q); |
| 2471 | printf("--------------------------------------------------\n" |
| 2472 | "Type Rate (1/ps) Time (ps) DG (kJ/mol) Chi^2\n"); |
| 2473 | printf("Forward %10.3f %8.3f %10.3f %10g\n", |
| 2474 | k, 1/k, calc_dg(1/k, temp), chi2); |
| 2475 | printf("Backward %10.3f %8.3f %10.3f\n", |
| 2476 | kp, 1/kp, calc_dg(1/kp, temp)); |
| 2477 | } |
| 2478 | } |
| 2479 | if (sc2 > 0) |
| 2480 | { |
| 2481 | kow = 2*sck/sc2; |
| 2482 | printf("One-way %10.3f %s%8.3f %10.3f\n", |
| 2483 | kow, bError ? " " : "", 1/kow, calc_dg(1/kow, temp)); |
| 2484 | } |
| 2485 | else |
| 2486 | { |
| 2487 | printf(" - Numerical problems computing HB thermodynamics:\n" |
| 2488 | "sc2 = %g sn2 = %g sk2 = %g sck = %g snk = %g scn = %g\n", |
| 2489 | sc2, sn2, sk2, sck, snk, scn); |
| 2490 | } |
| 2491 | /* Determine integral of the correlation function */ |
| 2492 | tau_hb = evaluate_integral(n, t, ct, NULL((void*)0), (t[n-1]-t[0])/2, &dtau); |
| 2493 | printf("Integral %10.3f %s%8.3f %10.3f\n", 1/tau_hb, |
| 2494 | bError ? " " : "", tau_hb, calc_dg(tau_hb, temp)); |
| 2495 | e_1 = exp(-1); |
| 2496 | for (i = 0; (i < n-2); i++) |
| 2497 | { |
| 2498 | if ((ct[i] > e_1) && (ct[i+1] <= e_1)) |
| 2499 | { |
| 2500 | break; |
| 2501 | } |
| 2502 | } |
| 2503 | if (i < n-2) |
| 2504 | { |
| 2505 | /* Determine tau_relax from linear interpolation */ |
| 2506 | tau_rlx = t[i]-t[0] + (e_1-ct[i])*(t[i+1]-t[i])/(ct[i+1]-ct[i]); |
| 2507 | printf("Relaxation %10.3f %8.3f %s%10.3f\n", 1/tau_rlx, |
| 2508 | tau_rlx, bError ? " " : "", |
| 2509 | calc_dg(tau_rlx, temp)); |
| 2510 | } |
| 2511 | } |
| 2512 | else |
| 2513 | { |
| 2514 | printf("Correlation functions too short to compute thermodynamics\n"); |
| 2515 | } |
| 2516 | } |
| 2517 | |
| 2518 | void compute_derivative(int nn, real x[], real y[], real dydx[]) |
| 2519 | { |
| 2520 | int j; |
| 2521 | |
| 2522 | /* Compute k(t) = dc(t)/dt */ |
| 2523 | for (j = 1; (j < nn-1); j++) |
| 2524 | { |
| 2525 | dydx[j] = (y[j+1]-y[j-1])/(x[j+1]-x[j-1]); |
| 2526 | } |
| 2527 | /* Extrapolate endpoints */ |
| 2528 | dydx[0] = 2*dydx[1] - dydx[2]; |
| 2529 | dydx[nn-1] = 2*dydx[nn-2] - dydx[nn-3]; |
| 2530 | } |
| 2531 | |
| 2532 | static void parallel_print(int *data, int nThreads) |
| 2533 | { |
| 2534 | /* This prints the donors on which each tread is currently working. */ |
| 2535 | int i; |
| 2536 | |
| 2537 | fprintf(stderrstderr, "\r"); |
| 2538 | for (i = 0; i < nThreads; i++) |
| 2539 | { |
| 2540 | fprintf(stderrstderr, "%-7i", data[i]); |
| 2541 | } |
| 2542 | } |
| 2543 | |
| 2544 | static void normalizeACF(real *ct, real *gt, int nhb, int len) |
| 2545 | { |
| 2546 | real ct_fac, gt_fac; |
| 2547 | int i; |
| 2548 | |
| 2549 | /* Xu and Berne use the same normalization constant */ |
| 2550 | |
| 2551 | ct_fac = 1.0/ct[0]; |
| 2552 | gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb; |
| 2553 | |
| 2554 | printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac); |
| 2555 | for (i = 0; i < len; i++) |
| 2556 | { |
| 2557 | ct[i] *= ct_fac; |
| 2558 | if (gt != NULL((void*)0)) |
| 2559 | { |
| 2560 | gt[i] *= gt_fac; |
| 2561 | } |
| 2562 | } |
| 2563 | } |
| 2564 | |
| 2565 | /* Added argument bContact in do_hbac(...). Also |
| 2566 | * added support for -contact in the actual code. |
| 2567 | * - Erik Marklund, May 31, 2006 */ |
| 2568 | /* Changed contact code and added argument R2 |
| 2569 | * - Erik Marklund, June 29, 2006 |
| 2570 | */ |
| 2571 | static void do_hbac(const char *fn, t_hbdata *hb, |
| 2572 | int nDump, gmx_bool bMerge, gmx_bool bContact, real fit_start, |
| 2573 | real temp, gmx_bool R2, real smooth_tail_start, const output_env_t oenv, |
| 2574 | const char *gemType, int nThreads, |
| 2575 | const int NN, const gmx_bool bBallistic, const gmx_bool bGemFit) |
| 2576 | { |
| 2577 | FILE *fp; |
| 2578 | int i, j, k, m, n, o, nd, ihb, idist, n2, nn, iter, nSets; |
| 2579 | const char *legNN[] = { |
| 2580 | "Ac(t)", |
| 2581 | "Ac'(t)" |
| 2582 | }; |
| 2583 | static char **legGem; |
| 2584 | |
| 2585 | const char *legLuzar[] = { |
| 2586 | "Ac\\sfin sys\\v{}\\z{}(t)", |
| 2587 | "Ac(t)", |
| 2588 | "Cc\\scontact,hb\\v{}\\z{}(t)", |
| 2589 | "-dAc\\sfs\\v{}\\z{}/dt" |
| 2590 | }; |
| 2591 | gmx_bool bNorm = FALSE0, bOMP = FALSE0; |
| 2592 | double nhb = 0; |
| 2593 | int nhbi = 0; |
| 2594 | real *rhbex = NULL((void*)0), *ht, *gt, *ght, *dght, *kt; |
| 2595 | real *ct, *p_ct, tail, tail2, dtail, ct_fac, ght_fac, *cct; |
| 2596 | const real tol = 1e-3; |
| 2597 | int nframes = hb->nframes, nf; |
| 2598 | unsigned int **h = NULL((void*)0), **g = NULL((void*)0); |
| 2599 | int nh, nhbonds, nhydro, ngh; |
| 2600 | t_hbond *hbh; |
| 2601 | PSTYPEint p, *pfound = NULL((void*)0), np; |
| 2602 | t_pShift *pHist; |
| 2603 | int *ptimes = NULL((void*)0), *poff = NULL((void*)0), anhb, n0, mMax = INT_MIN(-2147483647 -1); |
| 2604 | real **rHbExGem = NULL((void*)0); |
| 2605 | gmx_bool c; |
| 2606 | int acType; |
| 2607 | t_E *E; |
| 2608 | double *ctdouble, *timedouble, *fittedct; |
| 2609 | double fittolerance = 0.1; |
| 2610 | int *dondata = NULL((void*)0), thisThread; |
| 2611 | |
| 2612 | enum { |
| 2613 | AC_NONE, AC_NN, AC_GEM, AC_LUZAR |
| 2614 | }; |
| 2615 | |
| 2616 | #ifdef GMX_OPENMP |
| 2617 | bOMP = TRUE1; |
| 2618 | #else |
| 2619 | bOMP = FALSE0; |
| 2620 | #endif |
| 2621 | |
| 2622 | printf("Doing autocorrelation "); |
| 2623 | |
| 2624 | /* Decide what kind of ACF calculations to do. */ |
| 2625 | if (NN > NN_NONE && NN < NN_NR) |
| 2626 | { |
| 2627 | #ifdef HAVE_NN_LOOPS |
| 2628 | acType = AC_NN; |
| 2629 | printf("using the energy estimate.\n"); |
| 2630 | #else |
| 2631 | acType = AC_NONE; |
| 2632 | printf("Can't do the NN-loop. Yet.\n"); |
| 2633 | #endif |
| 2634 | } |
| 2635 | else if (hb->bGem) |
| 2636 | { |
| 2637 | acType = AC_GEM; |
| 2638 | printf("according to the reversible geminate recombination model by Omer Markowitch.\n"); |
| 2639 | |
| 2640 | nSets = 1 + (bBallistic ? 1 : 0) + (bGemFit ? 1 : 0); |
| 2641 | snew(legGem, nSets)(legGem) = save_calloc("legGem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2641, (nSets), sizeof(*(legGem))); |
| 2642 | for (i = 0; i < nSets; i++) |
| 2643 | { |
| 2644 | snew(legGem[i], 128)(legGem[i]) = save_calloc("legGem[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2644, (128), sizeof(*(legGem[i]))); |
| 2645 | } |
| 2646 | sprintf(legGem[0], "Ac\\s%s\\v{}\\z{}(t)", gemType); |
| 2647 | if (bBallistic) |
| 2648 | { |
| 2649 | sprintf(legGem[1], "Ac'(t)"); |
| 2650 | } |
| 2651 | if (bGemFit) |
| 2652 | { |
| 2653 | sprintf(legGem[(bBallistic ? 3 : 2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType); |
| 2654 | } |
| 2655 | |
| 2656 | } |
| 2657 | else |
| 2658 | { |
| 2659 | acType = AC_LUZAR; |
| 2660 | printf("according to the theory of Luzar and Chandler.\n"); |
| 2661 | } |
| 2662 | fflush(stdoutstdout); |
| 2663 | |
| 2664 | /* build hbexist matrix in reals for autocorr */ |
| 2665 | /* Allocate memory for computing ACF (rhbex) and aggregating the ACF (ct) */ |
| 2666 | n2 = 1; |
| 2667 | while (n2 < nframes) |
| 2668 | { |
| 2669 | n2 *= 2; |
| 2670 | } |
| 2671 | |
| 2672 | nn = nframes/2; |
| 2673 | |
| 2674 | if (acType != AC_NN || bOMP) |
| 2675 | { |
| 2676 | snew(h, hb->maxhydro)(h) = save_calloc("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2676, (hb->maxhydro), sizeof(*(h))); |
| 2677 | snew(g, hb->maxhydro)(g) = save_calloc("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2677, (hb->maxhydro), sizeof(*(g))); |
| 2678 | } |
| 2679 | |
| 2680 | /* Dump hbonds for debugging */ |
| 2681 | dump_ac(hb, bMerge || bContact, nDump); |
| 2682 | |
| 2683 | /* Total number of hbonds analyzed here */ |
| 2684 | nhbonds = 0; |
| 2685 | ngh = 0; |
| 2686 | anhb = 0; |
Value stored to 'anhb' is never read | |
| 2687 | |
| 2688 | if (acType != AC_LUZAR && bOMP) |
| 2689 | { |
| 2690 | nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads())((((nThreads <= 0) ? 2147483647 : nThreads) < (gmx_omp_get_max_threads ())) ? ((nThreads <= 0) ? 2147483647 : nThreads) : (gmx_omp_get_max_threads ()) ); |
| 2691 | |
| 2692 | gmx_omp_set_num_threads(nThreads); |
| 2693 | snew(dondata, nThreads)(dondata) = save_calloc("dondata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2693, (nThreads), sizeof(*(dondata))); |
| 2694 | for (i = 0; i < nThreads; i++) |
| 2695 | { |
| 2696 | dondata[i] = -1; |
| 2697 | } |
| 2698 | printf("ACF calculations parallelized with OpenMP using %i threads.\n" |
| 2699 | "Expect close to linear scaling over this donor-loop.\n", nThreads); |
| 2700 | fflush(stdoutstdout); |
| 2701 | fprintf(stderrstderr, "Donors: [thread no]\n"); |
| 2702 | { |
| 2703 | char tmpstr[7]; |
| 2704 | for (i = 0; i < nThreads; i++) |
| 2705 | { |
| 2706 | snprintf(tmpstr, 7, "[%i]", i); |
| 2707 | fprintf(stderrstderr, "%-7s", tmpstr); |
| 2708 | } |
| 2709 | } |
| 2710 | fprintf(stderrstderr, "\n"); |
| 2711 | } |
| 2712 | |
| 2713 | |
| 2714 | /* Build the ACF according to acType */ |
| 2715 | switch (acType) |
| 2716 | { |
| 2717 | |
| 2718 | case AC_NN: |
| 2719 | #ifdef HAVE_NN_LOOPS |
| 2720 | /* Here we're using the estimated energy for the hydrogen bonds. */ |
| 2721 | snew(ct, nn)(ct) = save_calloc("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2721, (nn), sizeof(*(ct))); |
| 2722 | |
| 2723 | #pragma omp parallel \ |
| 2724 | private(i, j, k, nh, E, rhbex, thisThread) \ |
| 2725 | default(shared) |
| 2726 | { |
| 2727 | #pragma omp barrier |
| 2728 | thisThread = gmx_omp_get_thread_num(); |
| 2729 | rhbex = NULL((void*)0); |
| 2730 | |
| 2731 | snew(rhbex, n2)(rhbex) = save_calloc("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2731, (n2), sizeof(*(rhbex))); |
| 2732 | memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */ |
| 2733 | |
| 2734 | #pragma omp barrier |
| 2735 | #pragma omp for schedule (dynamic) |
| 2736 | for (i = 0; i < hb->d.nrd; i++) /* loop over donors */ |
| 2737 | { |
| 2738 | if (bOMP) |
| 2739 | { |
| 2740 | #pragma omp critical |
| 2741 | { |
| 2742 | dondata[thisThread] = i; |
| 2743 | parallel_print(dondata, nThreads); |
| 2744 | } |
| 2745 | } |
| 2746 | else |
| 2747 | { |
| 2748 | fprintf(stderrstderr, "\r %i", i); |
| 2749 | } |
| 2750 | |
| 2751 | for (j = 0; j < hb->a.nra; j++) /* loop over acceptors */ |
| 2752 | { |
| 2753 | for (nh = 0; nh < hb->d.nhydro[i]; nh++) /* loop over donors' hydrogens */ |
| 2754 | { |
| 2755 | E = hb->hbE.E[i][j][nh]; |
| 2756 | if (E != NULL((void*)0)) |
| 2757 | { |
| 2758 | for (k = 0; k < nframes; k++) |
| 2759 | { |
| 2760 | if (E[k] != NONSENSE_E) |
| 2761 | { |
| 2762 | rhbex[k] = (real)E[k]; |
| 2763 | } |
| 2764 | } |
| 2765 | |
| 2766 | low_do_autocorr(NULL((void*)0), oenv, NULL((void*)0), nframes, 1, -1, &(rhbex), hb->time[1]-hb->time[0], |
| 2767 | eacNormal(1<<0), 1, FALSE0, bNorm, FALSE0, 0, -1, 0, 1); |
| 2768 | #pragma omp critical |
| 2769 | { |
| 2770 | for (k = 0; (k < nn); k++) |
| 2771 | { |
| 2772 | ct[k] += rhbex[k]; |
| 2773 | } |
| 2774 | } |
| 2775 | } |
| 2776 | } /* k loop */ |
| 2777 | } /* j loop */ |
| 2778 | } /* i loop */ |
| 2779 | sfree(rhbex)save_free("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2779, (rhbex)); |
| 2780 | #pragma omp barrier |
| 2781 | } |
| 2782 | |
| 2783 | if (bOMP) |
| 2784 | { |
| 2785 | sfree(dondata)save_free("dondata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2785, (dondata)); |
| 2786 | } |
| 2787 | normalizeACF(ct, NULL((void*)0), 0, nn); |
| 2788 | snew(ctdouble, nn)(ctdouble) = save_calloc("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2788, (nn), sizeof(*(ctdouble))); |
| 2789 | snew(timedouble, nn)(timedouble) = save_calloc("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2789, (nn), sizeof(*(timedouble))); |
| 2790 | for (j = 0; j < nn; j++) |
| 2791 | { |
| 2792 | timedouble[j] = (double)(hb->time[j]); |
| 2793 | ctdouble[j] = (double)(ct[j]); |
| 2794 | } |
| 2795 | |
| 2796 | /* Remove ballistic term */ |
| 2797 | /* Ballistic component removal and fitting to the reversible geminate recombination model |
| 2798 | * will be taken out for the time being. First of all, one can remove the ballistic |
| 2799 | * component with g_analyze afterwards. Secondly, and more importantly, there are still |
| 2800 | * problems with the robustness of the fitting to the model. More work is needed. |
| 2801 | * A third reason is that we're currently using gsl for this and wish to reduce dependence |
| 2802 | * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with |
| 2803 | * a BSD-licence that can do the job. |
| 2804 | * |
| 2805 | * - Erik Marklund, June 18 2010. |
| 2806 | */ |
| 2807 | /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ |
| 2808 | /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ |
| 2809 | /* else */ |
| 2810 | /* printf("\nNumber of data points is less than the number of parameters to fit\n." */ |
| 2811 | /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ |
| 2812 | |
| 2813 | fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)"); |
| 2814 | xvgr_legend(fp, asize(legNN)((int)(sizeof(legNN)/sizeof((legNN)[0]))), legNN); |
| 2815 | |
| 2816 | for (j = 0; (j < nn); j++) |
| 2817 | { |
| 2818 | fprintf(fp, "%10g %10g %10g\n", |
| 2819 | hb->time[j]-hb->time[0], |
| 2820 | ct[j], |
| 2821 | ctdouble[j]); |
| 2822 | } |
| 2823 | xvgrclose(fp); |
| 2824 | sfree(ct)save_free("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2824, (ct)); |
| 2825 | sfree(ctdouble)save_free("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2825, (ctdouble)); |
| 2826 | sfree(timedouble)save_free("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2826, (timedouble)); |
| 2827 | #endif /* HAVE_NN_LOOPS */ |
| 2828 | break; /* case AC_NN */ |
| 2829 | |
| 2830 | case AC_GEM: |
| 2831 | snew(ct, 2*n2)(ct) = save_calloc("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2831, (2*n2), sizeof(*(ct))); |
| 2832 | memset(ct, 0, 2*n2*sizeof(real)); |
| 2833 | #ifndef GMX_OPENMP |
| 2834 | fprintf(stderrstderr, "Donor:\n"); |
| 2835 | #define __ACDATAct ct |
| 2836 | #else |
| 2837 | #define __ACDATAct p_ct |
| 2838 | #endif |
| 2839 | |
| 2840 | #pragma omp parallel \ |
| 2841 | private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \ |
| 2842 | pfound, poff, rHbExGem, p, ihb, mMax, \ |
| 2843 | thisThread, p_ct) \ |
| 2844 | default(shared) |
| 2845 | { /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */ |
| 2846 | h = NULL((void*)0); |
| 2847 | g = NULL((void*)0); |
| 2848 | thisThread = gmx_omp_get_thread_num(); |
| 2849 | snew(h, hb->maxhydro)(h) = save_calloc("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2849, (hb->maxhydro), sizeof(*(h))); |
| 2850 | snew(g, hb->maxhydro)(g) = save_calloc("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2850, (hb->maxhydro), sizeof(*(g))); |
| 2851 | mMax = INT_MIN(-2147483647 -1); |
| 2852 | rHbExGem = NULL((void*)0); |
| 2853 | poff = NULL((void*)0); |
| 2854 | pfound = NULL((void*)0); |
| 2855 | p_ct = NULL((void*)0); |
| 2856 | snew(p_ct, 2*n2)(p_ct) = save_calloc("p_ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2856, (2*n2), sizeof(*(p_ct))); |
| 2857 | memset(p_ct, 0, 2*n2*sizeof(real)); |
| 2858 | |
| 2859 | /* I'm using a chunk size of 1, since I expect \ |
| 2860 | * the overhead to be really small compared \ |
| 2861 | * to the actual calculations \ */ |
| 2862 | #pragma omp for schedule(dynamic,1) nowait |
| 2863 | for (i = 0; i < hb->d.nrd; i++) |
| 2864 | { |
| 2865 | |
| 2866 | if (bOMP) |
| 2867 | { |
| 2868 | #pragma omp critical |
| 2869 | { |
| 2870 | dondata[thisThread] = i; |
| 2871 | parallel_print(dondata, nThreads); |
| 2872 | } |
| 2873 | } |
| 2874 | else |
| 2875 | { |
| 2876 | fprintf(stderrstderr, "\r %i", i); |
| 2877 | } |
| 2878 | for (k = 0; k < hb->a.nra; k++) |
| 2879 | { |
| 2880 | for (nh = 0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) |
| 2881 | { |
| 2882 | hbh = hb->hbmap[i][k]; |
| 2883 | if (hbh) |
| 2884 | { |
| 2885 | /* Note that if hb->per->gemtype==gemDD, then distances will be stored in |
| 2886 | * hb->hbmap[d][a].h array anyway, because the contact flag will be set. |
| 2887 | * hence, it's only with the gemAD mode that hb->hbmap[d][a].g will be used. */ |
| 2888 | pHist = &(hb->per->pHist[i][k]); |
| 2889 | if (ISHB(hbh->history[nh])(((hbh->history[nh]) & 2) == 2) && pHist->len != 0) |
| 2890 | { |
| 2891 | |
| 2892 | { |
| 2893 | h[nh] = hbh->h[nh]; |
| 2894 | g[nh] = hb->per->gemtype == gemAD ? hbh->g[nh] : NULL((void*)0); |
| 2895 | } |
| 2896 | n0 = hbh->n0; |
| 2897 | nf = hbh->nframes; |
| 2898 | /* count the number of periodic shifts encountered and store |
| 2899 | * them in separate arrays. */ |
| 2900 | np = 0; |
| 2901 | for (j = 0; j < pHist->len; j++) |
| 2902 | { |
| 2903 | p = pHist->p[j]; |
| 2904 | for (m = 0; m <= np; m++) |
| 2905 | { |
| 2906 | if (m == np) /* p not recognized in list. Add it and set up new array. */ |
| 2907 | { |
| 2908 | np++; |
| 2909 | if (np > hb->per->nper) |
| 2910 | { |
| 2911 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2911, "Too many pshifts. Something's utterly wrong here."); |
| 2912 | } |
| 2913 | if (m >= mMax) /* Extend the arrays. |
| 2914 | * Doing it like this, using mMax to keep track of the sizes, |
| 2915 | * eleviates the need for freeing and re-allocating the arrays |
| 2916 | * when taking on the next donor-acceptor pair */ |
| 2917 | { |
| 2918 | mMax = m; |
| 2919 | srenew(pfound, np)(pfound) = save_realloc("pfound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2919, (pfound), (np), sizeof(*(pfound))); /* The list of found periodic shifts. */ |
| 2920 | srenew(rHbExGem, np)(rHbExGem) = save_realloc("rHbExGem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2920, (rHbExGem), (np), sizeof(*(rHbExGem))); /* The hb existence functions (-aver_hb). */ |
| 2921 | snew(rHbExGem[m], 2*n2)(rHbExGem[m]) = save_calloc("rHbExGem[m]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2921, (2*n2), sizeof(*(rHbExGem[m]))); |
| 2922 | srenew(poff, np)(poff) = save_realloc("poff", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2922, (poff), (np), sizeof(*(poff))); |
| 2923 | } |
| 2924 | |
| 2925 | { |
| 2926 | if (rHbExGem != NULL((void*)0) && rHbExGem[m] != NULL((void*)0)) |
| 2927 | { |
| 2928 | /* This must be done, as this array was most likey |
| 2929 | * used to store stuff in some previous iteration. */ |
| 2930 | memset(rHbExGem[m], 0, (sizeof(real)) * (2*n2)); |
| 2931 | } |
| 2932 | else |
| 2933 | { |
| 2934 | fprintf(stderrstderr, "rHbExGem not initialized! m = %i\n", m); |
| 2935 | } |
| 2936 | } |
| 2937 | pfound[m] = p; |
| 2938 | poff[m] = -1; |
| 2939 | |
| 2940 | break; |
| 2941 | } /* m==np */ |
| 2942 | if (p == pfound[m]) |
| 2943 | { |
| 2944 | break; |
| 2945 | } |
| 2946 | } /* m: Loop over found shifts */ |
| 2947 | } /* j: Loop over shifts */ |
| 2948 | |
| 2949 | /* Now unpack and disentangle the existence funtions. */ |
| 2950 | for (j = 0; j < nf; j++) |
| 2951 | { |
| 2952 | /* i: donor, |
| 2953 | * k: acceptor |
| 2954 | * nh: hydrogen |
| 2955 | * j: time |
| 2956 | * p: periodic shift |
| 2957 | * pfound: list of periodic shifts found for this pair. |
| 2958 | * poff: list of frame offsets; that is, the first |
| 2959 | * frame a hbond has a particular periodic shift. */ |
| 2960 | p = getPshift(*pHist, j+n0); |
| 2961 | if (p != -1) |
| 2962 | { |
| 2963 | for (m = 0; m < np; m++) |
| 2964 | { |
| 2965 | if (pfound[m] == p) |
| 2966 | { |
| 2967 | break; |
| 2968 | } |
| 2969 | if (m == (np-1)) |
| 2970 | { |
| 2971 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2971, "Shift not found, but must be there."); |
| 2972 | } |
| 2973 | } |
| 2974 | |
| 2975 | ihb = is_hb(h[nh], j) || ((hb->per->gemtype != gemAD || j == 0) ? FALSE0 : is_hb(g[nh], j)); |
| 2976 | if (ihb) |
| 2977 | { |
| 2978 | if (poff[m] == -1) |
| 2979 | { |
| 2980 | poff[m] = j; /* Here's where the first hbond with shift p is, |
| 2981 | * relative to the start of h[0].*/ |
| 2982 | } |
| 2983 | if (j < poff[m]) |
| 2984 | { |
| 2985 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 2985, "j<poff[m]"); |
| 2986 | } |
| 2987 | rHbExGem[m][j-poff[m]] += 1; |
| 2988 | } |
| 2989 | } |
| 2990 | } |
| 2991 | |
| 2992 | /* Now, build ac. */ |
| 2993 | for (m = 0; m < np; m++) |
| 2994 | { |
| 2995 | if (rHbExGem[m][0] > 0 && n0+poff[m] < nn /* && m==0 */) |
| 2996 | { |
| 2997 | low_do_autocorr(NULL((void*)0), oenv, NULL((void*)0), nframes, 1, -1, &(rHbExGem[m]), hb->time[1]-hb->time[0], |
| 2998 | eacNormal(1<<0), 1, FALSE0, bNorm, FALSE0, 0, -1, 0); |
| 2999 | for (j = 0; (j < nn); j++) |
| 3000 | { |
| 3001 | __ACDATAct[j] += rHbExGem[m][j]; |
| 3002 | } |
| 3003 | } |
| 3004 | } /* Building of ac. */ |
| 3005 | } /* if (ISHB(...*/ |
| 3006 | } /* if (hbh) */ |
| 3007 | } /* hydrogen loop */ |
| 3008 | } /* acceptor loop */ |
| 3009 | } /* donor loop */ |
| 3010 | |
| 3011 | for (m = 0; m <= mMax; m++) |
| 3012 | { |
| 3013 | sfree(rHbExGem[m])save_free("rHbExGem[m]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3013, (rHbExGem[m])); |
| 3014 | } |
| 3015 | sfree(pfound)save_free("pfound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3015, (pfound)); |
| 3016 | sfree(poff)save_free("poff", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3016, (poff)); |
| 3017 | sfree(rHbExGem)save_free("rHbExGem", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3017, (rHbExGem)); |
| 3018 | |
| 3019 | sfree(h)save_free("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3019, (h)); |
| 3020 | sfree(g)save_free("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3020, (g)); |
| 3021 | |
| 3022 | if (bOMP) |
| 3023 | { |
| 3024 | #pragma omp critical |
| 3025 | { |
| 3026 | for (i = 0; i < nn; i++) |
| 3027 | { |
| 3028 | ct[i] += p_ct[i]; |
| 3029 | } |
| 3030 | } |
| 3031 | sfree(p_ct)save_free("p_ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3031, (p_ct)); |
| 3032 | } |
| 3033 | |
| 3034 | } /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */ |
| 3035 | if (bOMP) |
| 3036 | { |
| 3037 | sfree(dondata)save_free("dondata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3037, (dondata)); |
| 3038 | } |
| 3039 | |
| 3040 | normalizeACF(ct, NULL((void*)0), 0, nn); |
| 3041 | |
| 3042 | fprintf(stderrstderr, "\n\nACF successfully calculated.\n"); |
| 3043 | |
| 3044 | /* Use this part to fit to geminate recombination - JCP 129, 84505 (2008) */ |
| 3045 | |
| 3046 | snew(ctdouble, nn)(ctdouble) = save_calloc("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3046, (nn), sizeof(*(ctdouble))); |
| 3047 | snew(timedouble, nn)(timedouble) = save_calloc("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3047, (nn), sizeof(*(timedouble))); |
| 3048 | snew(fittedct, nn)(fittedct) = save_calloc("fittedct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3048, (nn), sizeof(*(fittedct))); |
| 3049 | |
| 3050 | for (j = 0; j < nn; j++) |
| 3051 | { |
| 3052 | timedouble[j] = (double)(hb->time[j]); |
| 3053 | ctdouble[j] = (double)(ct[j]); |
| 3054 | } |
| 3055 | |
| 3056 | /* Remove ballistic term */ |
| 3057 | /* Ballistic component removal and fitting to the reversible geminate recombination model |
| 3058 | * will be taken out for the time being. First of all, one can remove the ballistic |
| 3059 | * component with g_analyze afterwards. Secondly, and more importantly, there are still |
| 3060 | * problems with the robustness of the fitting to the model. More work is needed. |
| 3061 | * A third reason is that we're currently using gsl for this and wish to reduce dependence |
| 3062 | * on external libraries. There are Levenberg-Marquardt and nsimplex solvers that come with |
| 3063 | * a BSD-licence that can do the job. |
| 3064 | * |
| 3065 | * - Erik Marklund, June 18 2010. |
| 3066 | */ |
| 3067 | /* if (bBallistic) { */ |
| 3068 | /* if (params->ballistic/params->tDelta >= params->nExpFit*2+1) */ |
| 3069 | /* takeAwayBallistic(ctdouble, timedouble, nn, params->ballistic, params->nExpFit, params->bDt); */ |
| 3070 | /* else */ |
| 3071 | /* printf("\nNumber of data points is less than the number of parameters to fit\n." */ |
| 3072 | /* "The system is underdetermined, hence no ballistic term can be found.\n\n"); */ |
| 3073 | /* } */ |
| 3074 | /* if (bGemFit) */ |
| 3075 | /* fitGemRecomb(ctdouble, timedouble, &fittedct, nn, params); */ |
| 3076 | |
| 3077 | |
| 3078 | if (bContact) |
| 3079 | { |
| 3080 | fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
| 3081 | } |
| 3082 | else |
| 3083 | { |
| 3084 | fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
| 3085 | } |
| 3086 | xvgr_legend(fp, asize(legGem)((int)(sizeof(legGem)/sizeof((legGem)[0]))), (const char**)legGem, oenv); |
| 3087 | |
| 3088 | for (j = 0; (j < nn); j++) |
| 3089 | { |
| 3090 | fprintf(fp, "%10g %10g", hb->time[j]-hb->time[0], ct[j]); |
| 3091 | if (bBallistic) |
| 3092 | { |
| 3093 | fprintf(fp, " %10g", ctdouble[j]); |
| 3094 | } |
| 3095 | if (bGemFit) |
| 3096 | { |
| 3097 | fprintf(fp, " %10g", fittedct[j]); |
| 3098 | } |
| 3099 | fprintf(fp, "\n"); |
| 3100 | } |
| 3101 | xvgrclose(fp); |
| 3102 | |
| 3103 | sfree(ctdouble)save_free("ctdouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3103, (ctdouble)); |
| 3104 | sfree(timedouble)save_free("timedouble", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3104, (timedouble)); |
| 3105 | sfree(fittedct)save_free("fittedct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3105, (fittedct)); |
| 3106 | sfree(ct)save_free("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3106, (ct)); |
| 3107 | |
| 3108 | break; /* case AC_GEM */ |
| 3109 | |
| 3110 | case AC_LUZAR: |
| 3111 | snew(rhbex, 2*n2)(rhbex) = save_calloc("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3111, (2*n2), sizeof(*(rhbex))); |
| 3112 | snew(ct, 2*n2)(ct) = save_calloc("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3112, (2*n2), sizeof(*(ct))); |
| 3113 | snew(gt, 2*n2)(gt) = save_calloc("gt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3113, (2*n2), sizeof(*(gt))); |
| 3114 | snew(ht, 2*n2)(ht) = save_calloc("ht", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3114, (2*n2), sizeof(*(ht))); |
| 3115 | snew(ght, 2*n2)(ght) = save_calloc("ght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3115, (2*n2), sizeof(*(ght))); |
| 3116 | snew(dght, 2*n2)(dght) = save_calloc("dght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3116, (2*n2), sizeof(*(dght))); |
| 3117 | |
| 3118 | snew(kt, nn)(kt) = save_calloc("kt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3118, (nn), sizeof(*(kt))); |
| 3119 | snew(cct, nn)(cct) = save_calloc("cct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3119, (nn), sizeof(*(cct))); |
| 3120 | |
| 3121 | for (i = 0; (i < hb->d.nrd); i++) |
| 3122 | { |
| 3123 | for (k = 0; (k < hb->a.nra); k++) |
| 3124 | { |
| 3125 | nhydro = 0; |
| 3126 | hbh = hb->hbmap[i][k]; |
| 3127 | |
| 3128 | if (hbh) |
| 3129 | { |
| 3130 | if (bMerge || bContact) |
| 3131 | { |
| 3132 | if (ISHB(hbh->history[0])(((hbh->history[0]) & 2) == 2)) |
| 3133 | { |
| 3134 | h[0] = hbh->h[0]; |
| 3135 | g[0] = hbh->g[0]; |
| 3136 | nhydro = 1; |
| 3137 | } |
| 3138 | } |
| 3139 | else |
| 3140 | { |
| 3141 | for (m = 0; (m < hb->maxhydro); m++) |
| 3142 | { |
| 3143 | if (bContact ? ISDIST(hbh->history[m])(((hbh->history[m]) & 1) == 1) : ISHB(hbh->history[m])(((hbh->history[m]) & 2) == 2)) |
| 3144 | { |
| 3145 | g[nhydro] = hbh->g[m]; |
| 3146 | h[nhydro] = hbh->h[m]; |
| 3147 | nhydro++; |
| 3148 | } |
| 3149 | } |
| 3150 | } |
| 3151 | |
| 3152 | nf = hbh->nframes; |
| 3153 | for (nh = 0; (nh < nhydro); nh++) |
| 3154 | { |
| 3155 | int nrint = bContact ? hb->nrdist : hb->nrhb; |
| 3156 | if ((((nhbonds+1) % 10) == 0) || (nhbonds+1 == nrint)) |
| 3157 | { |
| 3158 | fprintf(stderrstderr, "\rACF %d/%d", nhbonds+1, nrint); |
| 3159 | } |
| 3160 | nhbonds++; |
| 3161 | for (j = 0; (j < nframes); j++) |
| 3162 | { |
| 3163 | /* Changed '<' into '<=' below, just like I did in |
| 3164 | the hbm-output-loop in the gmx_hbond() block. |
| 3165 | - Erik Marklund, May 31, 2006 */ |
| 3166 | if (j <= nf) |
| 3167 | { |
| 3168 | ihb = is_hb(h[nh], j); |
| 3169 | idist = is_hb(g[nh], j); |
| 3170 | } |
| 3171 | else |
| 3172 | { |
| 3173 | ihb = idist = 0; |
| 3174 | } |
| 3175 | rhbex[j] = ihb; |
| 3176 | /* For contacts: if a second cut-off is provided, use it, |
| 3177 | * otherwise use g(t) = 1-h(t) */ |
| 3178 | if (!R2 && bContact) |
| 3179 | { |
| 3180 | gt[j] = 1-ihb; |
| 3181 | } |
| 3182 | else |
| 3183 | { |
| 3184 | gt[j] = idist*(1-ihb); |
| 3185 | } |
| 3186 | ht[j] = rhbex[j]; |
| 3187 | nhb += ihb; |
| 3188 | } |
| 3189 | |
| 3190 | |
| 3191 | /* The autocorrelation function is normalized after summation only */ |
| 3192 | low_do_autocorr(NULL((void*)0), oenv, NULL((void*)0), nframes, 1, -1, &rhbex, hb->time[1]-hb->time[0], |
| 3193 | eacNormal(1<<0), 1, FALSE0, bNorm, FALSE0, 0, -1, 0); |
| 3194 | |
| 3195 | /* Cross correlation analysis for thermodynamics */ |
| 3196 | for (j = nframes; (j < n2); j++) |
| 3197 | { |
| 3198 | ht[j] = 0; |
| 3199 | gt[j] = 0; |
| 3200 | } |
| 3201 | |
| 3202 | cross_corr(n2, ht, gt, dght); |
| 3203 | |
| 3204 | for (j = 0; (j < nn); j++) |
| 3205 | { |
| 3206 | ct[j] += rhbex[j]; |
| 3207 | ght[j] += dght[j]; |
| 3208 | } |
| 3209 | } |
| 3210 | } |
| 3211 | } |
| 3212 | } |
| 3213 | fprintf(stderrstderr, "\n"); |
| 3214 | sfree(h)save_free("h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3214, (h)); |
| 3215 | sfree(g)save_free("g", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3215, (g)); |
| 3216 | normalizeACF(ct, ght, nhb, nn); |
| 3217 | |
| 3218 | /* Determine tail value for statistics */ |
| 3219 | tail = 0; |
| 3220 | tail2 = 0; |
| 3221 | for (j = nn/2; (j < nn); j++) |
| 3222 | { |
| 3223 | tail += ct[j]; |
| 3224 | tail2 += ct[j]*ct[j]; |
| 3225 | } |
| 3226 | tail /= (nn - nn/2); |
| 3227 | tail2 /= (nn - nn/2); |
| 3228 | dtail = sqrt(tail2-tail*tail); |
| 3229 | |
| 3230 | /* Check whether the ACF is long enough */ |
| 3231 | if (dtail > tol) |
| 3232 | { |
| 3233 | printf("\nWARNING: Correlation function is probably not long enough\n" |
| 3234 | "because the standard deviation in the tail of C(t) > %g\n" |
| 3235 | "Tail value (average C(t) over second half of acf): %g +/- %g\n", |
| 3236 | tol, tail, dtail); |
| 3237 | } |
| 3238 | for (j = 0; (j < nn); j++) |
| 3239 | { |
| 3240 | cct[j] = ct[j]; |
| 3241 | ct[j] = (cct[j]-tail)/(1-tail); |
| 3242 | } |
| 3243 | /* Compute negative derivative k(t) = -dc(t)/dt */ |
| 3244 | compute_derivative(nn, hb->time, ct, kt); |
| 3245 | for (j = 0; (j < nn); j++) |
| 3246 | { |
| 3247 | kt[j] = -kt[j]; |
| 3248 | } |
| 3249 | |
| 3250 | |
| 3251 | if (bContact) |
| 3252 | { |
| 3253 | fp = xvgropen(fn, "Contact Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
| 3254 | } |
| 3255 | else |
| 3256 | { |
| 3257 | fp = xvgropen(fn, "Hydrogen Bond Autocorrelation", output_env_get_xvgr_tlabel(oenv), "C(t)", oenv); |
| 3258 | } |
| 3259 | xvgr_legend(fp, asize(legLuzar)((int)(sizeof(legLuzar)/sizeof((legLuzar)[0]))), legLuzar, oenv); |
| 3260 | |
| 3261 | |
| 3262 | for (j = 0; (j < nn); j++) |
| 3263 | { |
| 3264 | fprintf(fp, "%10g %10g %10g %10g %10g\n", |
| 3265 | hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]); |
| 3266 | } |
| 3267 | gmx_ffclose(fp); |
| 3268 | |
| 3269 | analyse_corr(nn, hb->time, ct, ght, kt, NULL((void*)0), NULL((void*)0), NULL((void*)0), |
| 3270 | fit_start, temp, smooth_tail_start, oenv); |
| 3271 | |
| 3272 | do_view(oenv, fn, NULL((void*)0)); |
| 3273 | sfree(rhbex)save_free("rhbex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3273, (rhbex)); |
| 3274 | sfree(ct)save_free("ct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3274, (ct)); |
| 3275 | sfree(gt)save_free("gt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3275, (gt)); |
| 3276 | sfree(ht)save_free("ht", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3276, (ht)); |
| 3277 | sfree(ght)save_free("ght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3277, (ght)); |
| 3278 | sfree(dght)save_free("dght", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3278, (dght)); |
| 3279 | sfree(cct)save_free("cct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3279, (cct)); |
| 3280 | sfree(kt)save_free("kt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3280, (kt)); |
| 3281 | /* sfree(h); */ |
| 3282 | /* sfree(g); */ |
| 3283 | |
| 3284 | break; /* case AC_LUZAR */ |
| 3285 | |
| 3286 | default: |
| 3287 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3287, "Unrecognized type of ACF-calulation. acType = %i.", acType); |
| 3288 | } /* switch (acType) */ |
| 3289 | } |
| 3290 | |
| 3291 | static void init_hbframe(t_hbdata *hb, int nframes, real t) |
| 3292 | { |
| 3293 | int i, j, m; |
| 3294 | |
| 3295 | hb->time[nframes] = t; |
| 3296 | hb->nhb[nframes] = 0; |
| 3297 | hb->ndist[nframes] = 0; |
| 3298 | for (i = 0; (i < max_hx7); i++) |
| 3299 | { |
| 3300 | hb->nhx[nframes][i] = 0; |
| 3301 | } |
| 3302 | /* Loop invalidated */ |
| 3303 | if (hb->bHBmap && 0) |
| 3304 | { |
| 3305 | for (i = 0; (i < hb->d.nrd); i++) |
| 3306 | { |
| 3307 | for (j = 0; (j < hb->a.nra); j++) |
| 3308 | { |
| 3309 | for (m = 0; (m < hb->maxhydro); m++) |
| 3310 | { |
| 3311 | if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[m]) |
| 3312 | { |
| 3313 | set_hb(hb, i, m, j, nframes, HB_NO0); |
| 3314 | } |
| 3315 | } |
| 3316 | } |
| 3317 | } |
| 3318 | } |
| 3319 | /*set_hb(hb->hbmap[i][j]->h[m],nframes-hb->hbmap[i][j]->n0,HB_NO);*/ |
| 3320 | } |
| 3321 | |
| 3322 | static void analyse_donor_props(const char *fn, t_hbdata *hb, int nframes, real t, |
| 3323 | const output_env_t oenv) |
| 3324 | { |
| 3325 | static FILE *fp = NULL((void*)0); |
| 3326 | const char *leg[] = { "Nbound", "Nfree" }; |
| 3327 | int i, j, k, nbound, nb, nhtot; |
| 3328 | |
| 3329 | if (!fn) |
| 3330 | { |
| 3331 | return; |
| 3332 | } |
| 3333 | if (!fp) |
| 3334 | { |
| 3335 | fp = xvgropen(fn, "Donor properties", output_env_get_xvgr_tlabel(oenv), "Number", oenv); |
| 3336 | xvgr_legend(fp, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); |
| 3337 | } |
| 3338 | nbound = 0; |
| 3339 | nhtot = 0; |
| 3340 | for (i = 0; (i < hb->d.nrd); i++) |
| 3341 | { |
| 3342 | for (k = 0; (k < hb->d.nhydro[i]); k++) |
| 3343 | { |
| 3344 | nb = 0; |
| 3345 | nhtot++; |
| 3346 | for (j = 0; (j < hb->a.nra) && (nb == 0); j++) |
| 3347 | { |
| 3348 | if (hb->hbmap[i][j] && hb->hbmap[i][j]->h[k] && |
| 3349 | is_hb(hb->hbmap[i][j]->h[k], nframes)) |
| 3350 | { |
| 3351 | nb = 1; |
| 3352 | } |
| 3353 | } |
| 3354 | nbound += nb; |
| 3355 | } |
| 3356 | } |
| 3357 | fprintf(fp, "%10.3e %6d %6d\n", t, nbound, nhtot-nbound); |
| 3358 | } |
| 3359 | |
| 3360 | static void dump_hbmap(t_hbdata *hb, |
| 3361 | int nfile, t_filenm fnm[], gmx_bool bTwo, |
| 3362 | gmx_bool bContact, int isize[], int *index[], char *grpnames[], |
| 3363 | t_atoms *atoms) |
| 3364 | { |
| 3365 | FILE *fp, *fplog; |
| 3366 | int ddd, hhh, aaa, i, j, k, m, grp; |
| 3367 | char ds[32], hs[32], as[32]; |
| 3368 | gmx_bool first; |
| 3369 | |
| 3370 | fp = opt2FILE("-hbn", nfile, fnm, "w")gmx_ffopen(opt2fn("-hbn", nfile, fnm), "w"); |
| 3371 | if (opt2bSet("-g", nfile, fnm)) |
| 3372 | { |
| 3373 | fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w"); |
| 3374 | fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor"); |
| 3375 | } |
| 3376 | else |
| 3377 | { |
| 3378 | fplog = NULL((void*)0); |
| 3379 | } |
| 3380 | for (grp = gr0; grp <= (bTwo ? gr1 : gr0); grp++) |
| 3381 | { |
| 3382 | fprintf(fp, "[ %s ]", grpnames[grp]); |
| 3383 | for (i = 0; i < isize[grp]; i++) |
| 3384 | { |
| 3385 | fprintf(fp, (i%15) ? " " : "\n"); |
| 3386 | fprintf(fp, " %4u", index[grp][i]+1); |
| 3387 | } |
| 3388 | fprintf(fp, "\n"); |
| 3389 | /* |
| 3390 | Added -contact support below. |
| 3391 | - Erik Marklund, May 29, 2006 |
| 3392 | */ |
| 3393 | if (!bContact) |
| 3394 | { |
| 3395 | fprintf(fp, "[ donors_hydrogens_%s ]\n", grpnames[grp]); |
| 3396 | for (i = 0; (i < hb->d.nrd); i++) |
| 3397 | { |
| 3398 | if (hb->d.grp[i] == grp) |
| 3399 | { |
| 3400 | for (j = 0; (j < hb->d.nhydro[i]); j++) |
| 3401 | { |
| 3402 | fprintf(fp, " %4u %4u", hb->d.don[i]+1, |
| 3403 | hb->d.hydro[i][j]+1); |
| 3404 | } |
| 3405 | fprintf(fp, "\n"); |
| 3406 | } |
| 3407 | } |
| 3408 | first = TRUE1; |
| 3409 | fprintf(fp, "[ acceptors_%s ]", grpnames[grp]); |
| 3410 | for (i = 0; (i < hb->a.nra); i++) |
| 3411 | { |
| 3412 | if (hb->a.grp[i] == grp) |
| 3413 | { |
| 3414 | fprintf(fp, (i%15 && !first) ? " " : "\n"); |
| 3415 | fprintf(fp, " %4u", hb->a.acc[i]+1); |
| 3416 | first = FALSE0; |
| 3417 | } |
| 3418 | } |
| 3419 | fprintf(fp, "\n"); |
| 3420 | } |
| 3421 | } |
| 3422 | if (bTwo) |
| 3423 | { |
| 3424 | fprintf(fp, bContact ? "[ contacts_%s-%s ]\n" : |
| 3425 | "[ hbonds_%s-%s ]\n", grpnames[0], grpnames[1]); |
| 3426 | } |
| 3427 | else |
| 3428 | { |
| 3429 | fprintf(fp, bContact ? "[ contacts_%s ]" : "[ hbonds_%s ]\n", grpnames[0]); |
| 3430 | } |
| 3431 | |
| 3432 | for (i = 0; (i < hb->d.nrd); i++) |
| 3433 | { |
| 3434 | ddd = hb->d.don[i]; |
| 3435 | for (k = 0; (k < hb->a.nra); k++) |
| 3436 | { |
| 3437 | aaa = hb->a.acc[k]; |
| 3438 | for (m = 0; (m < hb->d.nhydro[i]); m++) |
| 3439 | { |
| 3440 | if (hb->hbmap[i][k] && ISHB(hb->hbmap[i][k]->history[m])(((hb->hbmap[i][k]->history[m]) & 2) == 2)) |
| 3441 | { |
| 3442 | sprintf(ds, "%s", mkatomname(atoms, ddd)); |
| 3443 | sprintf(as, "%s", mkatomname(atoms, aaa)); |
| 3444 | if (bContact) |
| 3445 | { |
| 3446 | fprintf(fp, " %6u %6u\n", ddd+1, aaa+1); |
| 3447 | if (fplog) |
| 3448 | { |
| 3449 | fprintf(fplog, "%12s %12s\n", ds, as); |
| 3450 | } |
| 3451 | } |
| 3452 | else |
| 3453 | { |
| 3454 | hhh = hb->d.hydro[i][m]; |
| 3455 | sprintf(hs, "%s", mkatomname(atoms, hhh)); |
| 3456 | fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1); |
| 3457 | if (fplog) |
| 3458 | { |
| 3459 | fprintf(fplog, "%12s %12s %12s\n", ds, hs, as); |
| 3460 | } |
| 3461 | } |
| 3462 | } |
| 3463 | } |
| 3464 | } |
| 3465 | } |
| 3466 | gmx_ffclose(fp); |
| 3467 | if (fplog) |
| 3468 | { |
| 3469 | gmx_ffclose(fplog); |
| 3470 | } |
| 3471 | } |
| 3472 | |
| 3473 | /* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template. |
| 3474 | * It mimics add_frames() and init_frame() to some extent. */ |
| 3475 | static void sync_hbdata(t_hbdata *p_hb, int nframes) |
| 3476 | { |
| 3477 | int i; |
| 3478 | if (nframes >= p_hb->max_frames) |
| 3479 | { |
| 3480 | p_hb->max_frames += 4096; |
| 3481 | srenew(p_hb->nhb, p_hb->max_frames)(p_hb->nhb) = save_realloc("p_hb->nhb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3481, (p_hb->nhb), (p_hb->max_frames), sizeof(*(p_hb-> nhb))); |
| 3482 | srenew(p_hb->ndist, p_hb->max_frames)(p_hb->ndist) = save_realloc("p_hb->ndist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3482, (p_hb->ndist), (p_hb->max_frames), sizeof(*(p_hb ->ndist))); |
| 3483 | srenew(p_hb->n_bound, p_hb->max_frames)(p_hb->n_bound) = save_realloc("p_hb->n_bound", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3483, (p_hb->n_bound), (p_hb->max_frames), sizeof(*(p_hb ->n_bound))); |
| 3484 | srenew(p_hb->nhx, p_hb->max_frames)(p_hb->nhx) = save_realloc("p_hb->nhx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3484, (p_hb->nhx), (p_hb->max_frames), sizeof(*(p_hb-> nhx))); |
| 3485 | if (p_hb->bDAnr) |
| 3486 | { |
| 3487 | srenew(p_hb->danr, p_hb->max_frames)(p_hb->danr) = save_realloc("p_hb->danr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3487, (p_hb->danr), (p_hb->max_frames), sizeof(*(p_hb ->danr))); |
| 3488 | } |
| 3489 | memset(&(p_hb->nhb[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); |
| 3490 | memset(&(p_hb->ndist[nframes]), 0, sizeof(int) * (p_hb->max_frames-nframes)); |
| 3491 | p_hb->nhb[nframes] = 0; |
| 3492 | p_hb->ndist[nframes] = 0; |
| 3493 | |
| 3494 | } |
| 3495 | p_hb->nframes = nframes; |
| 3496 | /* for (i=0;) */ |
| 3497 | /* { */ |
| 3498 | /* p_hb->nhx[nframes][i] */ |
| 3499 | /* } */ |
| 3500 | memset(&(p_hb->nhx[nframes]), 0, sizeof(int)*max_hx7); /* zero the helix count for this frame */ |
| 3501 | |
| 3502 | /* hb->per will remain constant througout the frame loop, |
| 3503 | * even though the data its members point to will change, |
| 3504 | * hence no need for re-syncing. */ |
| 3505 | } |
| 3506 | |
| 3507 | int gmx_hbond(int argc, char *argv[]) |
| 3508 | { |
| 3509 | const char *desc[] = { |
| 3510 | "[THISMODULE] computes and analyzes hydrogen bonds. Hydrogen bonds are", |
| 3511 | "determined based on cutoffs for the angle Hydrogen - Donor - Acceptor", |
| 3512 | "(zero is extended) and the distance Donor - Acceptor", |
| 3513 | "(or Hydrogen - Acceptor using [TT]-noda[tt]).", |
| 3514 | "OH and NH groups are regarded as donors, O is an acceptor always,", |
| 3515 | "N is an acceptor by default, but this can be switched using", |
| 3516 | "[TT]-nitacc[tt]. Dummy hydrogen atoms are assumed to be connected", |
| 3517 | "to the first preceding non-hydrogen atom.[PAR]", |
| 3518 | |
| 3519 | "You need to specify two groups for analysis, which must be either", |
| 3520 | "identical or non-overlapping. All hydrogen bonds between the two", |
| 3521 | "groups are analyzed.[PAR]", |
| 3522 | |
| 3523 | "If you set [TT]-shell[tt], you will be asked for an additional index group", |
| 3524 | "which should contain exactly one atom. In this case, only hydrogen", |
| 3525 | "bonds between atoms within the shell distance from the one atom are", |
| 3526 | "considered.[PAR]", |
| 3527 | |
| 3528 | "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler", |
| 3529 | "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).", |
| 3530 | "If contact kinetics are analyzed by using the -contact option, then", |
| 3531 | "n(t) can be defined as either all pairs that are not within contact distance r at time t", |
| 3532 | "(corresponding to leaving the -r2 option at the default value 0) or all pairs that", |
| 3533 | "are within distance r2 (corresponding to setting a second cut-off value with option -r2).", |
| 3534 | "See mentioned literature for more details and definitions." |
| 3535 | "[PAR]", |
| 3536 | |
| 3537 | /* "It is also possible to analyse specific hydrogen bonds with", |
| 3538 | "[TT]-sel[tt]. This index file must contain a group of atom triplets", |
| 3539 | "Donor Hydrogen Acceptor, in the following way:[PAR]", |
| 3540 | */ |
| 3541 | "[TT]", |
| 3542 | "[ selected ][BR]", |
| 3543 | " 20 21 24[BR]", |
| 3544 | " 25 26 29[BR]", |
| 3545 | " 1 3 6[BR]", |
| 3546 | "[tt][BR]", |
| 3547 | "Note that the triplets need not be on separate lines.", |
| 3548 | "Each atom triplet specifies a hydrogen bond to be analyzed,", |
| 3549 | "note also that no check is made for the types of atoms.[PAR]", |
| 3550 | |
| 3551 | "[BB]Output:[bb][BR]", |
| 3552 | "[TT]-num[tt]: number of hydrogen bonds as a function of time.[BR]", |
| 3553 | "[TT]-ac[tt]: average over all autocorrelations of the existence", |
| 3554 | "functions (either 0 or 1) of all hydrogen bonds.[BR]", |
| 3555 | "[TT]-dist[tt]: distance distribution of all hydrogen bonds.[BR]", |
| 3556 | "[TT]-ang[tt]: angle distribution of all hydrogen bonds.[BR]", |
| 3557 | "[TT]-hx[tt]: the number of n-n+i hydrogen bonds as a function of time", |
| 3558 | "where n and n+i stand for residue numbers and i ranges from 0 to 6.", |
| 3559 | "This includes the n-n+3, n-n+4 and n-n+5 hydrogen bonds associated", |
| 3560 | "with helices in proteins.[BR]", |
| 3561 | "[TT]-hbn[tt]: all selected groups, donors, hydrogens and acceptors", |
| 3562 | "for selected groups, all hydrogen bonded atoms from all groups and", |
| 3563 | "all solvent atoms involved in insertion.[BR]", |
| 3564 | "[TT]-hbm[tt]: existence matrix for all hydrogen bonds over all", |
| 3565 | "frames, this also contains information on solvent insertion", |
| 3566 | "into hydrogen bonds. Ordering is identical to that in [TT]-hbn[tt]", |
| 3567 | "index file.[BR]", |
| 3568 | "[TT]-dan[tt]: write out the number of donors and acceptors analyzed for", |
| 3569 | "each timeframe. This is especially useful when using [TT]-shell[tt].[BR]", |
| 3570 | "[TT]-nhbdist[tt]: compute the number of HBonds per hydrogen in order to", |
| 3571 | "compare results to Raman Spectroscopy.", |
| 3572 | "[PAR]", |
| 3573 | "Note: options [TT]-ac[tt], [TT]-life[tt], [TT]-hbn[tt] and [TT]-hbm[tt]", |
| 3574 | "require an amount of memory proportional to the total numbers of donors", |
| 3575 | "times the total number of acceptors in the selected group(s)." |
| 3576 | }; |
| 3577 | |
| 3578 | static real acut = 30, abin = 1, rcut = 0.35, r2cut = 0, rbin = 0.005, rshell = -1; |
| 3579 | static real maxnhb = 0, fit_start = 1, fit_end = 60, temp = 298.15, smooth_tail_start = -1, D = -1; |
| 3580 | static gmx_bool bNitAcc = TRUE1, bDA = TRUE1, bMerge = TRUE1; |
| 3581 | static int nDump = 0, nFitPoints = 100; |
| 3582 | static int nThreads = 0, nBalExp = 4; |
| 3583 | |
| 3584 | static gmx_bool bContact = FALSE0, bBallistic = FALSE0, bGemFit = FALSE0; |
| 3585 | static real logAfterTime = 10, gemBallistic = 0.2; /* ps */ |
| 3586 | static const char *NNtype[] = {NULL((void*)0), "none", "binary", "oneOverR3", "dipole", NULL((void*)0)}; |
| 3587 | |
| 3588 | /* options */ |
| 3589 | t_pargs pa [] = { |
| 3590 | { "-a", FALSE0, etREAL, {&acut}, |
| 3591 | "Cutoff angle (degrees, Hydrogen - Donor - Acceptor)" }, |
| 3592 | { "-r", FALSE0, etREAL, {&rcut}, |
| 3593 | "Cutoff radius (nm, X - Acceptor, see next option)" }, |
| 3594 | { "-da", FALSE0, etBOOL, {&bDA}, |
| 3595 | "Use distance Donor-Acceptor (if TRUE) or Hydrogen-Acceptor (FALSE)" }, |
| 3596 | { "-r2", FALSE0, etREAL, {&r2cut}, |
| 3597 | "Second cutoff radius. Mainly useful with [TT]-contact[tt] and [TT]-ac[tt]"}, |
| 3598 | { "-abin", FALSE0, etREAL, {&abin}, |
| 3599 | "Binwidth angle distribution (degrees)" }, |
| 3600 | { "-rbin", FALSE0, etREAL, {&rbin}, |
| 3601 | "Binwidth distance distribution (nm)" }, |
| 3602 | { "-nitacc", FALSE0, etBOOL, {&bNitAcc}, |
| 3603 | "Regard nitrogen atoms as acceptors" }, |
| 3604 | { "-contact", FALSE0, etBOOL, {&bContact}, |
| 3605 | "Do not look for hydrogen bonds, but merely for contacts within the cut-off distance" }, |
| 3606 | { "-shell", FALSE0, etREAL, {&rshell}, |
| 3607 | "when > 0, only calculate hydrogen bonds within # nm shell around " |
| 3608 | "one particle" }, |
| 3609 | { "-fitstart", FALSE0, etREAL, {&fit_start}, |
| 3610 | "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]" }, |
| 3611 | { "-fitend", FALSE0, etREAL, {&fit_end}, |
| 3612 | "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])" }, |
| 3613 | { "-temp", FALSE0, etREAL, {&temp}, |
| 3614 | "Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" }, |
| 3615 | { "-smooth", FALSE0, etREAL, {&smooth_tail_start}, |
| 3616 | "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/[GRK]tau[grk])" }, |
| 3617 | { "-dump", FALSE0, etINT, {&nDump}, |
| 3618 | "Dump the first N hydrogen bond ACFs in a single [TT].xvg[tt] file for debugging" }, |
| 3619 | { "-max_hb", FALSE0, etREAL, {&maxnhb}, |
| 3620 | "Theoretical maximum number of hydrogen bonds used for normalizing HB autocorrelation function. Can be useful in case the program estimates it wrongly" }, |
| 3621 | { "-merge", FALSE0, etBOOL, {&bMerge}, |
| 3622 | "H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." }, |
| 3623 | { "-geminate", FALSE0, etENUM, {gemType}, |
| 3624 | "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."}, |
| 3625 | { "-diff", FALSE0, etREAL, {&D}, |
| 3626 | "HIDDENDffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."}, |
| 3627 | #ifdef GMX_OPENMP |
| 3628 | { "-nthreads", FALSE0, etINT, {&nThreads}, |
| 3629 | "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 processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"}, |
| 3630 | #endif |
| 3631 | }; |
| 3632 | const char *bugs[] = { |
| 3633 | "The option [TT]-sel[tt] that used to work on selected hbonds is out of order, and therefore not available for the time being." |
| 3634 | }; |
| 3635 | t_filenm fnm[] = { |
| 3636 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, |
| 3637 | { efTPX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
| 3638 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
| 3639 | /* { efNDX, "-sel", "select", ffOPTRD },*/ |
| 3640 | { efXVG, "-num", "hbnum", ffWRITE1<<2 }, |
| 3641 | { efLOG, "-g", "hbond", ffOPTWR(1<<2| 1<<3) }, |
| 3642 | { efXVG, "-ac", "hbac", ffOPTWR(1<<2| 1<<3) }, |
| 3643 | { efXVG, "-dist", "hbdist", ffOPTWR(1<<2| 1<<3) }, |
| 3644 | { efXVG, "-ang", "hbang", ffOPTWR(1<<2| 1<<3) }, |
| 3645 | { efXVG, "-hx", "hbhelix", ffOPTWR(1<<2| 1<<3) }, |
| 3646 | { efNDX, "-hbn", "hbond", ffOPTWR(1<<2| 1<<3) }, |
| 3647 | { efXPM, "-hbm", "hbmap", ffOPTWR(1<<2| 1<<3) }, |
| 3648 | { efXVG, "-don", "donor", ffOPTWR(1<<2| 1<<3) }, |
| 3649 | { efXVG, "-dan", "danum", ffOPTWR(1<<2| 1<<3) }, |
| 3650 | { efXVG, "-life", "hblife", ffOPTWR(1<<2| 1<<3) }, |
| 3651 | { efXVG, "-nhbdist", "nhbdist", ffOPTWR(1<<2| 1<<3) } |
| 3652 | |
| 3653 | }; |
| 3654 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
| 3655 | |
| 3656 | char hbmap [HB_NR(1<<2)] = { ' ', 'o', '-', '*' }; |
| 3657 | const char *hbdesc[HB_NR(1<<2)] = { "None", "Present", "Inserted", "Present & Inserted" }; |
| 3658 | t_rgb hbrgb [HB_NR(1<<2)] = { {1, 1, 1}, {1, 0, 0}, {0, 0, 1}, {1, 0, 1} }; |
| 3659 | |
| 3660 | t_trxstatus *status; |
| 3661 | int trrStatus = 1; |
| 3662 | t_topology top; |
| 3663 | t_inputrec ir; |
| 3664 | t_pargs *ppa; |
| 3665 | int npargs, natoms, nframes = 0, shatom; |
| 3666 | int *isize; |
| 3667 | char **grpnames; |
| 3668 | atom_id **index; |
| 3669 | rvec *x, hbox; |
| 3670 | matrix box; |
| 3671 | real t, ccut, dist = 0.0, ang = 0.0; |
| 3672 | double max_nhb, aver_nhb, aver_dist; |
| 3673 | int h = 0, i = 0, j, k = 0, l, start, end, id, ja, ogrp, nsel; |
| 3674 | int xi, yi, zi, ai; |
| 3675 | int xj, yj, zj, aj, xjj, yjj, zjj; |
| 3676 | int xk, yk, zk, ak, xkk, ykk, zkk; |
| 3677 | gmx_bool bSelected, bHBmap, bStop, bTwo, was, bBox, bTric; |
| 3678 | int *adist, *rdist, *aptr, *rprt; |
| 3679 | int grp, nabin, nrbin, bin, resdist, ihb; |
| 3680 | char **leg; |
| 3681 | t_hbdata *hb, *hbptr; |
| 3682 | FILE *fp, *fpins = NULL((void*)0), *fpnhb = NULL((void*)0); |
| 3683 | t_gridcell ***grid; |
| 3684 | t_ncell *icell, *jcell, *kcell; |
| 3685 | ivec ngrid; |
| 3686 | unsigned char *datable; |
| 3687 | output_env_t oenv; |
| 3688 | int gemmode, NN; |
| 3689 | PSTYPEint peri = 0; |
| 3690 | t_E E; |
| 3691 | int ii, jj, hh, actual_nThreads; |
| 3692 | int threadNr = 0; |
| 3693 | gmx_bool bGem, bNN, bParallel; |
| 3694 | t_gemParams *params = NULL((void*)0); |
| 3695 | gmx_bool bEdge_yjj, bEdge_xjj, bOMP; |
| 3696 | |
| 3697 | t_hbdata **p_hb = NULL((void*)0); /* one per thread, then merge after the frame loop */ |
| 3698 | int **p_adist = NULL((void*)0), **p_rdist = NULL((void*)0); /* a histogram for each thread. */ |
| 3699 | |
| 3700 | #ifdef GMX_OPENMP |
| 3701 | bOMP = TRUE1; |
| 3702 | #else |
| 3703 | bOMP = FALSE0; |
| 3704 | #endif |
| 3705 | |
| 3706 | npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))); |
| 3707 | ppa = add_acf_pargs(&npargs, pa); |
| 3708 | |
| 3709 | if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13), NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, |
| 3710 | ppa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, asize(bugs)((int)(sizeof(bugs)/sizeof((bugs)[0]))), bugs, &oenv)) |
| 3711 | { |
| 3712 | return 0; |
| 3713 | } |
| 3714 | |
| 3715 | /* NN-loop? If so, what estimator to use ?*/ |
| 3716 | NN = 1; |
| 3717 | /* Outcommented for now DvdS 2010-07-13 |
| 3718 | while (NN < NN_NR && gmx_strcasecmp(NNtype[0], NNtype[NN])!=0) |
| 3719 | NN++; |
| 3720 | if (NN == NN_NR) |
| 3721 | gmx_fatal(FARGS, "Invalid NN-loop type."); |
| 3722 | */ |
| 3723 | bNN = FALSE0; |
| 3724 | for (i = 2; bNN == FALSE0 && i < NN_NR; i++) |
| 3725 | { |
| 3726 | bNN = bNN || NN == i; |
| 3727 | } |
| 3728 | |
| 3729 | if (NN > NN_NONE && bMerge) |
| 3730 | { |
| 3731 | bMerge = FALSE0; |
| 3732 | } |
| 3733 | |
| 3734 | /* geminate recombination? If so, which flavor? */ |
| 3735 | gemmode = 1; |
| 3736 | while (gemmode < gemNR && gmx_strcasecmp(gemType[0], gemType[gemmode]) != 0) |
| 3737 | { |
| 3738 | gemmode++; |
| 3739 | } |
| 3740 | if (gemmode == gemNR) |
| 3741 | { |
| 3742 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3742, "Invalid recombination type."); |
| 3743 | } |
| 3744 | |
| 3745 | bGem = FALSE0; |
| 3746 | for (i = 2; bGem == FALSE0 && i < gemNR; i++) |
| 3747 | { |
| 3748 | bGem = bGem || gemmode == i; |
| 3749 | } |
| 3750 | |
| 3751 | if (bGem) |
| 3752 | { |
| 3753 | printf("Geminate recombination: %s\n", gemType[gemmode]); |
| 3754 | if (bContact) |
| 3755 | { |
| 3756 | if (gemmode != gemDD) |
| 3757 | { |
| 3758 | printf("Turning off -contact option...\n"); |
| 3759 | bContact = FALSE0; |
| 3760 | } |
| 3761 | } |
| 3762 | else |
| 3763 | { |
| 3764 | if (gemmode == gemDD) |
| 3765 | { |
| 3766 | printf("Turning on -contact option...\n"); |
| 3767 | bContact = TRUE1; |
| 3768 | } |
| 3769 | } |
| 3770 | if (bMerge) |
| 3771 | { |
| 3772 | if (gemmode == gemAA) |
| 3773 | { |
| 3774 | printf("Turning off -merge option...\n"); |
| 3775 | bMerge = FALSE0; |
| 3776 | } |
| 3777 | } |
| 3778 | else |
| 3779 | { |
| 3780 | if (gemmode != gemAA) |
| 3781 | { |
| 3782 | printf("Turning on -merge option...\n"); |
| 3783 | bMerge = TRUE1; |
| 3784 | } |
| 3785 | } |
| 3786 | } |
| 3787 | |
| 3788 | /* process input */ |
| 3789 | bSelected = FALSE0; |
| 3790 | ccut = cos(acut*DEG2RAD(3.14159265358979323846/180.0)); |
| 3791 | |
| 3792 | if (bContact) |
| 3793 | { |
| 3794 | if (bSelected) |
| 3795 | { |
| 3796 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3796, "Can not analyze selected contacts."); |
| 3797 | } |
| 3798 | if (!bDA) |
| 3799 | { |
| 3800 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3800, "Can not analyze contact between H and A: turn off -noda"); |
| 3801 | } |
| 3802 | } |
| 3803 | |
| 3804 | /* Initiate main data structure! */ |
| 3805 | bHBmap = (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
| 3806 | opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
| 3807 | opt2bSet("-hbn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
| 3808 | opt2bSet("-hbm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || |
| 3809 | bGem); |
| 3810 | |
| 3811 | if (opt2bSet("-nhbdist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 3812 | { |
| 3813 | const char *leg[MAXHH4+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" }; |
| 3814 | fpnhb = xvgropen(opt2fn("-nhbdist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 3815 | "Number of donor-H with N HBs", output_env_get_xvgr_tlabel(oenv), "N", oenv); |
| 3816 | xvgr_legend(fpnhb, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); |
| 3817 | } |
| 3818 | |
| 3819 | hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bMerge || bContact, bGem, gemmode); |
| 3820 | |
| 3821 | /* get topology */ |
| 3822 | read_tpx_top(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &ir, box, &natoms, NULL((void*)0), NULL((void*)0), NULL((void*)0), &top); |
| 3823 | |
| 3824 | snew(grpnames, grNR)(grpnames) = save_calloc("grpnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3824, (grNR), sizeof(*(grpnames))); |
| 3825 | snew(index, grNR)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3825, (grNR), sizeof(*(index))); |
| 3826 | snew(isize, grNR)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3826, (grNR), sizeof(*(isize))); |
| 3827 | /* Make Donor-Acceptor table */ |
| 3828 | snew(datable, top.atoms.nr)(datable) = save_calloc("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3828, (top.atoms.nr), sizeof(*(datable))); |
| 3829 | gen_datable(index[0], isize[0], datable, top.atoms.nr); |
| 3830 | |
| 3831 | if (bSelected) |
| 3832 | { |
| 3833 | /* analyze selected hydrogen bonds */ |
| 3834 | printf("Select group with selected atoms:\n"); |
| 3835 | get_index(&(top.atoms), opt2fn("-sel", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 3836 | 1, &nsel, index, grpnames); |
| 3837 | if (nsel % 3) |
| 3838 | { |
| 3839 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3839, "Number of atoms in group '%s' not a multiple of 3\n" |
| 3840 | "and therefore cannot contain triplets of " |
| 3841 | "Donor-Hydrogen-Acceptor", grpnames[0]); |
| 3842 | } |
| 3843 | bTwo = FALSE0; |
| 3844 | |
| 3845 | for (i = 0; (i < nsel); i += 3) |
| 3846 | { |
| 3847 | int dd = index[0][i]; |
| 3848 | int aa = index[0][i+2]; |
| 3849 | /* int */ hh = index[0][i+1]; |
| 3850 | add_dh (&hb->d, dd, hh, i, datable); |
| 3851 | add_acc(&hb->a, aa, i); |
| 3852 | /* Should this be here ? */ |
| 3853 | snew(hb->d.dptr, top.atoms.nr)(hb->d.dptr) = save_calloc("hb->d.dptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3853, (top.atoms.nr), sizeof(*(hb->d.dptr))); |
| 3854 | snew(hb->a.aptr, top.atoms.nr)(hb->a.aptr) = save_calloc("hb->a.aptr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3854, (top.atoms.nr), sizeof(*(hb->a.aptr))); |
| 3855 | add_hbond(hb, dd, aa, hh, gr0, gr0, 0, bMerge, 0, bContact, peri); |
| 3856 | } |
| 3857 | printf("Analyzing %d selected hydrogen bonds from '%s'\n", |
| 3858 | isize[0], grpnames[0]); |
| 3859 | } |
| 3860 | else |
| 3861 | { |
| 3862 | /* analyze all hydrogen bonds: get group(s) */ |
| 3863 | printf("Specify 2 groups to analyze:\n"); |
| 3864 | get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 3865 | 2, isize, index, grpnames); |
| 3866 | |
| 3867 | /* check if we have two identical or two non-overlapping groups */ |
| 3868 | bTwo = isize[0] != isize[1]; |
| 3869 | for (i = 0; (i < isize[0]) && !bTwo; i++) |
| 3870 | { |
| 3871 | bTwo = index[0][i] != index[1][i]; |
| 3872 | } |
| 3873 | if (bTwo) |
| 3874 | { |
| 3875 | printf("Checking for overlap in atoms between %s and %s\n", |
| 3876 | grpnames[0], grpnames[1]); |
| 3877 | for (i = 0; i < isize[1]; i++) |
| 3878 | { |
| 3879 | if (ISINGRP(datable[index[1][i]])(((datable[index[1][i]]) & 4) == 4)) |
| 3880 | { |
| 3881 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3881, "Partial overlap between groups '%s' and '%s'", |
| 3882 | grpnames[0], grpnames[1]); |
| 3883 | } |
| 3884 | } |
| 3885 | /* |
| 3886 | printf("Checking for overlap in atoms between %s and %s\n", |
| 3887 | grpnames[0],grpnames[1]); |
| 3888 | for (i=0; i<isize[0]; i++) |
| 3889 | for (j=0; j<isize[1]; j++) |
| 3890 | if (index[0][i] == index[1][j]) |
| 3891 | gmx_fatal(FARGS,"Partial overlap between groups '%s' and '%s'", |
| 3892 | grpnames[0],grpnames[1]); |
| 3893 | */ |
| 3894 | } |
| 3895 | if (bTwo) |
| 3896 | { |
| 3897 | printf("Calculating %s " |
| 3898 | "between %s (%d atoms) and %s (%d atoms)\n", |
| 3899 | bContact ? "contacts" : "hydrogen bonds", |
| 3900 | grpnames[0], isize[0], grpnames[1], isize[1]); |
| 3901 | } |
| 3902 | else |
| 3903 | { |
| 3904 | fprintf(stderrstderr, "Calculating %s in %s (%d atoms)\n", |
| 3905 | bContact ? "contacts" : "hydrogen bonds", grpnames[0], isize[0]); |
| 3906 | } |
| 3907 | } |
| 3908 | sfree(datable)save_free("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3908, (datable)); |
| 3909 | |
| 3910 | /* search donors and acceptors in groups */ |
| 3911 | snew(datable, top.atoms.nr)(datable) = save_calloc("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3911, (top.atoms.nr), sizeof(*(datable))); |
| 3912 | for (i = 0; (i < grNR); i++) |
| 3913 | { |
| 3914 | if ( ((i == gr0) && !bSelected ) || |
| 3915 | ((i == gr1) && bTwo )) |
| 3916 | { |
| 3917 | gen_datable(index[i], isize[i], datable, top.atoms.nr); |
| 3918 | if (bContact) |
| 3919 | { |
| 3920 | search_acceptors(&top, isize[i], index[i], &hb->a, i, |
| 3921 | bNitAcc, TRUE1, (bTwo && (i == gr0)) || !bTwo, datable); |
| 3922 | search_donors (&top, isize[i], index[i], &hb->d, i, |
| 3923 | TRUE1, (bTwo && (i == gr1)) || !bTwo, datable); |
| 3924 | } |
| 3925 | else |
| 3926 | { |
| 3927 | search_acceptors(&top, isize[i], index[i], &hb->a, i, bNitAcc, FALSE0, TRUE1, datable); |
| 3928 | search_donors (&top, isize[i], index[i], &hb->d, i, FALSE0, TRUE1, datable); |
| 3929 | } |
| 3930 | if (bTwo) |
| 3931 | { |
| 3932 | clear_datable_grp(datable, top.atoms.nr); |
| 3933 | } |
| 3934 | } |
| 3935 | } |
| 3936 | sfree(datable)save_free("datable", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3936, (datable)); |
| 3937 | printf("Found %d donors and %d acceptors\n", hb->d.nrd, hb->a.nra); |
| 3938 | /*if (bSelected) |
| 3939 | snew(donors[gr0D], dons[gr0D].nrd);*/ |
| 3940 | |
| 3941 | if (bHBmap) |
| 3942 | { |
| 3943 | printf("Making hbmap structure..."); |
| 3944 | /* Generate hbond data structure */ |
| 3945 | mk_hbmap(hb); |
| 3946 | printf("done.\n"); |
| 3947 | } |
| 3948 | |
| 3949 | #ifdef HAVE_NN_LOOPS |
| 3950 | if (bNN) |
| 3951 | { |
| 3952 | mk_hbEmap(hb, 0); |
| 3953 | } |
| 3954 | #endif |
| 3955 | |
| 3956 | if (bGem) |
| 3957 | { |
| 3958 | printf("Making per structure..."); |
| 3959 | /* Generate hbond data structure */ |
| 3960 | mk_per(hb); |
| 3961 | printf("done.\n"); |
| 3962 | } |
| 3963 | |
| 3964 | /* check input */ |
| 3965 | bStop = FALSE0; |
| 3966 | if (hb->d.nrd + hb->a.nra == 0) |
| 3967 | { |
| 3968 | printf("No Donors or Acceptors found\n"); |
| 3969 | bStop = TRUE1; |
| 3970 | } |
| 3971 | if (!bStop) |
| 3972 | { |
| 3973 | if (hb->d.nrd == 0) |
| 3974 | { |
| 3975 | printf("No Donors found\n"); |
| 3976 | bStop = TRUE1; |
| 3977 | } |
| 3978 | if (hb->a.nra == 0) |
| 3979 | { |
| 3980 | printf("No Acceptors found\n"); |
| 3981 | bStop = TRUE1; |
| 3982 | } |
| 3983 | } |
| 3984 | if (bStop) |
| 3985 | { |
| 3986 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 3986, "Nothing to be done"); |
| 3987 | } |
| 3988 | |
| 3989 | shatom = 0; |
| 3990 | if (rshell > 0) |
| 3991 | { |
| 3992 | int shisz; |
| 3993 | atom_id *shidx; |
| 3994 | char *shgrpnm; |
| 3995 | /* get index group with atom for shell */ |
| 3996 | do |
| 3997 | { |
| 3998 | printf("Select atom for shell (1 atom):\n"); |
| 3999 | get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 4000 | 1, &shisz, &shidx, &shgrpnm); |
| 4001 | if (shisz != 1) |
| 4002 | { |
| 4003 | printf("group contains %d atoms, should be 1 (one)\n", shisz); |
| 4004 | } |
| 4005 | } |
| 4006 | while (shisz != 1); |
| 4007 | shatom = shidx[0]; |
| 4008 | printf("Will calculate hydrogen bonds within a shell " |
| 4009 | "of %g nm around atom %i\n", rshell, shatom+1); |
| 4010 | } |
| 4011 | |
| 4012 | /* Analyze trajectory */ |
| 4013 | natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box); |
| 4014 | if (natoms > top.atoms.nr) |
| 4015 | { |
| 4016 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4016, "Topology (%d atoms) does not match trajectory (%d atoms)", |
| 4017 | top.atoms.nr, natoms); |
| 4018 | } |
| 4019 | |
| 4020 | bBox = ir.ePBC != epbcNONE; |
| 4021 | grid = init_grid(bBox, box, (rcut > r2cut) ? rcut : r2cut, ngrid); |
| 4022 | nabin = acut/abin; |
| 4023 | nrbin = rcut/rbin; |
| 4024 | snew(adist, nabin+1)(adist) = save_calloc("adist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4024, (nabin+1), sizeof(*(adist))); |
| 4025 | snew(rdist, nrbin+1)(rdist) = save_calloc("rdist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4025, (nrbin+1), sizeof(*(rdist))); |
| 4026 | |
| 4027 | if (bGem && !bBox) |
| 4028 | { |
| 4029 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4029, "Can't do geminate recombination without periodic box."); |
| 4030 | } |
| 4031 | |
| 4032 | bParallel = FALSE0; |
| 4033 | |
| 4034 | #ifndef GMX_OPENMP |
| 4035 | #define __ADISTadist adist |
| 4036 | #define __RDISTrdist rdist |
| 4037 | #define __HBDATAhb hb |
| 4038 | #else /* GMX_OPENMP ================================================== \ |
| 4039 | * Set up the OpenMP stuff, | |
| 4040 | * like the number of threads and such | |
| 4041 | * Also start the parallel loop. | |
| 4042 | */ |
| 4043 | #define __ADISTadist p_adist[threadNr] |
| 4044 | #define __RDISTrdist p_rdist[threadNr] |
| 4045 | #define __HBDATAhb p_hb[threadNr] |
| 4046 | #endif |
| 4047 | if (bOMP) |
| 4048 | { |
| 4049 | bParallel = !bSelected; |
| 4050 | |
| 4051 | if (bParallel) |
| 4052 | { |
| 4053 | actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads())((((nThreads <= 0) ? 2147483647 : nThreads) < (gmx_omp_get_max_threads ())) ? ((nThreads <= 0) ? 2147483647 : nThreads) : (gmx_omp_get_max_threads ()) ); |
| 4054 | |
| 4055 | gmx_omp_set_num_threads(actual_nThreads); |
| 4056 | printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads); |
| 4057 | fflush(stdoutstdout); |
| 4058 | } |
| 4059 | else |
| 4060 | { |
| 4061 | actual_nThreads = 1; |
| 4062 | } |
| 4063 | |
| 4064 | snew(p_hb, actual_nThreads)(p_hb) = save_calloc("p_hb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4064, (actual_nThreads), sizeof(*(p_hb))); |
| 4065 | snew(p_adist, actual_nThreads)(p_adist) = save_calloc("p_adist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4065, (actual_nThreads), sizeof(*(p_adist))); |
| 4066 | snew(p_rdist, actual_nThreads)(p_rdist) = save_calloc("p_rdist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4066, (actual_nThreads), sizeof(*(p_rdist))); |
| 4067 | for (i = 0; i < actual_nThreads; i++) |
| 4068 | { |
| 4069 | snew(p_hb[i], 1)(p_hb[i]) = save_calloc("p_hb[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4069, (1), sizeof(*(p_hb[i]))); |
| 4070 | snew(p_adist[i], nabin+1)(p_adist[i]) = save_calloc("p_adist[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4070, (nabin+1), sizeof(*(p_adist[i]))); |
| 4071 | snew(p_rdist[i], nrbin+1)(p_rdist[i]) = save_calloc("p_rdist[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4071, (nrbin+1), sizeof(*(p_rdist[i]))); |
| 4072 | |
| 4073 | p_hb[i]->max_frames = 0; |
| 4074 | p_hb[i]->nhb = NULL((void*)0); |
| 4075 | p_hb[i]->ndist = NULL((void*)0); |
| 4076 | p_hb[i]->n_bound = NULL((void*)0); |
| 4077 | p_hb[i]->time = NULL((void*)0); |
| 4078 | p_hb[i]->nhx = NULL((void*)0); |
| 4079 | |
| 4080 | p_hb[i]->bHBmap = hb->bHBmap; |
| 4081 | p_hb[i]->bDAnr = hb->bDAnr; |
| 4082 | p_hb[i]->bGem = hb->bGem; |
| 4083 | p_hb[i]->wordlen = hb->wordlen; |
| 4084 | p_hb[i]->nframes = hb->nframes; |
| 4085 | p_hb[i]->maxhydro = hb->maxhydro; |
| 4086 | p_hb[i]->danr = hb->danr; |
| 4087 | p_hb[i]->d = hb->d; |
| 4088 | p_hb[i]->a = hb->a; |
| 4089 | p_hb[i]->hbmap = hb->hbmap; |
| 4090 | p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */ |
| 4091 | p_hb[i]->per = hb->per; |
| 4092 | |
| 4093 | #ifdef HAVE_NN_LOOPS |
| 4094 | p_hb[i]->hbE = hb->hbE; |
| 4095 | #endif |
| 4096 | |
| 4097 | p_hb[i]->nrhb = 0; |
| 4098 | p_hb[i]->nrdist = 0; |
| 4099 | } |
| 4100 | } |
| 4101 | |
| 4102 | /* Make a thread pool here, |
| 4103 | * instead of forking anew at every frame. */ |
| 4104 | |
| 4105 | #pragma omp parallel \ |
| 4106 | firstprivate(i) \ |
| 4107 | private(j, h, ii, jj, hh, E, \ |
| 4108 | xi, yi, zi, xj, yj, zj, threadNr, \ |
| 4109 | dist, ang, peri, icell, jcell, \ |
| 4110 | grp, ogrp, ai, aj, xjj, yjj, zjj, \ |
| 4111 | xk, yk, zk, ihb, id, resdist, \ |
| 4112 | xkk, ykk, zkk, kcell, ak, k, bTric, \ |
| 4113 | bEdge_xjj, bEdge_yjj) \ |
| 4114 | default(shared) |
| 4115 | { /* Start of parallel region */ |
| 4116 | threadNr = gmx_omp_get_thread_num(); |
| 4117 | |
| 4118 | do |
| 4119 | { |
| 4120 | |
| 4121 | bTric = bBox && TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0); |
| 4122 | |
| 4123 | if (bOMP) |
| 4124 | { |
| 4125 | sync_hbdata(p_hb[threadNr], nframes); |
| 4126 | } |
| 4127 | #pragma omp single |
| 4128 | { |
| 4129 | build_grid(hb, x, x[shatom], bBox, box, hbox, (rcut > r2cut) ? rcut : r2cut, |
| 4130 | rshell, ngrid, grid); |
| 4131 | reset_nhbonds(&(hb->d)); |
| 4132 | |
| 4133 | if (debug && bDebug) |
| 4134 | { |
| 4135 | dump_grid(debug, ngrid, grid); |
| 4136 | } |
| 4137 | |
| 4138 | add_frames(hb, nframes); |
| 4139 | init_hbframe(hb, nframes, output_env_conv_time(oenv, t)); |
| 4140 | |
| 4141 | if (hb->bDAnr) |
| 4142 | { |
| 4143 | count_da_grid(ngrid, grid, hb->danr[nframes]); |
| 4144 | } |
| 4145 | } /* omp single */ |
| 4146 | |
| 4147 | if (bOMP) |
| 4148 | { |
| 4149 | p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */ |
| 4150 | } |
| 4151 | |
| 4152 | if (bNN) |
| 4153 | { |
| 4154 | #ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */ |
| 4155 | /* Loop over all atom pairs and estimate interaction energy */ |
| 4156 | |
| 4157 | #pragma omp single |
| 4158 | { |
| 4159 | addFramesNN(hb, nframes); |
| 4160 | } |
| 4161 | |
| 4162 | #pragma omp barrier |
| 4163 | #pragma omp for schedule(dynamic) |
| 4164 | for (i = 0; i < hb->d.nrd; i++) |
| 4165 | { |
| 4166 | for (j = 0; j < hb->a.nra; j++) |
| 4167 | { |
| 4168 | for (h = 0; |
| 4169 | h < (bContact ? 1 : hb->d.nhydro[i]); |
| 4170 | h++) |
| 4171 | { |
| 4172 | if (i == hb->d.nrd || j == hb->a.nra) |
| 4173 | { |
| 4174 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4174, "out of bounds"); |
| 4175 | } |
| 4176 | |
| 4177 | /* Get the real atom ids */ |
| 4178 | ii = hb->d.don[i]; |
| 4179 | jj = hb->a.acc[j]; |
| 4180 | hh = hb->d.hydro[i][h]; |
| 4181 | |
| 4182 | /* Estimate the energy from the geometry */ |
| 4183 | E = calcHbEnergy(ii, jj, hh, x, NN, box, hbox, &(hb->d)); |
| 4184 | /* Store the energy */ |
| 4185 | storeHbEnergy(hb, i, j, h, E, nframes); |
| 4186 | } |
| 4187 | } |
| 4188 | } |
| 4189 | #endif /* HAVE_NN_LOOPS */ |
| 4190 | } /* if (bNN)*/ |
| 4191 | else |
| 4192 | { |
| 4193 | if (bSelected) |
| 4194 | { |
| 4195 | |
| 4196 | #pragma omp single |
| 4197 | { |
| 4198 | /* Do not parallelize this just yet. */ |
| 4199 | /* int ii; */ |
| 4200 | for (ii = 0; (ii < nsel); ii++) |
| 4201 | { |
| 4202 | int dd = index[0][i]; |
| 4203 | int aa = index[0][i+2]; |
| 4204 | /* int */ hh = index[0][i+1]; |
| 4205 | ihb = is_hbond(hb, ii, ii, dd, aa, rcut, r2cut, ccut, x, bBox, box, |
| 4206 | hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); |
| 4207 | |
| 4208 | if (ihb) |
| 4209 | { |
| 4210 | /* add to index if not already there */ |
| 4211 | /* Add a hbond */ |
| 4212 | add_hbond(hb, dd, aa, hh, ii, ii, nframes, bMerge, ihb, bContact, peri); |
| 4213 | } |
| 4214 | } |
| 4215 | } /* omp single */ |
| 4216 | } /* if (bSelected) */ |
| 4217 | else |
| 4218 | { |
| 4219 | |
| 4220 | #pragma omp single |
| 4221 | { |
| 4222 | if (bGem) |
| 4223 | { |
| 4224 | calcBoxProjection(box, hb->per->P); |
| 4225 | } |
| 4226 | |
| 4227 | /* loop over all gridcells (xi,yi,zi) */ |
| 4228 | /* Removed confusing macro, DvdS 27/12/98 */ |
| 4229 | |
| 4230 | } |
| 4231 | /* The outer grid loop will have to do for now. */ |
| 4232 | #pragma omp for schedule(dynamic) |
| 4233 | for (xi = 0; xi < ngrid[XX0]; xi++) |
| 4234 | { |
| 4235 | for (yi = 0; (yi < ngrid[YY1]); yi++) |
| 4236 | { |
| 4237 | for (zi = 0; (zi < ngrid[ZZ2]); zi++) |
| 4238 | { |
| 4239 | |
| 4240 | /* loop over donor groups gr0 (always) and gr1 (if necessary) */ |
| 4241 | for (grp = gr0; (grp <= (bTwo ? gr1 : gr0)); grp++) |
| 4242 | { |
| 4243 | icell = &(grid[zi][yi][xi].d[grp]); |
| 4244 | |
| 4245 | if (bTwo) |
| 4246 | { |
| 4247 | ogrp = 1-grp; |
| 4248 | } |
| 4249 | else |
| 4250 | { |
| 4251 | ogrp = grp; |
| 4252 | } |
| 4253 | |
| 4254 | /* loop over all hydrogen atoms from group (grp) |
| 4255 | * in this gridcell (icell) |
| 4256 | */ |
| 4257 | for (ai = 0; (ai < icell->nr); ai++) |
| 4258 | { |
| 4259 | i = icell->atoms[ai]; |
| 4260 | |
| 4261 | /* loop over all adjacent gridcells (xj,yj,zj) */ |
| 4262 | for (zjj = grid_loop_begin(ngrid[ZZ2], zi, bTric, FALSE0); |
| 4263 | zjj <= grid_loop_end(ngrid[ZZ2], zi, bTric, FALSE0); |
| 4264 | zjj++) |
| 4265 | { |
| 4266 | zj = grid_mod(zjj, ngrid[ZZ2]); |
| 4267 | bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ2] - 1); |
| 4268 | for (yjj = grid_loop_begin(ngrid[YY1], yi, bTric, bEdge_yjj); |
| 4269 | yjj <= grid_loop_end(ngrid[YY1], yi, bTric, bEdge_yjj); |
| 4270 | yjj++) |
| 4271 | { |
| 4272 | yj = grid_mod(yjj, ngrid[YY1]); |
| 4273 | bEdge_xjj = |
| 4274 | (yj == 0) || (yj == ngrid[YY1] - 1) || |
| 4275 | (zj == 0) || (zj == ngrid[ZZ2] - 1); |
| 4276 | for (xjj = grid_loop_begin(ngrid[XX0], xi, bTric, bEdge_xjj); |
| 4277 | xjj <= grid_loop_end(ngrid[XX0], xi, bTric, bEdge_xjj); |
| 4278 | xjj++) |
| 4279 | { |
| 4280 | xj = grid_mod(xjj, ngrid[XX0]); |
| 4281 | jcell = &(grid[zj][yj][xj].a[ogrp]); |
| 4282 | /* loop over acceptor atoms from other group (ogrp) |
| 4283 | * in this adjacent gridcell (jcell) |
| 4284 | */ |
| 4285 | for (aj = 0; (aj < jcell->nr); aj++) |
| 4286 | { |
| 4287 | j = jcell->atoms[aj]; |
| 4288 | |
| 4289 | /* check if this once was a h-bond */ |
| 4290 | peri = -1; |
| 4291 | ihb = is_hbond(__HBDATAhb, grp, ogrp, i, j, rcut, r2cut, ccut, x, bBox, box, |
| 4292 | hbox, &dist, &ang, bDA, &h, bContact, bMerge, &peri); |
| 4293 | |
| 4294 | if (ihb) |
| 4295 | { |
| 4296 | /* add to index if not already there */ |
| 4297 | /* Add a hbond */ |
| 4298 | add_hbond(__HBDATAhb, i, j, h, grp, ogrp, nframes, bMerge, ihb, bContact, peri); |
| 4299 | |
| 4300 | /* make angle and distance distributions */ |
| 4301 | if (ihb == hbHB && !bContact) |
| 4302 | { |
| 4303 | if (dist > rcut) |
| 4304 | { |
| 4305 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4305, "distance is higher than what is allowed for an hbond: %f", dist); |
| 4306 | } |
| 4307 | ang *= RAD2DEG(180.0/3.14159265358979323846); |
| 4308 | __ADISTadist[(int)( ang/abin)]++; |
| 4309 | __RDISTrdist[(int)(dist/rbin)]++; |
| 4310 | if (!bTwo) |
| 4311 | { |
| 4312 | int id, ia; |
| 4313 | if ((id = donor_index(&hb->d, grp, i)_donor_index(&hb->d, grp, i, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4313)) == NOTSET-12345) |
| 4314 | { |
| 4315 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4315, "Invalid donor %d", i); |
| 4316 | } |
| 4317 | if ((ia = acceptor_index(&hb->a, ogrp, j)_acceptor_index(&hb->a, ogrp, j, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4317)) == NOTSET-12345) |
| 4318 | { |
| 4319 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4319, "Invalid acceptor %d", j); |
| 4320 | } |
| 4321 | resdist = abs(top.atoms.atom[i].resind- |
| 4322 | top.atoms.atom[j].resind); |
| 4323 | if (resdist >= max_hx7) |
| 4324 | { |
| 4325 | resdist = max_hx7-1; |
| 4326 | } |
| 4327 | __HBDATAhb->nhx[nframes][resdist]++; |
| 4328 | } |
| 4329 | } |
| 4330 | |
| 4331 | } |
| 4332 | } /* for aj */ |
| 4333 | } /* for xjj */ |
| 4334 | } /* for yjj */ |
| 4335 | } /* for zjj */ |
| 4336 | } /* for ai */ |
| 4337 | } /* for grp */ |
| 4338 | } /* for xi,yi,zi */ |
| 4339 | } |
| 4340 | } |
| 4341 | } /* if (bSelected) {...} else */ |
| 4342 | |
| 4343 | |
| 4344 | /* Better wait for all threads to finnish using x[] before updating it. */ |
| 4345 | k = nframes; |
| 4346 | #pragma omp barrier |
| 4347 | #pragma omp critical |
| 4348 | { |
| 4349 | /* Sum up histograms and counts from p_hb[] into hb */ |
| 4350 | if (bOMP) |
| 4351 | { |
| 4352 | hb->nhb[k] += p_hb[threadNr]->nhb[k]; |
| 4353 | hb->ndist[k] += p_hb[threadNr]->ndist[k]; |
| 4354 | for (j = 0; j < max_hx7; j++) |
| 4355 | { |
| 4356 | hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j]; |
| 4357 | } |
| 4358 | } |
| 4359 | } |
| 4360 | |
| 4361 | /* Here are a handful of single constructs |
| 4362 | * to share the workload a bit. The most |
| 4363 | * important one is of course the last one, |
| 4364 | * where there's a potential bottleneck in form |
| 4365 | * of slow I/O. */ |
| 4366 | #pragma omp barrier |
| 4367 | #pragma omp single |
| 4368 | { |
| 4369 | if (hb != NULL((void*)0)) |
| 4370 | { |
| 4371 | analyse_donor_props(opt2fn_null("-don", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), hb, k, t, oenv); |
| 4372 | } |
| 4373 | } |
| 4374 | |
| 4375 | #pragma omp single |
| 4376 | { |
| 4377 | if (fpnhb) |
| 4378 | { |
| 4379 | do_nhb_dist(fpnhb, hb, t); |
| 4380 | } |
| 4381 | } |
| 4382 | } /* if (bNN) {...} else + */ |
| 4383 | |
| 4384 | #pragma omp single |
| 4385 | { |
| 4386 | trrStatus = (read_next_x(oenv, status, &t, x, box)); |
| 4387 | nframes++; |
| 4388 | } |
| 4389 | |
| 4390 | #pragma omp barrier |
| 4391 | } |
| 4392 | while (trrStatus); |
| 4393 | |
| 4394 | if (bOMP) |
| 4395 | { |
| 4396 | #pragma omp critical |
| 4397 | { |
| 4398 | hb->nrhb += p_hb[threadNr]->nrhb; |
| 4399 | hb->nrdist += p_hb[threadNr]->nrdist; |
| 4400 | } |
| 4401 | /* Free parallel datastructures */ |
| 4402 | sfree(p_hb[threadNr]->nhb)save_free("p_hb[threadNr]->nhb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4402, (p_hb[threadNr]->nhb)); |
| 4403 | sfree(p_hb[threadNr]->ndist)save_free("p_hb[threadNr]->ndist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4403, (p_hb[threadNr]->ndist)); |
| 4404 | sfree(p_hb[threadNr]->nhx)save_free("p_hb[threadNr]->nhx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4404, (p_hb[threadNr]->nhx)); |
| 4405 | |
| 4406 | #pragma omp for |
| 4407 | for (i = 0; i < nabin; i++) |
| 4408 | { |
| 4409 | for (j = 0; j < actual_nThreads; j++) |
| 4410 | { |
| 4411 | |
| 4412 | adist[i] += p_adist[j][i]; |
| 4413 | } |
| 4414 | } |
| 4415 | #pragma omp for |
| 4416 | for (i = 0; i <= nrbin; i++) |
| 4417 | { |
| 4418 | for (j = 0; j < actual_nThreads; j++) |
| 4419 | { |
| 4420 | rdist[i] += p_rdist[j][i]; |
| 4421 | } |
| 4422 | } |
| 4423 | |
| 4424 | sfree(p_adist[threadNr])save_free("p_adist[threadNr]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4424, (p_adist[threadNr])); |
| 4425 | sfree(p_rdist[threadNr])save_free("p_rdist[threadNr]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4425, (p_rdist[threadNr])); |
| 4426 | } |
| 4427 | } /* End of parallel region */ |
| 4428 | if (bOMP) |
| 4429 | { |
| 4430 | sfree(p_adist)save_free("p_adist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4430, (p_adist)); |
| 4431 | sfree(p_rdist)save_free("p_rdist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4431, (p_rdist)); |
| 4432 | } |
| 4433 | |
| 4434 | if (nframes < 2 && (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))) |
| 4435 | { |
| 4436 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4436, "Cannot calculate autocorrelation of life times with less than two frames"); |
| 4437 | } |
| 4438 | |
| 4439 | free_grid(ngrid, &grid); |
| 4440 | |
| 4441 | close_trj(status); |
| 4442 | if (fpnhb) |
| 4443 | { |
| 4444 | gmx_ffclose(fpnhb); |
| 4445 | } |
| 4446 | |
| 4447 | /* Compute maximum possible number of different hbonds */ |
| 4448 | if (maxnhb > 0) |
| 4449 | { |
| 4450 | max_nhb = maxnhb; |
| 4451 | } |
| 4452 | else |
| 4453 | { |
| 4454 | max_nhb = 0.5*(hb->d.nrd*hb->a.nra); |
| 4455 | } |
| 4456 | /* Added support for -contact below. |
| 4457 | * - Erik Marklund, May 29-31, 2006 */ |
| 4458 | /* Changed contact code. |
| 4459 | * - Erik Marklund, June 29, 2006 */ |
| 4460 | if (bHBmap && !bNN) |
| 4461 | { |
| 4462 | if (hb->nrhb == 0) |
| 4463 | { |
| 4464 | printf("No %s found!!\n", bContact ? "contacts" : "hydrogen bonds"); |
| 4465 | } |
| 4466 | else |
| 4467 | { |
| 4468 | printf("Found %d different %s in trajectory\n" |
| 4469 | "Found %d different atom-pairs within %s distance\n", |
| 4470 | hb->nrhb, bContact ? "contacts" : "hydrogen bonds", |
| 4471 | hb->nrdist, (r2cut > 0) ? "second cut-off" : "hydrogen bonding"); |
| 4472 | |
| 4473 | /*Control the pHist.*/ |
| 4474 | |
| 4475 | if (bMerge) |
| 4476 | { |
| 4477 | merge_hb(hb, bTwo, bContact); |
| 4478 | } |
| 4479 | |
| 4480 | if (opt2bSet("-hbn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4481 | { |
| 4482 | dump_hbmap(hb, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, bTwo, bContact, isize, index, grpnames, &top.atoms); |
| 4483 | } |
| 4484 | |
| 4485 | /* Moved the call to merge_hb() to a line BEFORE dump_hbmap |
| 4486 | * to make the -hbn and -hmb output match eachother. |
| 4487 | * - Erik Marklund, May 30, 2006 */ |
| 4488 | } |
| 4489 | } |
| 4490 | /* Print out number of hbonds and distances */ |
| 4491 | aver_nhb = 0; |
| 4492 | aver_dist = 0; |
| 4493 | fp = xvgropen(opt2fn("-num", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bContact ? "Contacts" : |
| 4494 | "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Number", oenv); |
| 4495 | snew(leg, 2)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4495, (2), sizeof(*(leg))); |
| 4496 | snew(leg[0], STRLEN)(leg[0]) = save_calloc("leg[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4496, (4096), sizeof(*(leg[0]))); |
| 4497 | snew(leg[1], STRLEN)(leg[1]) = save_calloc("leg[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4497, (4096), sizeof(*(leg[1]))); |
| 4498 | sprintf(leg[0], "%s", bContact ? "Contacts" : "Hydrogen bonds"); |
| 4499 | sprintf(leg[1], "Pairs within %g nm", (r2cut > 0) ? r2cut : rcut); |
| 4500 | xvgr_legend(fp, 2, (const char**)leg, oenv); |
| 4501 | sfree(leg[1])save_free("leg[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4501, (leg[1])); |
| 4502 | sfree(leg[0])save_free("leg[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4502, (leg[0])); |
| 4503 | sfree(leg)save_free("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4503, (leg)); |
| 4504 | for (i = 0; (i < nframes); i++) |
| 4505 | { |
| 4506 | fprintf(fp, "%10g %10d %10d\n", hb->time[i], hb->nhb[i], hb->ndist[i]); |
| 4507 | aver_nhb += hb->nhb[i]; |
| 4508 | aver_dist += hb->ndist[i]; |
| 4509 | } |
| 4510 | gmx_ffclose(fp); |
| 4511 | aver_nhb /= nframes; |
| 4512 | aver_dist /= nframes; |
| 4513 | /* Print HB distance distribution */ |
| 4514 | if (opt2bSet("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4515 | { |
| 4516 | long sum; |
| 4517 | |
| 4518 | sum = 0; |
| 4519 | for (i = 0; i < nrbin; i++) |
| 4520 | { |
| 4521 | sum += rdist[i]; |
| 4522 | } |
| 4523 | |
| 4524 | fp = xvgropen(opt2fn("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 4525 | "Hydrogen Bond Distribution", |
| 4526 | bDA ? |
| 4527 | "Donor - Acceptor Distance (nm)" : |
| 4528 | "Hydrogen - Acceptor Distance (nm)", "", oenv); |
| 4529 | for (i = 0; i < nrbin; i++) |
| 4530 | { |
| 4531 | fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum)); |
| 4532 | } |
| 4533 | gmx_ffclose(fp); |
| 4534 | } |
| 4535 | |
| 4536 | /* Print HB angle distribution */ |
| 4537 | if (opt2bSet("-ang", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4538 | { |
| 4539 | long sum; |
| 4540 | |
| 4541 | sum = 0; |
| 4542 | for (i = 0; i < nabin; i++) |
| 4543 | { |
| 4544 | sum += adist[i]; |
| 4545 | } |
| 4546 | |
| 4547 | fp = xvgropen(opt2fn("-ang", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 4548 | "Hydrogen Bond Distribution", |
| 4549 | "Hydrogen - Donor - Acceptor Angle (\\SO\\N)", "", oenv); |
| 4550 | for (i = 0; i < nabin; i++) |
| 4551 | { |
| 4552 | fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum)); |
| 4553 | } |
| 4554 | gmx_ffclose(fp); |
| 4555 | } |
| 4556 | |
| 4557 | /* Print HB in alpha-helix */ |
| 4558 | if (opt2bSet("-hx", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4559 | { |
| 4560 | fp = xvgropen(opt2fn("-hx", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 4561 | "Hydrogen Bonds", output_env_get_xvgr_tlabel(oenv), "Count", oenv); |
| 4562 | xvgr_legend(fp, NRHXTYPES7, hxtypenames, oenv); |
| 4563 | for (i = 0; i < nframes; i++) |
| 4564 | { |
| 4565 | fprintf(fp, "%10g", hb->time[i]); |
| 4566 | for (j = 0; j < max_hx7; j++) |
| 4567 | { |
| 4568 | fprintf(fp, " %6d", hb->nhx[i][j]); |
| 4569 | } |
| 4570 | fprintf(fp, "\n"); |
| 4571 | } |
| 4572 | gmx_ffclose(fp); |
| 4573 | } |
| 4574 | if (!bNN) |
| 4575 | { |
| 4576 | printf("Average number of %s per timeframe %.3f out of %g possible\n", |
| 4577 | bContact ? "contacts" : "hbonds", |
| 4578 | bContact ? aver_dist : aver_nhb, max_nhb); |
| 4579 | } |
| 4580 | |
| 4581 | /* Do Autocorrelation etc. */ |
| 4582 | if (hb->bHBmap) |
| 4583 | { |
| 4584 | /* |
| 4585 | Added support for -contact in ac and hbm calculations below. |
| 4586 | - Erik Marklund, May 29, 2006 |
| 4587 | */ |
| 4588 | ivec itmp; |
| 4589 | rvec rtmp; |
| 4590 | if (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4591 | { |
| 4592 | please_cite(stdoutstdout, "Spoel2006b"); |
| 4593 | } |
| 4594 | if (opt2bSet("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4595 | { |
| 4596 | char *gemstring = NULL((void*)0); |
| 4597 | |
| 4598 | if (bGem || bNN) |
| 4599 | { |
| 4600 | params = init_gemParams(rcut, D, hb->time, hb->nframes/2, nFitPoints, fit_start, fit_end, |
| 4601 | gemBallistic, nBalExp); |
| 4602 | if (params == NULL((void*)0)) |
| 4603 | { |
| 4604 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4604, "Could not initiate t_gemParams params."); |
| 4605 | } |
| 4606 | } |
| 4607 | gemstring = strdup(gemType[hb->per->gemtype])(__extension__ (__builtin_constant_p (gemType[hb->per-> gemtype]) && ((size_t)(const void *)((gemType[hb-> per->gemtype]) + 1) - (size_t)(const void *)(gemType[hb-> per->gemtype]) == 1) ? (((const char *) (gemType[hb->per ->gemtype]))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t ) 1) : ({ size_t __len = strlen (gemType[hb->per->gemtype ]) + 1; char *__retval = (char *) malloc (__len); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval, gemType [hb->per->gemtype], __len); __retval; })) : __strdup (gemType [hb->per->gemtype]))); |
| 4608 | do_hbac(opt2fn("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), hb, nDump, |
| 4609 | bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv, |
| 4610 | gemstring, nThreads, NN, bBallistic, bGemFit); |
| 4611 | } |
| 4612 | if (opt2bSet("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4613 | { |
| 4614 | do_hblife(opt2fn("-life", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), hb, bMerge, bContact, oenv); |
| 4615 | } |
| 4616 | if (opt2bSet("-hbm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
| 4617 | { |
| 4618 | t_matrix mat; |
| 4619 | int id, ia, hh, x, y; |
| 4620 | |
| 4621 | if ((nframes > 0) && (hb->nrhb > 0)) |
| 4622 | { |
| 4623 | mat.nx = nframes; |
| 4624 | mat.ny = hb->nrhb; |
| 4625 | |
| 4626 | snew(mat.matrix, mat.nx)(mat.matrix) = save_calloc("mat.matrix", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4626, (mat.nx), sizeof(*(mat.matrix))); |
| 4627 | for (x = 0; (x < mat.nx); x++) |
| 4628 | { |
| 4629 | snew(mat.matrix[x], mat.ny)(mat.matrix[x]) = save_calloc("mat.matrix[x]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4629, (mat.ny), sizeof(*(mat.matrix[x]))); |
| 4630 | } |
| 4631 | y = 0; |
| 4632 | for (id = 0; (id < hb->d.nrd); id++) |
| 4633 | { |
| 4634 | for (ia = 0; (ia < hb->a.nra); ia++) |
| 4635 | { |
| 4636 | for (hh = 0; (hh < hb->maxhydro); hh++) |
| 4637 | { |
| 4638 | if (hb->hbmap[id][ia]) |
| 4639 | { |
| 4640 | if (ISHB(hb->hbmap[id][ia]->history[hh])(((hb->hbmap[id][ia]->history[hh]) & 2) == 2)) |
| 4641 | { |
| 4642 | /* Changed '<' into '<=' in the for-statement below. |
| 4643 | * It fixed the previously undiscovered bug that caused |
| 4644 | * the last occurance of an hbond/contact to not be |
| 4645 | * set in mat.matrix. Have a look at any old -hbm-output |
| 4646 | * and you will notice that the last column is allways empty. |
| 4647 | * - Erik Marklund May 30, 2006 |
| 4648 | */ |
| 4649 | for (x = 0; (x <= hb->hbmap[id][ia]->nframes); x++) |
| 4650 | { |
| 4651 | int nn0 = hb->hbmap[id][ia]->n0; |
| 4652 | range_check(y, 0, mat.ny)_range_check(y, 0, mat.ny, ((void*)0),"y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4652); |
| 4653 | mat.matrix[x+nn0][y] = is_hb(hb->hbmap[id][ia]->h[hh], x); |
| 4654 | } |
| 4655 | y++; |
| 4656 | } |
| 4657 | } |
| 4658 | } |
| 4659 | } |
| 4660 | } |
| 4661 | mat.axis_x = hb->time; |
| 4662 | snew(mat.axis_y, mat.ny)(mat.axis_y) = save_calloc("mat.axis_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4662, (mat.ny), sizeof(*(mat.axis_y))); |
| 4663 | for (j = 0; j < mat.ny; j++) |
| 4664 | { |
| 4665 | mat.axis_y[j] = j; |
| 4666 | } |
| 4667 | sprintf(mat.title, bContact ? "Contact Existence Map" : |
| 4668 | "Hydrogen Bond Existence Map"); |
| 4669 | sprintf(mat.legend, bContact ? "Contacts" : "Hydrogen Bonds"); |
| 4670 | sprintf(mat.label_x, "%s", output_env_get_xvgr_tlabel(oenv)); |
| 4671 | sprintf(mat.label_y, bContact ? "Contact Index" : "Hydrogen Bond Index"); |
| 4672 | mat.bDiscrete = TRUE1; |
| 4673 | mat.nmap = 2; |
| 4674 | snew(mat.map, mat.nmap)(mat.map) = save_calloc("mat.map", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4674, (mat.nmap), sizeof(*(mat.map))); |
| 4675 | for (i = 0; i < mat.nmap; i++) |
| 4676 | { |
| 4677 | mat.map[i].code.c1 = hbmap[i]; |
| 4678 | mat.map[i].desc = hbdesc[i]; |
| 4679 | mat.map[i].rgb = hbrgb[i]; |
| 4680 | } |
| 4681 | fp = opt2FILE("-hbm", NFILE, fnm, "w")gmx_ffopen(opt2fn("-hbm", ((int)(sizeof(fnm)/sizeof((fnm)[0]) )), fnm), "w"); |
| 4682 | write_xpm_m(fp, mat); |
| 4683 | gmx_ffclose(fp); |
| 4684 | for (x = 0; x < mat.nx; x++) |
| 4685 | { |
| 4686 | sfree(mat.matrix[x])save_free("mat.matrix[x]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4686, (mat.matrix[x])); |
| 4687 | } |
| 4688 | sfree(mat.axis_y)save_free("mat.axis_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4688, (mat.axis_y)); |
| 4689 | sfree(mat.matrix)save_free("mat.matrix", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4689, (mat.matrix)); |
| 4690 | sfree(mat.map)save_free("mat.map", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4690, (mat.map)); |
| 4691 | } |
| 4692 | else |
| 4693 | { |
| 4694 | fprintf(stderrstderr, "No hydrogen bonds/contacts found. No hydrogen bond map will be printed.\n"); |
| 4695 | } |
| 4696 | } |
| 4697 | } |
| 4698 | |
| 4699 | if (bGem) |
| 4700 | { |
| 4701 | fprintf(stderrstderr, "There were %i periodic shifts\n", hb->per->nper); |
| 4702 | fprintf(stderrstderr, "Freeing pHist for all donors...\n"); |
| 4703 | for (i = 0; i < hb->d.nrd; i++) |
| 4704 | { |
| 4705 | fprintf(stderrstderr, "\r%i", i); |
| 4706 | if (hb->per->pHist[i] != NULL((void*)0)) |
| 4707 | { |
| 4708 | for (j = 0; j < hb->a.nra; j++) |
| 4709 | { |
| 4710 | clearPshift(&(hb->per->pHist[i][j])); |
| 4711 | } |
| 4712 | sfree(hb->per->pHist[i])save_free("hb->per->pHist[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4712, (hb->per->pHist[i])); |
| 4713 | } |
| 4714 | } |
| 4715 | sfree(hb->per->pHist)save_free("hb->per->pHist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4715, (hb->per->pHist)); |
| 4716 | sfree(hb->per->p2i)save_free("hb->per->p2i", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4716, (hb->per->p2i)); |
| 4717 | sfree(hb->per)save_free("hb->per", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4717, (hb->per)); |
| 4718 | fprintf(stderrstderr, "...done.\n"); |
| 4719 | } |
| 4720 | |
| 4721 | #ifdef HAVE_NN_LOOPS |
| 4722 | if (bNN) |
| 4723 | { |
| 4724 | free_hbEmap(hb); |
| 4725 | } |
| 4726 | #endif |
| 4727 | |
| 4728 | if (hb->bDAnr) |
| 4729 | { |
| 4730 | int i, j, nleg; |
| 4731 | char **legnames; |
| 4732 | char buf[STRLEN4096]; |
| 4733 | |
| 4734 | #define USE_THIS_GROUP(j)( (j == gr0) || (bTwo && (j == gr1)) ) ( (j == gr0) || (bTwo && (j == gr1)) ) |
| 4735 | |
| 4736 | fp = xvgropen(opt2fn("-dan", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
| 4737 | "Donors and Acceptors", output_env_get_xvgr_tlabel(oenv), "Count", oenv); |
| 4738 | nleg = (bTwo ? 2 : 1)*2; |
| 4739 | snew(legnames, nleg)(legnames) = save_calloc("legnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4739, (nleg), sizeof(*(legnames))); |
| 4740 | i = 0; |
| 4741 | for (j = 0; j < grNR; j++) |
| 4742 | { |
| 4743 | if (USE_THIS_GROUP(j)( (j == gr0) || (bTwo && (j == gr1)) ) ) |
| 4744 | { |
| 4745 | sprintf(buf, "Donors %s", grpnames[j]); |
| 4746 | legnames[i++] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); |
| 4747 | sprintf(buf, "Acceptors %s", grpnames[j]); |
| 4748 | legnames[i++] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); |
| 4749 | } |
| 4750 | } |
| 4751 | if (i != nleg) |
| 4752 | { |
| 4753 | gmx_incons("number of legend entries")_gmx_error("incons", "number of legend entries", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_hbond.c" , 4753); |
| 4754 | } |
| 4755 | xvgr_legend(fp, nleg, (const char**)legnames, oenv); |
| 4756 | for (i = 0; i < nframes; i++) |
| 4757 | { |
| 4758 | fprintf(fp, "%10g", hb->time[i]); |
| 4759 | for (j = 0; (j < grNR); j++) |
| 4760 | { |
| 4761 | if (USE_THIS_GROUP(j)( (j == gr0) || (bTwo && (j == gr1)) ) ) |
| 4762 | { |
| 4763 | fprintf(fp, " %6d", hb->danr[i][j]); |
| 4764 | } |
| 4765 | } |
| 4766 | fprintf(fp, "\n"); |
| 4767 | } |
| 4768 | gmx_ffclose(fp); |
| 4769 | } |
| 4770 | |
| 4771 | return 0; |
| 4772 | } |