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");