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
43 #include "gromacs/fileio/confio.h"
44 #include "gmx_fatal.h"
50 #include "gromacs/fileio/trxio.h"
53 static void rot_conf(t_atoms *atoms, rvec x[], rvec v[], real trans, real angle,
54 rvec head, rvec tail, int isize, atom_id index[],
55 rvec xout[], rvec vout[])
58 real theta, phi, arrow_len;
59 mat4 Rx, Ry, Rz, Rinvy, Rinvz, Mtot;
60 mat4 temp1, temp2, temp3;
64 rvec_sub(tail, head, arrow);
65 arrow_len = norm(arrow);
68 fprintf(debug, "Arrow vector: %10.4f %10.4f %10.4f\n",
69 arrow[XX], arrow[YY], arrow[ZZ]);
70 fprintf(debug, "Effective translation %g nm\n", trans);
74 gmx_fatal(FARGS, "Arrow vector not given");
77 /* Copy all aoms to output */
78 for (i = 0; (i < atoms->nr); i++)
80 copy_rvec(x[i], xout[i]);
81 copy_rvec(v[i], vout[i]);
84 /* Compute center of mass and move atoms there */
86 for (i = 0; (i < isize); i++)
88 rvec_inc(xcm, x[index[i]]);
90 for (i = 0; (i < DIM); i++)
96 fprintf(debug, "Center of mass: %10.4f %10.4f %10.4f\n",
97 xcm[XX], xcm[YY], xcm[ZZ]);
99 for (i = 0; (i < isize); i++)
101 rvec_sub(x[index[i]], xcm, xout[index[i]]);
104 /* Compute theta and phi that describe the arrow */
105 theta = acos(arrow[ZZ]/arrow_len);
106 phi = atan2(arrow[YY]/arrow_len, arrow[XX]/arrow_len);
109 fprintf(debug, "Phi = %.1f, Theta = %.1f\n", RAD2DEG*phi, RAD2DEG*theta);
112 /* Now the total rotation matrix: */
113 /* Rotate a couple of times */
114 rotate(ZZ, -phi, Rz);
115 rotate(YY, M_PI/2-theta, Ry);
116 rotate(XX, angle*DEG2RAD, Rx);
118 rotate(YY, theta-M_PI/2, Rinvy);
119 rotate(ZZ, phi, Rinvz);
121 mult_matrix(temp1, Ry, Rz);
122 mult_matrix(temp2, Rinvy, Rx);
123 mult_matrix(temp3, temp2, temp1);
124 mult_matrix(Mtot, Rinvz, temp3);
126 print_m4(debug, "Rz", Rz);
127 print_m4(debug, "Ry", Ry);
128 print_m4(debug, "Rx", Rx);
129 print_m4(debug, "Rinvy", Rinvy);
130 print_m4(debug, "Rinvz", Rinvz);
131 print_m4(debug, "Mtot", Mtot);
133 for (i = 0; (i < isize); i++)
136 m4_op(Mtot, xout[ai], xv);
137 rvec_add(xv, xcm, xout[ai]);
138 m4_op(Mtot, v[ai], xv);
139 copy_rvec(xv, vout[ai]);
143 int gmx_dyndom(int argc, char *argv[])
145 const char *desc[] = {
146 "[THISMODULE] reads a [TT].pdb[tt] file output from DynDom",
147 "(http://www.cmp.uea.ac.uk/dyndom/).",
148 "It reads the coordinates, the coordinates of the rotation axis,",
149 "and an index file containing the domains.",
150 "Furthermore, it takes the first and last atom of the arrow file",
151 "as command line arguments (head and tail) and",
152 "finally it takes the translation vector (given in DynDom info file)",
153 "and the angle of rotation (also as command line arguments). If the angle",
154 "determined by DynDom is given, one should be able to recover the",
155 "second structure used for generating the DynDom output.",
156 "Because of limited numerical accuracy this should be verified by",
157 "computing an all-atom RMSD (using [gmx-confrms]) rather than by file",
158 "comparison (using diff).[PAR]",
159 "The purpose of this program is to interpolate and extrapolate the",
160 "rotation as found by DynDom. As a result unphysical structures with",
161 "long or short bonds, or overlapping atoms may be produced. Visual",
162 "inspection, and energy minimization may be necessary to",
163 "validate the structure."
165 static real trans0 = 0;
166 static rvec head = { 0, 0, 0 };
167 static rvec tail = { 0, 0, 0 };
168 static real angle0 = 0, angle1 = 0, maxangle = 0;
169 static int label = 0, nframes = 11;
171 { "-firstangle", FALSE, etREAL, {&angle0},
172 "Angle of rotation about rotation vector" },
173 { "-lastangle", FALSE, etREAL, {&angle1},
174 "Angle of rotation about rotation vector" },
175 { "-nframe", FALSE, etINT, {&nframes},
176 "Number of steps on the pathway" },
177 { "-maxangle", FALSE, etREAL, {&maxangle},
178 "DymDom dtermined angle of rotation about rotation vector" },
179 { "-trans", FALSE, etREAL, {&trans0},
180 "Translation (Angstrom) along rotation vector (see DynDom info file)" },
181 { "-head", FALSE, etRVEC, {head},
182 "First atom of the arrow vector" },
183 { "-tail", FALSE, etRVEC, {tail},
184 "Last atom of the arrow vector" }
186 int i, j, natoms, isize;
188 atom_id *index = NULL, *index_all;
189 char title[256], *grpname;
192 rvec *x, *v, *xout, *vout;
197 { efPDB, "-f", "dyndom", ffREAD },
198 { efTRO, "-o", "rotated", ffWRITE },
199 { efNDX, "-n", "domains", ffREAD }
201 #define NFILE asize(fnm)
203 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
204 asize(desc), desc, 0, NULL, &oenv))
209 get_stx_coordnum (opt2fn("-f", NFILE, fnm), &natoms);
210 init_t_atoms(&atoms, natoms, TRUE);
213 read_stx_conf(opt2fn("-f", NFILE, fnm), title, &atoms, x, v, NULL, box);
217 printf("Select group to rotate:\n");
218 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
219 printf("Going to rotate %s containg %d atoms\n", grpname, isize);
221 snew(index_all, atoms.nr);
222 for (i = 0; (i < atoms.nr); i++)
227 status = open_trx(opt2fn("-o", NFILE, fnm), "w");
230 for (i = 0; (i < nframes); i++, label++)
232 angle = angle0 + (i*(angle1-angle0))/(nframes-1);
233 trans = trans0*0.1*angle/maxangle;
234 printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
235 i, label, angle, trans);
236 rot_conf(&atoms, x, v, trans, angle, head, tail, isize, index, xout, vout);
242 for (j = 0; (j < atoms.nr); j++)
244 atoms.resinfo[atoms.atom[j].resind].chainid = label;
247 write_trx(status, atoms.nr, index_all, &atoms, i, angle, box, xout, vout, NULL);