| File: | gromacs/gmxana/gmx_clustsize.c |
| Location: | line 148, column 18 |
| Description: | Access to field 'nr' results in a dereference of a null pointer (loaded from variable 'mols') |
| 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-2007, 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 | ||||
| 41 | #include <math.h> | |||
| 42 | ||||
| 43 | #include "typedefs.h" | |||
| 44 | #include "macros.h" | |||
| 45 | #include "gromacs/math/vec.h" | |||
| 46 | #include "pbc.h" | |||
| 47 | #include "rmpbc.h" | |||
| 48 | #include "gromacs/commandline/pargs.h" | |||
| 49 | #include "gromacs/fileio/xvgr.h" | |||
| 50 | #include "gromacs/utility/futil.h" | |||
| 51 | #include "gromacs/fileio/tpxio.h" | |||
| 52 | #include "gromacs/fileio/trxio.h" | |||
| 53 | #include "index.h" | |||
| 54 | #include "gromacs/utility/smalloc.h" | |||
| 55 | #include "calcgrid.h" | |||
| 56 | #include "nrnb.h" | |||
| 57 | #include "physics.h" | |||
| 58 | #include "coulomb.h" | |||
| 59 | #include "pme.h" | |||
| 60 | #include "gstat.h" | |||
| 61 | #include "gromacs/fileio/matio.h" | |||
| 62 | #include "mtop_util.h" | |||
| 63 | #include "gmx_ana.h" | |||
| 64 | ||||
| 65 | #include "gromacs/utility/fatalerror.h" | |||
| 66 | ||||
| 67 | static void clust_size(const char *ndx, const char *trx, const char *xpm, | |||
| 68 | const char *xpmw, const char *ncl, const char *acl, | |||
| 69 | const char *mcl, const char *histo, const char *tempf, | |||
| 70 | const char *mcn, gmx_bool bMol, gmx_bool bPBC, const char *tpr, | |||
| 71 | real cut, int nskip, int nlevels, | |||
| 72 | t_rgb rmid, t_rgb rhi, int ndf, | |||
| 73 | const output_env_t oenv) | |||
| 74 | { | |||
| 75 | FILE *fp, *gp, *hp, *tp; | |||
| 76 | atom_id *index = NULL((void*)0); | |||
| 77 | int nindex, natoms; | |||
| 78 | t_trxstatus *status; | |||
| 79 | rvec *x = NULL((void*)0), *v = NULL((void*)0), dx; | |||
| 80 | t_pbc pbc; | |||
| 81 | char *gname; | |||
| 82 | char timebuf[32]; | |||
| 83 | gmx_bool bSame, bTPRwarn = TRUE1; | |||
| 84 | /* Topology stuff */ | |||
| 85 | t_trxframe fr; | |||
| 86 | t_tpxheader tpxh; | |||
| 87 | gmx_mtop_t *mtop = NULL((void*)0); | |||
| 88 | int ePBC = -1; | |||
| 89 | t_block *mols = NULL((void*)0); | |||
| 90 | gmx_mtop_atomlookup_t alook; | |||
| 91 | t_atom *atom; | |||
| 92 | int version, generation, ii, jj, nsame; | |||
| 93 | real temp, tfac; | |||
| 94 | /* Cluster size distribution (matrix) */ | |||
| 95 | real **cs_dist = NULL((void*)0); | |||
| 96 | real tf, dx2, cut2, *t_x = NULL((void*)0), *t_y, cmid, cmax, cav, ekin; | |||
| 97 | int i, j, k, ai, aj, ak, ci, cj, nframe, nclust, n_x, n_y, max_size = 0; | |||
| 98 | int *clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto; | |||
| 99 | t_rgb rlo = { 1.0, 1.0, 1.0 }; | |||
| 100 | ||||
| 101 | clear_trxframe(&fr, TRUE1); | |||
| 102 | sprintf(timebuf, "Time (%s)", output_env_get_time_unit(oenv)); | |||
| 103 | tf = output_env_get_time_factor(oenv); | |||
| 104 | fp = xvgropen(ncl, "Number of clusters", timebuf, "N", oenv); | |||
| 105 | gp = xvgropen(acl, "Average cluster size", timebuf, "#molecules", oenv); | |||
| 106 | hp = xvgropen(mcl, "Max cluster size", timebuf, "#molecules", oenv); | |||
| 107 | tp = xvgropen(tempf, "Temperature of largest cluster", timebuf, "T (K)", | |||
| 108 | oenv); | |||
| 109 | ||||
| 110 | if (!read_first_frame(oenv, &status, trx, &fr, TRX_NEED_X(1<<1) | TRX_READ_V(1<<2))) | |||
| ||||
| 111 | { | |||
| 112 | gmx_file(trx)_gmx_error("file", trx, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 112); | |||
| 113 | } | |||
| 114 | ||||
| 115 | natoms = fr.natoms; | |||
| 116 | x = fr.x; | |||
| 117 | ||||
| 118 | if (tpr) | |||
| 119 | { | |||
| 120 | snew(mtop, 1)(mtop) = save_calloc("mtop", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 120, (1), sizeof(*(mtop))); | |||
| 121 | read_tpxheader(tpr, &tpxh, TRUE1, &version, &generation); | |||
| 122 | if (tpxh.natoms != natoms) | |||
| 123 | { | |||
| 124 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 124, "tpr (%d atoms) and trajectory (%d atoms) do not match!", | |||
| 125 | tpxh.natoms, natoms); | |||
| 126 | } | |||
| 127 | ePBC = read_tpx(tpr, NULL((void*)0), NULL((void*)0), &natoms, NULL((void*)0), NULL((void*)0), NULL((void*)0), mtop); | |||
| 128 | } | |||
| 129 | if (ndf <= -1) | |||
| 130 | { | |||
| 131 | tfac = 1; | |||
| 132 | } | |||
| 133 | else | |||
| 134 | { | |||
| 135 | tfac = ndf/(3.0*natoms); | |||
| 136 | } | |||
| 137 | ||||
| 138 | if (bMol) | |||
| 139 | { | |||
| 140 | if (ndx) | |||
| 141 | { | |||
| 142 | printf("Using molecules rather than atoms. Not reading index file %s\n", | |||
| 143 | ndx); | |||
| 144 | } | |||
| 145 | mols = &(mtop->mols); | |||
| 146 | ||||
| 147 | /* Make dummy index */ | |||
| 148 | nindex = mols->nr; | |||
| ||||
| 149 | snew(index, nindex)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 149, (nindex), sizeof(*(index))); | |||
| 150 | for (i = 0; (i < nindex); i++) | |||
| 151 | { | |||
| 152 | index[i] = i; | |||
| 153 | } | |||
| 154 | gname = strdup("mols")(__extension__ (__builtin_constant_p ("mols") && ((size_t )(const void *)(("mols") + 1) - (size_t)(const void *)("mols" ) == 1) ? (((const char *) ("mols"))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len = strlen ("mols") + 1; char *__retval = (char *) malloc (__len); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval, "mols", __len ); __retval; })) : __strdup ("mols"))); | |||
| 155 | } | |||
| 156 | else | |||
| 157 | { | |||
| 158 | rd_index(ndx, 1, &nindex, &index, &gname); | |||
| 159 | } | |||
| 160 | ||||
| 161 | alook = gmx_mtop_atomlookup_init(mtop); | |||
| 162 | ||||
| 163 | snew(clust_index, nindex)(clust_index) = save_calloc("clust_index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 163, (nindex), sizeof(*(clust_index))); | |||
| 164 | snew(clust_size, nindex)(clust_size) = save_calloc("clust_size", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 164, (nindex), sizeof(*(clust_size))); | |||
| 165 | cut2 = cut*cut; | |||
| 166 | nframe = 0; | |||
| 167 | n_x = 0; | |||
| 168 | snew(t_y, nindex)(t_y) = save_calloc("t_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 168, (nindex), sizeof(*(t_y))); | |||
| 169 | for (i = 0; (i < nindex); i++) | |||
| 170 | { | |||
| 171 | t_y[i] = i+1; | |||
| 172 | } | |||
| 173 | max_clust_size = 1; | |||
| 174 | max_clust_ind = -1; | |||
| 175 | do | |||
| 176 | { | |||
| 177 | if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0))) | |||
| 178 | { | |||
| 179 | if (bPBC) | |||
| 180 | { | |||
| 181 | set_pbc(&pbc, ePBC, fr.box); | |||
| 182 | } | |||
| 183 | max_clust_size = 1; | |||
| 184 | max_clust_ind = -1; | |||
| 185 | ||||
| 186 | /* Put all atoms/molecules in their own cluster, with size 1 */ | |||
| 187 | for (i = 0; (i < nindex); i++) | |||
| 188 | { | |||
| 189 | /* Cluster index is indexed with atom index number */ | |||
| 190 | clust_index[i] = i; | |||
| 191 | /* Cluster size is indexed with cluster number */ | |||
| 192 | clust_size[i] = 1; | |||
| 193 | } | |||
| 194 | ||||
| 195 | /* Loop over atoms */ | |||
| 196 | for (i = 0; (i < nindex); i++) | |||
| 197 | { | |||
| 198 | ai = index[i]; | |||
| 199 | ci = clust_index[i]; | |||
| 200 | ||||
| 201 | /* Loop over atoms (only half a matrix) */ | |||
| 202 | for (j = i+1; (j < nindex); j++) | |||
| 203 | { | |||
| 204 | cj = clust_index[j]; | |||
| 205 | ||||
| 206 | /* If they are not in the same cluster already */ | |||
| 207 | if (ci != cj) | |||
| 208 | { | |||
| 209 | aj = index[j]; | |||
| 210 | ||||
| 211 | /* Compute distance */ | |||
| 212 | if (bMol) | |||
| 213 | { | |||
| 214 | bSame = FALSE0; | |||
| 215 | for (ii = mols->index[ai]; !bSame && (ii < mols->index[ai+1]); ii++) | |||
| 216 | { | |||
| 217 | for (jj = mols->index[aj]; !bSame && (jj < mols->index[aj+1]); jj++) | |||
| 218 | { | |||
| 219 | if (bPBC) | |||
| 220 | { | |||
| 221 | pbc_dx(&pbc, x[ii], x[jj], dx); | |||
| 222 | } | |||
| 223 | else | |||
| 224 | { | |||
| 225 | rvec_sub(x[ii], x[jj], dx); | |||
| 226 | } | |||
| 227 | dx2 = iprod(dx, dx); | |||
| 228 | bSame = (dx2 < cut2); | |||
| 229 | } | |||
| 230 | } | |||
| 231 | } | |||
| 232 | else | |||
| 233 | { | |||
| 234 | if (bPBC) | |||
| 235 | { | |||
| 236 | pbc_dx(&pbc, x[ai], x[aj], dx); | |||
| 237 | } | |||
| 238 | else | |||
| 239 | { | |||
| 240 | rvec_sub(x[ai], x[aj], dx); | |||
| 241 | } | |||
| 242 | dx2 = iprod(dx, dx); | |||
| 243 | bSame = (dx2 < cut2); | |||
| 244 | } | |||
| 245 | /* If distance less than cut-off */ | |||
| 246 | if (bSame) | |||
| 247 | { | |||
| 248 | /* Merge clusters: check for all atoms whether they are in | |||
| 249 | * cluster cj and if so, put them in ci | |||
| 250 | */ | |||
| 251 | for (k = 0; (k < nindex); k++) | |||
| 252 | { | |||
| 253 | if (clust_index[k] == cj) | |||
| 254 | { | |||
| 255 | if (clust_size[cj] <= 0) | |||
| 256 | { | |||
| 257 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 257, "negative cluster size %d for element %d", | |||
| 258 | clust_size[cj], cj); | |||
| 259 | } | |||
| 260 | clust_size[cj]--; | |||
| 261 | clust_index[k] = ci; | |||
| 262 | clust_size[ci]++; | |||
| 263 | } | |||
| 264 | } | |||
| 265 | } | |||
| 266 | } | |||
| 267 | } | |||
| 268 | } | |||
| 269 | n_x++; | |||
| 270 | srenew(t_x, n_x)(t_x) = save_realloc("t_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 270, (t_x), (n_x), sizeof(*(t_x))); | |||
| 271 | t_x[n_x-1] = fr.time*tf; | |||
| 272 | srenew(cs_dist, n_x)(cs_dist) = save_realloc("cs_dist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 272, (cs_dist), (n_x), sizeof(*(cs_dist))); | |||
| 273 | snew(cs_dist[n_x-1], nindex)(cs_dist[n_x-1]) = save_calloc("cs_dist[n_x-1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 273, (nindex), sizeof(*(cs_dist[n_x-1]))); | |||
| 274 | nclust = 0; | |||
| 275 | cav = 0; | |||
| 276 | nav = 0; | |||
| 277 | for (i = 0; (i < nindex); i++) | |||
| 278 | { | |||
| 279 | ci = clust_size[i]; | |||
| 280 | if (ci > max_clust_size) | |||
| 281 | { | |||
| 282 | max_clust_size = ci; | |||
| 283 | max_clust_ind = i; | |||
| 284 | } | |||
| 285 | if (ci > 0) | |||
| 286 | { | |||
| 287 | nclust++; | |||
| 288 | cs_dist[n_x-1][ci-1] += 1.0; | |||
| 289 | max_size = max(max_size, ci)(((max_size) > (ci)) ? (max_size) : (ci) ); | |||
| 290 | if (ci > 1) | |||
| 291 | { | |||
| 292 | cav += ci; | |||
| 293 | nav++; | |||
| 294 | } | |||
| 295 | } | |||
| 296 | } | |||
| 297 | fprintf(fp, "%14.6e %10d\n", fr.time, nclust); | |||
| 298 | if (nav > 0) | |||
| 299 | { | |||
| 300 | fprintf(gp, "%14.6e %10.3f\n", fr.time, cav/nav); | |||
| 301 | } | |||
| 302 | fprintf(hp, "%14.6e %10d\n", fr.time, max_clust_size); | |||
| 303 | } | |||
| 304 | /* Analyse velocities, if present */ | |||
| 305 | if (fr.bV) | |||
| 306 | { | |||
| 307 | if (!tpr) | |||
| 308 | { | |||
| 309 | if (bTPRwarn) | |||
| 310 | { | |||
| 311 | printf("You need a [TT].tpr[tt] file to analyse temperatures\n"); | |||
| 312 | bTPRwarn = FALSE0; | |||
| 313 | } | |||
| 314 | } | |||
| 315 | else | |||
| 316 | { | |||
| 317 | v = fr.v; | |||
| 318 | /* Loop over clusters and for each cluster compute 1/2 m v^2 */ | |||
| 319 | if (max_clust_ind >= 0) | |||
| 320 | { | |||
| 321 | ekin = 0; | |||
| 322 | for (i = 0; (i < nindex); i++) | |||
| 323 | { | |||
| 324 | if (clust_index[i] == max_clust_ind) | |||
| 325 | { | |||
| 326 | ai = index[i]; | |||
| 327 | gmx_mtop_atomnr_to_atom(alook, ai, &atom); | |||
| 328 | ekin += 0.5*atom->m*iprod(v[ai], v[ai]); | |||
| 329 | } | |||
| 330 | } | |||
| 331 | temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))); | |||
| 332 | fprintf(tp, "%10.3f %10.3f\n", fr.time, temp); | |||
| 333 | } | |||
| 334 | } | |||
| 335 | } | |||
| 336 | nframe++; | |||
| 337 | } | |||
| 338 | while (read_next_frame(oenv, status, &fr)); | |||
| 339 | close_trx(status); | |||
| 340 | gmx_ffclose(fp); | |||
| 341 | gmx_ffclose(gp); | |||
| 342 | gmx_ffclose(hp); | |||
| 343 | gmx_ffclose(tp); | |||
| 344 | ||||
| 345 | gmx_mtop_atomlookup_destroy(alook); | |||
| 346 | ||||
| 347 | if (max_clust_ind >= 0) | |||
| 348 | { | |||
| 349 | fp = gmx_ffopen(mcn, "w"); | |||
| 350 | fprintf(fp, "[ max_clust ]\n"); | |||
| 351 | for (i = 0; (i < nindex); i++) | |||
| 352 | { | |||
| 353 | if (clust_index[i] == max_clust_ind) | |||
| 354 | { | |||
| 355 | if (bMol) | |||
| 356 | { | |||
| 357 | for (j = mols->index[i]; (j < mols->index[i+1]); j++) | |||
| 358 | { | |||
| 359 | fprintf(fp, "%d\n", j+1); | |||
| 360 | } | |||
| 361 | } | |||
| 362 | else | |||
| 363 | { | |||
| 364 | fprintf(fp, "%d\n", index[i]+1); | |||
| 365 | } | |||
| 366 | } | |||
| 367 | } | |||
| 368 | gmx_ffclose(fp); | |||
| 369 | } | |||
| 370 | ||||
| 371 | /* Print the real distribution cluster-size/numer, averaged over the trajectory. */ | |||
| 372 | fp = xvgropen(histo, "Cluster size distribution", "Cluster size", "()", oenv); | |||
| 373 | nhisto = 0; | |||
| 374 | fprintf(fp, "%5d %8.3f\n", 0, 0.0); | |||
| 375 | for (j = 0; (j < max_size); j++) | |||
| 376 | { | |||
| 377 | real nelem = 0; | |||
| 378 | for (i = 0; (i < n_x); i++) | |||
| 379 | { | |||
| 380 | nelem += cs_dist[i][j]; | |||
| 381 | } | |||
| 382 | fprintf(fp, "%5d %8.3f\n", j+1, nelem/n_x); | |||
| 383 | nhisto += (int)((j+1)*nelem/n_x); | |||
| 384 | } | |||
| 385 | fprintf(fp, "%5d %8.3f\n", j+1, 0.0); | |||
| 386 | gmx_ffclose(fp); | |||
| 387 | ||||
| 388 | fprintf(stderrstderr, "Total number of atoms in clusters = %d\n", nhisto); | |||
| 389 | ||||
| 390 | /* Look for the smallest entry that is not zero | |||
| 391 | * This will make that zero is white, and not zero is coloured. | |||
| 392 | */ | |||
| 393 | cmid = 100.0; | |||
| 394 | cmax = 0.0; | |||
| 395 | for (i = 0; (i < n_x); i++) | |||
| 396 | { | |||
| 397 | for (j = 0; (j < max_size); j++) | |||
| 398 | { | |||
| 399 | if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid)) | |||
| 400 | { | |||
| 401 | cmid = cs_dist[i][j]; | |||
| 402 | } | |||
| 403 | cmax = max(cs_dist[i][j], cmax)(((cs_dist[i][j]) > (cmax)) ? (cs_dist[i][j]) : (cmax) ); | |||
| 404 | } | |||
| 405 | } | |||
| 406 | fprintf(stderrstderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size); | |||
| 407 | cmid = 1; | |||
| 408 | fp = gmx_ffopen(xpm, "w"); | |||
| 409 | write_xpm3(fp, 0, "Cluster size distribution", "# clusters", timebuf, "Size", | |||
| 410 | n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax, | |||
| 411 | rlo, rmid, rhi, &nlevels); | |||
| 412 | gmx_ffclose(fp); | |||
| 413 | cmid = 100.0; | |||
| 414 | cmax = 0.0; | |||
| 415 | for (i = 0; (i < n_x); i++) | |||
| 416 | { | |||
| 417 | for (j = 0; (j < max_size); j++) | |||
| 418 | { | |||
| 419 | cs_dist[i][j] *= (j+1); | |||
| 420 | if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid)) | |||
| 421 | { | |||
| 422 | cmid = cs_dist[i][j]; | |||
| 423 | } | |||
| 424 | cmax = max(cs_dist[i][j], cmax)(((cs_dist[i][j]) > (cmax)) ? (cs_dist[i][j]) : (cmax) ); | |||
| 425 | } | |||
| 426 | } | |||
| 427 | fprintf(stderrstderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size); | |||
| 428 | fp = gmx_ffopen(xpmw, "w"); | |||
| 429 | write_xpm3(fp, 0, "Weighted cluster size distribution", "Fraction", timebuf, | |||
| 430 | "Size", n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax, | |||
| 431 | rlo, rmid, rhi, &nlevels); | |||
| 432 | gmx_ffclose(fp); | |||
| 433 | ||||
| 434 | sfree(clust_index)save_free("clust_index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 434, (clust_index)); | |||
| 435 | sfree(clust_size)save_free("clust_size", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 435, (clust_size)); | |||
| 436 | sfree(index)save_free("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 436, (index)); | |||
| 437 | } | |||
| 438 | ||||
| 439 | int gmx_clustsize(int argc, char *argv[]) | |||
| 440 | { | |||
| 441 | const char *desc[] = { | |||
| 442 | "[THISMODULE] computes the size distributions of molecular/atomic clusters in", | |||
| 443 | "the gas phase. The output is given in the form of an [TT].xpm[tt] file.", | |||
| 444 | "The total number of clusters is written to an [TT].xvg[tt] file.[PAR]", | |||
| 445 | "When the [TT]-mol[tt] option is given clusters will be made out of", | |||
| 446 | "molecules rather than atoms, which allows clustering of large molecules.", | |||
| 447 | "In this case an index file would still contain atom numbers", | |||
| 448 | "or your calculation will die with a SEGV.[PAR]", | |||
| 449 | "When velocities are present in your trajectory, the temperature of", | |||
| 450 | "the largest cluster will be printed in a separate [TT].xvg[tt] file assuming", | |||
| 451 | "that the particles are free to move. If you are using constraints,", | |||
| 452 | "please correct the temperature. For instance water simulated with SHAKE", | |||
| 453 | "or SETTLE will yield a temperature that is 1.5 times too low. You can", | |||
| 454 | "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal", | |||
| 455 | "of center of mass motion into account.[PAR]", | |||
| 456 | "The [TT]-mc[tt] option will produce an index file containing the", | |||
| 457 | "atom numbers of the largest cluster." | |||
| 458 | }; | |||
| 459 | ||||
| 460 | static real cutoff = 0.35; | |||
| 461 | static int nskip = 0; | |||
| 462 | static int nlevels = 20; | |||
| 463 | static int ndf = -1; | |||
| 464 | static gmx_bool bMol = FALSE0; | |||
| 465 | static gmx_bool bPBC = TRUE1; | |||
| 466 | static rvec rlo = { 1.0, 1.0, 0.0 }; | |||
| 467 | static rvec rhi = { 0.0, 0.0, 1.0 }; | |||
| 468 | ||||
| 469 | output_env_t oenv; | |||
| 470 | ||||
| 471 | t_pargs pa[] = { | |||
| 472 | { "-cut", FALSE0, etREAL, {&cutoff}, | |||
| 473 | "Largest distance (nm) to be considered in a cluster" }, | |||
| 474 | { "-mol", FALSE0, etBOOL, {&bMol}, | |||
| 475 | "Cluster molecules rather than atoms (needs [TT].tpr[tt] file)" }, | |||
| 476 | { "-pbc", FALSE0, etBOOL, {&bPBC}, | |||
| 477 | "Use periodic boundary conditions" }, | |||
| 478 | { "-nskip", FALSE0, etINT, {&nskip}, | |||
| 479 | "Number of frames to skip between writing" }, | |||
| 480 | { "-nlevels", FALSE0, etINT, {&nlevels}, | |||
| 481 | "Number of levels of grey in [TT].xpm[tt] output" }, | |||
| 482 | { "-ndf", FALSE0, etINT, {&ndf}, | |||
| 483 | "Number of degrees of freedom of the entire system for temperature calculation. If not set, the number of atoms times three is used." }, | |||
| 484 | { "-rgblo", FALSE0, etRVEC, {rlo}, | |||
| 485 | "RGB values for the color of the lowest occupied cluster size" }, | |||
| 486 | { "-rgbhi", FALSE0, etRVEC, {rhi}, | |||
| 487 | "RGB values for the color of the highest occupied cluster size" } | |||
| 488 | }; | |||
| 489 | #define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))) | |||
| 490 | const char *fnNDX, *fnTPR; | |||
| 491 | gmx_bool bSQ, bRDF; | |||
| 492 | t_rgb rgblo, rgbhi; | |||
| 493 | ||||
| 494 | t_filenm fnm[] = { | |||
| 495 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, | |||
| 496 | { efTPR, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
| 497 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
| 498 | { efXPM, "-o", "csize", ffWRITE1<<2 }, | |||
| 499 | { efXPM, "-ow", "csizew", ffWRITE1<<2 }, | |||
| 500 | { efXVG, "-nc", "nclust", ffWRITE1<<2 }, | |||
| 501 | { efXVG, "-mc", "maxclust", ffWRITE1<<2 }, | |||
| 502 | { efXVG, "-ac", "avclust", ffWRITE1<<2 }, | |||
| 503 | { efXVG, "-hc", "histo-clust", ffWRITE1<<2 }, | |||
| 504 | { efXVG, "-temp", "temp", ffOPTWR(1<<2| 1<<3) }, | |||
| 505 | { efNDX, "-mcn", "maxclust", ffOPTWR(1<<2| 1<<3) } | |||
| 506 | }; | |||
| 507 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) | |||
| 508 | ||||
| 509 | if (!parse_common_args(&argc, argv, | |||
| 510 | PCA_CAN_VIEW(1<<5) | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13), | |||
| 511 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) | |||
| 512 | { | |||
| 513 | return 0; | |||
| 514 | } | |||
| 515 | ||||
| 516 | fnNDX = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 517 | rgblo.r = rlo[XX0], rgblo.g = rlo[YY1], rgblo.b = rlo[ZZ2]; | |||
| 518 | rgbhi.r = rhi[XX0], rgbhi.g = rhi[YY1], rgbhi.b = rhi[ZZ2]; | |||
| 519 | ||||
| 520 | fnTPR = ftp2fn_null(efTPR, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); | |||
| 521 | if (bMol && !fnTPR) | |||
| 522 | { | |||
| 523 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_clustsize.c" , 523, "You need a tpr file for the -mol option"); | |||
| 524 | } | |||
| 525 | ||||
| 526 | clust_size(fnNDX, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
| 527 | opt2fn("-ow", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
| 528 | opt2fn("-nc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-ac", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
| 529 | opt2fn("-mc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-hc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
| 530 | opt2fn("-temp", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-mcn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
| 531 | bMol, bPBC, fnTPR, | |||
| 532 | cutoff, nskip, nlevels, rgblo, rgbhi, ndf, oenv); | |||
| 533 | ||||
| 534 | return 0; | |||
| 535 | } |