3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
60 static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
61 real *rmin, real *rmax, int *min_ind)
64 int sx, sy, sz, i, j, s;
65 real sqr_box, r2min, r2max, r2;
66 rvec shift[NSHIFT], d0, d;
68 sqr_box = sqr(min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ])));
71 for (sz = -1; sz <= 1; sz++)
73 for (sy = -1; sy <= 1; sy++)
75 for (sx = -1; sx <= 1; sx++)
77 if (sx != 0 || sy != 0 || sz != 0)
79 for (i = 0; i < DIM; i++)
81 shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
92 for (i = 0; i < n; i++)
94 for (j = i+1; j < n; j++)
96 rvec_sub(x[index[i]], x[index[j]], d0);
102 for (s = 0; s < NSHIFT; s++)
104 rvec_add(d0, shift[s], d);
120 static void periodic_mindist_plot(const char *trxfn, const char *outfn,
121 t_topology *top, int ePBC,
122 int n, atom_id index[], gmx_bool bSplit,
123 const output_env_t oenv)
126 const char *leg[5] = { "min per.", "max int.", "box1", "box2", "box3" };
131 int natoms, ind_min[2] = {0, 0}, ind_mini = 0, ind_minj = 0;
132 real r, rmin, rmax, rmint, tmint;
134 gmx_rmpbc_t gpbc = NULL;
136 natoms = read_first_x(oenv, &status, trxfn, &t, &x, box);
138 check_index(NULL, n, index, NULL, natoms);
140 out = xvgropen(outfn, "Minimum distance to periodic image",
141 output_env_get_time_label(oenv), "Distance (nm)", oenv);
142 if (output_env_get_print_xvgr_codes(oenv))
144 fprintf(out, "@ subtitle \"and maximum internal distance\"\n");
146 xvgr_legend(out, 5, leg, oenv);
153 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
161 gmx_rmpbc(gpbc, natoms, box, x);
164 periodic_dist(box, x, n, index, &rmin, &rmax, ind_min);
169 ind_mini = ind_min[0];
170 ind_minj = ind_min[1];
172 if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
176 fprintf(out, "\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
177 output_env_conv_time(oenv, t), rmin, rmax, norm(box[0]), norm(box[1]), norm(box[2]));
180 while (read_next_x(oenv, status, &t, natoms, x, box));
184 gmx_rmpbc_done(gpbc);
190 "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
191 "between atoms %d and %d\n",
192 rmint, output_env_conv_time(oenv, tmint), output_env_get_time_unit(oenv),
193 index[ind_mini]+1, index[ind_minj]+1);
196 static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[],
197 int nx1, int nx2, atom_id index1[], atom_id index2[],
199 real *rmin, real *rmax, int *nmin, int *nmax,
200 int *ixmin, int *jxmin, int *ixmax, int *jxmax)
202 int i, j, i0 = 0, j1;
206 real r2, rmin2, rmax2, rcut2;
219 /* Must init pbc every step because of pressure coupling */
222 set_pbc(&pbc, ePBC, box);
239 for (j = 0; (j < j1); j++)
248 for (i = i0; (i < nx1); i++)
255 pbc_dx(&pbc, x[ix], x[jx], dx);
259 rvec_sub(x[ix], x[jx], dx);
305 void dist_plot(const char *fn, const char *afile, const char *dfile,
306 const char *nfile, const char *rfile, const char *xfile,
307 real rcut, gmx_bool bMat, t_atoms *atoms,
308 int ng, atom_id *index[], int gnx[], char *grpn[], gmx_bool bSplit,
309 gmx_bool bMin, int nres, atom_id *residue, gmx_bool bPBC, int ePBC,
310 gmx_bool bGroup, gmx_bool bEachResEachTime, gmx_bool bPrintResName,
311 const output_env_t oenv)
313 FILE *atm, *dist, *num;
317 real t, dmin, dmax, **mindres = NULL, **maxdres = NULL;
320 int i = -1, j, k, natoms;
321 int min1, min2, max1, max2, min1r, min2r, max1r, max2r;
327 FILE *respertime = NULL;
329 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
331 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
334 sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
335 dist = xvgropen(dfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
336 sprintf(buf, "Number of Contacts %s %g nm", bMin ? "<" : ">", rcut);
337 num = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : NULL;
338 atm = afile ? ffopen(afile, "w") : NULL;
339 trxout = xfile ? open_trx(xfile, "w") : NULL;
346 sprintf(buf, "Internal in %s", grpn[0]);
347 leg[0] = strdup(buf);
348 xvgr_legend(dist, 0, (const char**)leg, oenv);
351 xvgr_legend(num, 0, (const char**)leg, oenv);
356 snew(leg, (ng*(ng-1))/2);
357 for (i = j = 0; (i < ng-1); i++)
359 for (k = i+1; (k < ng); k++, j++)
361 sprintf(buf, "%s-%s", grpn[i], grpn[k]);
362 leg[j] = strdup(buf);
365 xvgr_legend(dist, j, (const char**)leg, oenv);
368 xvgr_legend(num, j, (const char**)leg, oenv);
375 for (i = 0; (i < ng-1); i++)
377 sprintf(buf, "%s-%s", grpn[0], grpn[i+1]);
378 leg[i] = strdup(buf);
380 xvgr_legend(dist, ng-1, (const char**)leg, oenv);
383 xvgr_legend(num, ng-1, (const char**)leg, oenv);
387 if (bEachResEachTime)
389 sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
390 respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
391 xvgr_legend(respertime, ng-1, (const char**)leg, oenv);
394 fprintf(respertime, "# ");
396 for (j = 0; j < nres; j++)
398 fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind);
400 fprintf(respertime, "\n");
409 for (i = 1; i < ng; i++)
411 snew(mindres[i-1], nres);
412 snew(maxdres[i-1], nres);
413 for (j = 0; j < nres; j++)
415 mindres[i-1][j] = 1e6;
417 /* maxdres[*][*] is already 0 */
423 if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
425 fprintf(dist, "&\n");
435 fprintf(dist, "%12e", output_env_conv_time(oenv, t));
438 fprintf(num, "%12e", output_env_conv_time(oenv, t));
445 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[0], index[0], index[0], bGroup,
446 &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
447 fprintf(dist, " %12e", bMin ? dmin : dmax);
450 fprintf(num, " %8d", bMin ? nmin : nmax);
455 for (i = 0; (i < ng-1); i++)
457 for (k = i+1; (k < ng); k++)
459 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[i], gnx[k], index[i], index[k],
460 bGroup, &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
461 fprintf(dist, " %12e", bMin ? dmin : dmax);
464 fprintf(num, " %8d", bMin ? nmin : nmax);
472 for (i = 1; (i < ng); i++)
474 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[i], index[0], index[i], bGroup,
475 &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
476 fprintf(dist, " %12e", bMin ? dmin : dmax);
479 fprintf(num, " %8d", bMin ? nmin : nmax);
483 for (j = 0; j < nres; j++)
485 calc_dist(rcut, bPBC, ePBC, box, x0, residue[j+1]-residue[j], gnx[i],
486 &(index[0][residue[j]]), index[i], bGroup,
487 &dmin, &dmax, &nmin, &nmax, &min1r, &min2r, &max1r, &max2r);
488 mindres[i-1][j] = min(mindres[i-1][j], dmin);
489 maxdres[i-1][j] = max(maxdres[i-1][j], dmax);
499 if ( (bMin ? min1 : max1) != -1)
503 fprintf(atm, "%12e %12d %12d\n",
504 output_env_conv_time(oenv, t), 1+(bMin ? min1 : max1),
505 1+(bMin ? min2 : max2));
511 oindex[0] = bMin ? min1 : max1;
512 oindex[1] = bMin ? min2 : max2;
513 write_trx(trxout, 2, oindex, atoms, i, t, box, x0, NULL, NULL);
516 /*dmin should be minimum distance for residue and group*/
517 if (bEachResEachTime)
519 fprintf(respertime, "%12e", t);
520 for (i = 1; i < ng; i++)
522 for (j = 0; j < nres; j++)
524 fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
525 /*reset distances for next time point*/
526 mindres[i-1][j] = 1e6;
530 fprintf(respertime, "\n");
533 while (read_next_x(oenv, status, &t, natoms, x0, box));
550 if (nres && !bEachResEachTime)
554 sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
555 res = xvgropen(rfile, buf, "Residue (#)", "Distance (nm)", oenv);
556 xvgr_legend(res, ng-1, (const char**)leg, oenv);
557 for (j = 0; j < nres; j++)
559 fprintf(res, "%4d", j+1);
560 for (i = 1; i < ng; i++)
562 fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
571 int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
574 int nres = 0, resnr, presnr;
577 /* build index of first atom numbers for each residue */
579 snew(residx, atoms->nres+1);
580 for (i = 0; i < n; i++)
582 resnr = atoms->atom[index[i]].resind;
592 printf("Found %d residues out of %d (%d/%d atoms)\n",
593 nres, atoms->nres, atoms->nr, n);
595 srenew(residx, nres+1);
596 /* mark end of last residue */
602 void dump_res(FILE *out, int nres, atom_id *resindex, int n, atom_id index[])
606 for (i = 0; i < nres-1; i++)
608 fprintf(out, "Res %d (%d):", i, resindex[i+1]-resindex[i]);
609 for (j = resindex[i]; j < resindex[i+1]; j++)
611 fprintf(out, " %d(%d)", j, index[j]);
617 int gmx_mindist(int argc, char *argv[])
619 const char *desc[] = {
620 "[TT]g_mindist[tt] computes the distance between one group and a number of",
621 "other groups. Both the minimum distance",
622 "(between any pair of atoms from the respective groups)",
623 "and the number of contacts within a given",
624 "distance are written to two separate output files.",
625 "With the [TT]-group[tt] option a contact of an atom in another group",
626 "with multiple atoms in the first group is counted as one contact",
627 "instead of as multiple contacts.",
628 "With [TT]-or[tt], minimum distances to each residue in the first",
629 "group are determined and plotted as a function of residue number.[PAR]",
630 "With option [TT]-pi[tt] the minimum distance of a group to its",
631 "periodic image is plotted. This is useful for checking if a protein",
632 "has seen its periodic image during a simulation. Only one shift in",
633 "each direction is considered, giving a total of 26 shifts.",
634 "It also plots the maximum distance within the group and the lengths",
635 "of the three box vectors.[PAR]",
636 "Other programs that calculate distances are [TT]g_dist[tt]",
637 "and [TT]g_bond[tt]."
639 const char *bugs[] = {
640 "The [TT]-pi[tt] option is very slow."
643 static gmx_bool bMat = FALSE, bPI = FALSE, bSplit = FALSE, bMax = FALSE, bPBC = TRUE;
644 static gmx_bool bGroup = FALSE;
645 static real rcutoff = 0.6;
647 static gmx_bool bEachResEachTime = FALSE, bPrintResName = FALSE;
649 { "-matrix", FALSE, etBOOL, {&bMat},
650 "Calculate half a matrix of group-group distances" },
651 { "-max", FALSE, etBOOL, {&bMax},
652 "Calculate *maximum* distance instead of minimum" },
653 { "-d", FALSE, etREAL, {&rcutoff},
654 "Distance for contacts" },
655 { "-group", FALSE, etBOOL, {&bGroup},
656 "Count contacts with multiple atoms in the first group as one" },
657 { "-pi", FALSE, etBOOL, {&bPI},
658 "Calculate minimum distance with periodic images" },
659 { "-split", FALSE, etBOOL, {&bSplit},
660 "Split graph where time is zero" },
661 { "-ng", FALSE, etINT, {&ng},
662 "Number of secondary groups to compute distance to a central group" },
663 { "-pbc", FALSE, etBOOL, {&bPBC},
664 "Take periodic boundary conditions into account" },
665 { "-respertime", FALSE, etBOOL, {&bEachResEachTime},
666 "When writing per-residue distances, write distance for each time point" },
667 { "-printresname", FALSE, etBOOL, {&bPrintResName},
668 "Write residue names" }
671 t_topology *top = NULL;
677 gmx_bool bTop = FALSE;
681 const char *trxfnm, *tpsfnm, *ndxfnm, *distfnm, *numfnm, *atmfnm, *oxfnm, *resfnm;
684 atom_id **index, *residues = NULL;
686 { efTRX, "-f", NULL, ffREAD },
687 { efTPS, NULL, NULL, ffOPTRD },
688 { efNDX, NULL, NULL, ffOPTRD },
689 { efXVG, "-od", "mindist", ffWRITE },
690 { efXVG, "-on", "numcont", ffOPTWR },
691 { efOUT, "-o", "atm-pair", ffOPTWR },
692 { efTRO, "-ox", "mindist", ffOPTWR },
693 { efXVG, "-or", "mindistres", ffOPTWR }
695 #define NFILE asize(fnm)
697 parse_common_args(&argc, argv,
698 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
699 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
701 trxfnm = ftp2fn(efTRX, NFILE, fnm);
702 ndxfnm = ftp2fn_null(efNDX, NFILE, fnm);
703 distfnm = opt2fn("-od", NFILE, fnm);
704 numfnm = opt2fn_null("-on", NFILE, fnm);
705 atmfnm = ftp2fn_null(efOUT, NFILE, fnm);
706 oxfnm = opt2fn_null("-ox", NFILE, fnm);
707 resfnm = opt2fn_null("-or", NFILE, fnm);
708 if (bPI || resfnm != NULL)
710 /* We need a tps file */
711 tpsfnm = ftp2fn(efTPS, NFILE, fnm);
715 tpsfnm = ftp2fn_null(efTPS, NFILE, fnm);
718 if (!tpsfnm && !ndxfnm)
720 gmx_fatal(FARGS, "You have to specify either the index file or a tpr file");
726 fprintf(stderr, "Choose a group for distance calculation\n");
737 if (tpsfnm || resfnm || !ndxfnm)
740 bTop = read_tps_conf(tpsfnm, title, top, &ePBC, &x, NULL, box, FALSE);
743 printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
746 get_index(top ? &(top->atoms) : NULL, ndxfnm, ng, gnx, index, grpname);
748 if (bMat && (ng == 1))
751 printf("Special case: making distance matrix between all atoms in group %s\n",
756 for (i = 1; (i < ng); i++)
759 grpname[i] = grpname[0];
761 index[i][0] = index[0][i];
768 nres = find_residues(top ? &(top->atoms) : NULL,
769 gnx[0], index[0], &residues);
772 dump_res(debug, nres, residues, gnx[0], index[0]);
778 periodic_mindist_plot(trxfnm, distfnm, top, ePBC, gnx[0], index[0], bSplit, oenv);
782 dist_plot(trxfnm, atmfnm, distfnm, numfnm, resfnm, oxfnm,
783 rcutoff, bMat, top ? &(top->atoms) : NULL,
784 ng, index, gnx, grpname, bSplit, !bMax, nres, residues, bPBC, ePBC,
785 bGroup, bEachResEachTime, bPrintResName, oenv);
788 do_view(oenv, distfnm, "-nxy");
791 do_view(oenv, numfnm, "-nxy");