?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Add replacements for pbc enumerations
[gromacs.git] / src / gromacs / gmxana / gmx_clustsize.cpp
blob7dbf95759ff96495474d3a4ae435a5dbcd4aa0e9
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,2015,2016,2017,2018,2019, 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 #include "gmxpre.h"
38
39 #include <cmath>
40
41 #include <algorithm>
42
43 #include "gromacs/commandline/filenm.h"
44 #include "gromacs/commandline/pargs.h"
45 #include "gromacs/fileio/matio.h"
46 #include "gromacs/fileio/tpxio.h"
47 #include "gromacs/fileio/trxio.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/gmxana/gmx_ana.h"
50 #include "gromacs/gmxana/gstat.h"
51 #include "gromacs/gmxlib/nrnb.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/pbcutil/pbc.h"
55 #include "gromacs/topology/index.h"
56 #include "gromacs/topology/mtop_lookup.h"
57 #include "gromacs/topology/mtop_util.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/trajectory/trajectoryframe.h"
60 #include "gromacs/utility/arraysize.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/gmxassert.h"
65 #include "gromacs/utility/smalloc.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 gmx_output_env_t *oenv)
74 {
75 FILE *fp, *gp, *hp, *tp;
76 int *index = nullptr;
77 int nindex, natoms;
78 t_trxstatus *status;
79 rvec *x = nullptr, *v = nullptr, dx;
80 t_pbc pbc;
81 gmx_bool bSame, bTPRwarn = TRUE;
82 /* Topology stuff */
83 t_trxframe fr;
84 TpxFileHeader tpxh;
85 gmx_mtop_t *mtop = nullptr;
86 int ePBC = -1;
87 int ii, jj;
88 real temp, tfac;
89 /* Cluster size distribution (matrix) */
90 real **cs_dist = nullptr;
91 real tf, dx2, cut2, *t_x = nullptr, *t_y, cmid, cmax, cav, ekin;
92 int i, j, k, ai, aj, ci, cj, nframe, nclust, n_x, max_size = 0;
93 int *clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto;
94 t_rgb rlo = { 1.0, 1.0, 1.0 };
95 int frameCounter = 0;
96 real frameTime;
97
98 clear_trxframe(&fr, TRUE);
99 auto timeLabel = output_env_get_time_label(oenv);
100 tf = output_env_get_time_factor(oenv);
101 fp = xvgropen(ncl, "Number of clusters", timeLabel, "N", oenv);
102 gp = xvgropen(acl, "Average cluster size", timeLabel, "#molecules", oenv);
103 hp = xvgropen(mcl, "Max cluster size", timeLabel, "#molecules", oenv);
104 tp = xvgropen(tempf, "Temperature of largest cluster", timeLabel, "T (K)",
105 oenv);
106
107 if (!read_first_frame(oenv, &status, trx, &fr, TRX_NEED_X | TRX_READ_V))
108 {
109 gmx_file(trx);
110 }
111
112 natoms = fr.natoms;
113 x = fr.x;
114
115 if (tpr)
116 {
117 mtop = new gmx_mtop_t;
118 tpxh = readTpxHeader(tpr, true);
119 if (tpxh.natoms != natoms)
120 {
121 gmx_fatal(FARGS, "tpr (%d atoms) and trajectory (%d atoms) do not match!",
122 tpxh.natoms, natoms);
123 }
124 ePBC = read_tpx(tpr, nullptr, nullptr, &natoms, nullptr, nullptr, mtop);
125 }
126 if (ndf <= -1)
127 {
128 tfac = 1;
129 }
130 else
131 {
132 tfac = ndf/(3.0*natoms);
133 }
134
135 gmx::RangePartitioning mols;
136 if (bMol)
137 {
138 if (ndx)
139 {
140 printf("Using molecules rather than atoms. Not reading index file %s\n",
141 ndx);
142 }
143 GMX_RELEASE_ASSERT(mtop != nullptr, "Trying to access mtop->mols from NULL mtop pointer");
144 mols = gmx_mtop_molecules(*mtop);
145
146 /* Make dummy index */
147 nindex = mols.numBlocks();
148 snew(index, nindex);
149 for (i = 0; (i < nindex); i++)
150 {
151 index[i] = i;
152 }
153 }
154 else
155 {
156 char *gname;
157 rd_index(ndx, 1, &nindex, &index, &gname);
158 sfree(gname);
159 }
160
161 snew(clust_index, nindex);
162 snew(clust_size, nindex);
163 cut2 = cut*cut;
164 nframe = 0;
165 n_x = 0;
166 snew(t_y, nindex);
167 for (i = 0; (i < nindex); i++)
168 {
169 t_y[i] = i+1;
170 }
171 max_clust_size = 1;
172 max_clust_ind = -1;
173 int molb = 0;
174 do
175 {
176 if ((nskip == 0) || ((nskip > 0) && ((nframe % nskip) == 0)))
177 {
178 if (bPBC)
179 {
180 set_pbc(&pbc, ePBC, fr.box);
181 }
182 max_clust_size = 1;
183 max_clust_ind = -1;
184
185 /* Put all atoms/molecules in their own cluster, with size 1 */
186 for (i = 0; (i < nindex); i++)
187 {
188 /* Cluster index is indexed with atom index number */
189 clust_index[i] = i;
190 /* Cluster size is indexed with cluster number */
191 clust_size[i] = 1;
192 }
193
194 /* Loop over atoms */
195 for (i = 0; (i < nindex); i++)
196 {
197 ai = index[i];
198 ci = clust_index[i];
199
200 /* Loop over atoms (only half a matrix) */
201 for (j = i+1; (j < nindex); j++)
202 {
203 cj = clust_index[j];
204
205 /* If they are not in the same cluster already */
206 if (ci != cj)
207 {
208 aj = index[j];
209
210 /* Compute distance */
211 if (bMol)
212 {
213 GMX_RELEASE_ASSERT(mols.numBlocks() > 0, "Cannot access index[] from empty mols");
214 bSame = FALSE;
215 for (ii = mols.block(ai).begin(); !bSame && ii < mols.block(ai).end(); ii++)
216 {
217 for (jj = mols.block(aj).begin(); !bSame && jj < mols.block(aj).end(); 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(FARGS, "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);
271 if (fr.bTime)
272 {
273 frameTime = fr.time;
274 }
275 else if (fr.bStep)
276 {
277 frameTime = fr.step;
278 }
279 else
280 {
281 frameTime = ++frameCounter;
282 }
283 t_x[n_x-1] = frameTime*tf;
284 srenew(cs_dist, n_x);
285 snew(cs_dist[n_x-1], nindex);
286 nclust = 0;
287 cav = 0;
288 nav = 0;
289 for (i = 0; (i < nindex); i++)
290 {
291 ci = clust_size[i];
292 if (ci > max_clust_size)
293 {
294 max_clust_size = ci;
295 max_clust_ind = i;
296 }
297 if (ci > 0)
298 {
299 nclust++;
300 cs_dist[n_x-1][ci-1] += 1.0;
301 max_size = std::max(max_size, ci);
302 if (ci > 1)
303 {
304 cav += ci;
305 nav++;
306 }
307 }
308 }
309 fprintf(fp, "%14.6e %10d\n", frameTime, nclust);
310 if (nav > 0)
311 {
312 fprintf(gp, "%14.6e %10.3f\n", frameTime, cav/nav);
313 }
314 fprintf(hp, "%14.6e %10d\n", frameTime, max_clust_size);
315 }
316 /* Analyse velocities, if present */
317 if (fr.bV)
318 {
319 if (!tpr)
320 {
321 if (bTPRwarn)
322 {
323 printf("You need a [REF].tpr[ref] file to analyse temperatures\n");
324 bTPRwarn = FALSE;
325 }
326 }
327 else
328 {
329 v = fr.v;
330 /* Loop over clusters and for each cluster compute 1/2 m v^2 */
331 if (max_clust_ind >= 0)
332 {
333 ekin = 0;
334 for (i = 0; (i < nindex); i++)
335 {
336 if (clust_index[i] == max_clust_ind)
337 {
338 ai = index[i];
339 real m = mtopGetAtomMass(mtop, ai, &molb);
340 ekin += 0.5*m*iprod(v[ai], v[ai]);
341 }
342 }
343 temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
344 fprintf(tp, "%10.3f %10.3f\n", frameTime, temp);
345 }
346 }
347 }
348 nframe++;
349 }
350 while (read_next_frame(oenv, status, &fr));
351 close_trx(status);
352 done_frame(&fr);
353 xvgrclose(fp);
354 xvgrclose(gp);
355 xvgrclose(hp);
356 xvgrclose(tp);
357
358 if (max_clust_ind >= 0)
359 {
360 fp = gmx_ffopen(mcn, "w");
361 fprintf(fp, "[ max_clust ]\n");
362 for (i = 0; (i < nindex); i++)
363 {
364 if (clust_index[i] == max_clust_ind)
365 {
366 if (bMol)
367 {
368 GMX_RELEASE_ASSERT(mols.numBlocks() > 0, "Cannot access index[] from empty mols");
369 for (int j : mols.block(i))
370 {
371 fprintf(fp, "%d\n", j+1);
372 }
373 }
374 else
375 {
376 fprintf(fp, "%d\n", index[i]+1);
377 }
378 }
379 }
380 gmx_ffclose(fp);
381 }
382
383 /* Print the real distribution cluster-size/numer, averaged over the trajectory. */
384 fp = xvgropen(histo, "Cluster size distribution", "Cluster size", "()", oenv);
385 nhisto = 0;
386 fprintf(fp, "%5d %8.3f\n", 0, 0.0);
387 for (j = 0; (j < max_size); j++)
388 {
389 real nelem = 0;
390 for (i = 0; (i < n_x); i++)
391 {
392 nelem += cs_dist[i][j];
393 }
394 fprintf(fp, "%5d %8.3f\n", j+1, nelem/n_x);
395 nhisto += static_cast<int>((j+1)*nelem/n_x);
396 }
397 fprintf(fp, "%5d %8.3f\n", j+1, 0.0);
398 xvgrclose(fp);
399
400 fprintf(stderr, "Total number of atoms in clusters = %d\n", nhisto);
401
402 /* Look for the smallest entry that is not zero
403 * This will make that zero is white, and not zero is coloured.
404 */
405 cmid = 100.0;
406 cmax = 0.0;
407 for (i = 0; (i < n_x); i++)
408 {
409 for (j = 0; (j < max_size); j++)
410 {
411 if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
412 {
413 cmid = cs_dist[i][j];
414 }
415 cmax = std::max(cs_dist[i][j], cmax);
416 }
417 }
418 fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
419 cmid = 1;
420 fp = gmx_ffopen(xpm, "w");
421 write_xpm3(fp, 0, "Cluster size distribution", "# clusters", timeLabel, "Size",
422 n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax,
423 rlo, rmid, rhi, &nlevels);
424 gmx_ffclose(fp);
425 cmid = 100.0;
426 cmax = 0.0;
427 for (i = 0; (i < n_x); i++)
428 {
429 for (j = 0; (j < max_size); j++)
430 {
431 cs_dist[i][j] *= (j+1);
432 if ((cs_dist[i][j] > 0) && (cs_dist[i][j] < cmid))
433 {
434 cmid = cs_dist[i][j];
435 }
436 cmax = std::max(cs_dist[i][j], cmax);
437 }
438 }
439 fprintf(stderr, "cmid: %g, cmax: %g, max_size: %d\n", cmid, cmax, max_size);
440 fp = gmx_ffopen(xpmw, "w");
441 write_xpm3(fp, 0, "Weighted cluster size distribution", "Fraction", timeLabel,
442 "Size", n_x, max_size, t_x, t_y, cs_dist, 0, cmid, cmax,
443 rlo, rmid, rhi, &nlevels);
444 gmx_ffclose(fp);
445 delete mtop;
446 sfree(t_x);
447 sfree(t_y);
448 for (i = 0; (i < n_x); i++)
449 {
450 sfree(cs_dist[i]);
451 }
452 sfree(cs_dist);
453 sfree(clust_index);
454 sfree(clust_size);
455 sfree(index);
456 }
457
458 int gmx_clustsize(int argc, char *argv[])
459 {
460 const char *desc[] = {
461 "[THISMODULE] computes the size distributions of molecular/atomic clusters in",
462 "the gas phase. The output is given in the form of an [REF].xpm[ref] file.",
463 "The total number of clusters is written to an [REF].xvg[ref] file.[PAR]",
464 "When the [TT]-mol[tt] option is given clusters will be made out of",
465 "molecules rather than atoms, which allows clustering of large molecules.",
466 "In this case an index file would still contain atom numbers",
467 "or your calculation will die with a SEGV.[PAR]",
468 "When velocities are present in your trajectory, the temperature of",
469 "the largest cluster will be printed in a separate [REF].xvg[ref] file assuming",
470 "that the particles are free to move. If you are using constraints,",
471 "please correct the temperature. For instance water simulated with SHAKE",
472 "or SETTLE will yield a temperature that is 1.5 times too low. You can",
473 "compensate for this with the [TT]-ndf[tt] option. Remember to take the removal",
474 "of center of mass motion into account.[PAR]",
475 "The [TT]-mc[tt] option will produce an index file containing the",
476 "atom numbers of the largest cluster."
477 };
478
479 real cutoff = 0.35;
480 int nskip = 0;
481 int nlevels = 20;
482 int ndf = -1;
483 gmx_bool bMol = FALSE;
484 gmx_bool bPBC = TRUE;
485 rvec rlo = { 1.0, 1.0, 0.0 };
486 rvec rhi = { 0.0, 0.0, 1.0 };
487
488 gmx_output_env_t *oenv;
489
490 t_pargs pa[] = {
491 { "-cut", FALSE, etREAL, {&cutoff},
492 "Largest distance (nm) to be considered in a cluster" },
493 { "-mol", FALSE, etBOOL, {&bMol},
494 "Cluster molecules rather than atoms (needs [REF].tpr[ref] file)" },
495 { "-pbc", FALSE, etBOOL, {&bPBC},
496 "Use periodic boundary conditions" },
497 { "-nskip", FALSE, etINT, {&nskip},
498 "Number of frames to skip between writing" },
499 { "-nlevels", FALSE, etINT, {&nlevels},
500 "Number of levels of grey in [REF].xpm[ref] output" },
501 { "-ndf", FALSE, etINT, {&ndf},
502 "Number of degrees of freedom of the entire system for temperature calculation. If not set, the number of atoms times three is used." },
503 { "-rgblo", FALSE, etRVEC, {rlo},
504 "RGB values for the color of the lowest occupied cluster size" },
505 { "-rgbhi", FALSE, etRVEC, {rhi},
506 "RGB values for the color of the highest occupied cluster size" }
507 };
508 #define NPA asize(pa)
509 const char *fnNDX, *fnTPR;
510 t_rgb rgblo, rgbhi;
511
512 t_filenm fnm[] = {
513 { efTRX, "-f", nullptr, ffREAD },
514 { efTPR, nullptr, nullptr, ffOPTRD },
515 { efNDX, nullptr, nullptr, ffOPTRD },
516 { efXPM, "-o", "csize", ffWRITE },
517 { efXPM, "-ow", "csizew", ffWRITE },
518 { efXVG, "-nc", "nclust", ffWRITE },
519 { efXVG, "-mc", "maxclust", ffWRITE },
520 { efXVG, "-ac", "avclust", ffWRITE },
521 { efXVG, "-hc", "histo-clust", ffWRITE },
522 { efXVG, "-temp", "temp", ffOPTWR },
523 { efNDX, "-mcn", "maxclust", ffOPTWR }
524 };
525 #define NFILE asize(fnm)
526
527 if (!parse_common_args(&argc, argv,
528 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
529 NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
530 {
531 return 0;
532 }
533
534 fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
535 rgblo.r = rlo[XX]; rgblo.g = rlo[YY]; rgblo.b = rlo[ZZ];
536 rgbhi.r = rhi[XX]; rgbhi.g = rhi[YY]; rgbhi.b = rhi[ZZ];
537
538 fnTPR = ftp2fn_null(efTPR, NFILE, fnm);
539 if (bMol && !fnTPR)
540 {
541 gmx_fatal(FARGS, "You need a tpr file for the -mol option");
542 }
543
544 clust_size(fnNDX, ftp2fn(efTRX, NFILE, fnm), opt2fn("-o", NFILE, fnm),
545 opt2fn("-ow", NFILE, fnm),
546 opt2fn("-nc", NFILE, fnm), opt2fn("-ac", NFILE, fnm),
547 opt2fn("-mc", NFILE, fnm), opt2fn("-hc", NFILE, fnm),
548 opt2fn("-temp", NFILE, fnm), opt2fn("-mcn", NFILE, fnm),
549 bMol, bPBC, fnTPR,
550 cutoff, nskip, nlevels, rgblo, rgbhi, ndf, oenv);
551
552 output_env_done(oenv);
553
554 return 0;
555 }
The ultimate molecular dynamics simulation package