2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, 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.
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.
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.
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.
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.
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.
39 #include "gromacs/commandline/pargs.h"
40 #include "gromacs/fileio/confio.h"
41 #include "gromacs/fileio/trxio.h"
42 #include "gromacs/gmxana/gmx_ana.h"
43 #include "gromacs/legacyheaders/macros.h"
44 #include "gromacs/math/3dtransforms.h"
45 #include "gromacs/math/units.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/topology/atoms.h"
48 #include "gromacs/topology/index.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/utility/smalloc.h"
52 static void rot_conf(t_atoms *atoms, rvec x[], rvec v[], real trans, real angle,
53 rvec head, rvec tail, int isize, atom_id index[],
54 rvec xout[], rvec vout[])
57 real theta, phi, arrow_len;
58 mat4 Rx, Ry, Rz, Rinvy, Rinvz, Mtot;
59 mat4 temp1, temp2, temp3;
63 rvec_sub(tail, head, arrow);
64 arrow_len = norm(arrow);
67 fprintf(debug, "Arrow vector: %10.4f %10.4f %10.4f\n",
68 arrow[XX], arrow[YY], arrow[ZZ]);
69 fprintf(debug, "Effective translation %g nm\n", trans);
73 gmx_fatal(FARGS, "Arrow vector not given");
76 /* Copy all aoms to output */
77 for (i = 0; (i < atoms->nr); i++)
79 copy_rvec(x[i], xout[i]);
80 copy_rvec(v[i], vout[i]);
83 /* Compute center of mass and move atoms there */
85 for (i = 0; (i < isize); i++)
87 rvec_inc(xcm, x[index[i]]);
89 for (i = 0; (i < DIM); i++)
95 fprintf(debug, "Center of mass: %10.4f %10.4f %10.4f\n",
96 xcm[XX], xcm[YY], xcm[ZZ]);
98 for (i = 0; (i < isize); i++)
100 rvec_sub(x[index[i]], xcm, xout[index[i]]);
103 /* Compute theta and phi that describe the arrow */
104 theta = acos(arrow[ZZ]/arrow_len);
105 phi = atan2(arrow[YY]/arrow_len, arrow[XX]/arrow_len);
108 fprintf(debug, "Phi = %.1f, Theta = %.1f\n", RAD2DEG*phi, RAD2DEG*theta);
111 /* Now the total rotation matrix: */
112 /* Rotate a couple of times */
113 gmx_mat4_init_rotation(ZZ, -phi, Rz);
114 gmx_mat4_init_rotation(YY, M_PI/2-theta, Ry);
115 gmx_mat4_init_rotation(XX, angle*DEG2RAD, Rx);
117 gmx_mat4_init_rotation(YY, theta-M_PI/2, Rinvy);
118 gmx_mat4_init_rotation(ZZ, phi, Rinvz);
120 gmx_mat4_mmul(temp1, Ry, Rz);
121 gmx_mat4_mmul(temp2, Rinvy, Rx);
122 gmx_mat4_mmul(temp3, temp2, temp1);
123 gmx_mat4_mmul(Mtot, Rinvz, temp3);
127 gmx_mat4_print(debug, "Rz", Rz);
128 gmx_mat4_print(debug, "Ry", Ry);
129 gmx_mat4_print(debug, "Rx", Rx);
130 gmx_mat4_print(debug, "Rinvy", Rinvy);
131 gmx_mat4_print(debug, "Rinvz", Rinvz);
132 gmx_mat4_print(debug, "Mtot", Mtot);
135 for (i = 0; (i < isize); i++)
138 gmx_mat4_transform_point(Mtot, xout[ai], xv);
139 rvec_add(xv, xcm, xout[ai]);
140 gmx_mat4_transform_point(Mtot, v[ai], xv);
141 copy_rvec(xv, vout[ai]);
145 int gmx_dyndom(int argc, char *argv[])
147 const char *desc[] = {
148 "[THISMODULE] reads a [TT].pdb[tt] file output from DynDom",
149 "(http://www.cmp.uea.ac.uk/dyndom/).",
150 "It reads the coordinates, the coordinates of the rotation axis,",
151 "and an index file containing the domains.",
152 "Furthermore, it takes the first and last atom of the arrow file",
153 "as command line arguments (head and tail) and",
154 "finally it takes the translation vector (given in DynDom info file)",
155 "and the angle of rotation (also as command line arguments). If the angle",
156 "determined by DynDom is given, one should be able to recover the",
157 "second structure used for generating the DynDom output.",
158 "Because of limited numerical accuracy this should be verified by",
159 "computing an all-atom RMSD (using [gmx-confrms]) rather than by file",
160 "comparison (using diff).[PAR]",
161 "The purpose of this program is to interpolate and extrapolate the",
162 "rotation as found by DynDom. As a result unphysical structures with",
163 "long or short bonds, or overlapping atoms may be produced. Visual",
164 "inspection, and energy minimization may be necessary to",
165 "validate the structure."
167 static real trans0 = 0;
168 static rvec head = { 0, 0, 0 };
169 static rvec tail = { 0, 0, 0 };
170 static real angle0 = 0, angle1 = 0, maxangle = 0;
171 static int label = 0, nframes = 11;
173 { "-firstangle", FALSE, etREAL, {&angle0},
174 "Angle of rotation about rotation vector" },
175 { "-lastangle", FALSE, etREAL, {&angle1},
176 "Angle of rotation about rotation vector" },
177 { "-nframe", FALSE, etINT, {&nframes},
178 "Number of steps on the pathway" },
179 { "-maxangle", FALSE, etREAL, {&maxangle},
180 "DymDom dtermined angle of rotation about rotation vector" },
181 { "-trans", FALSE, etREAL, {&trans0},
182 "Translation (Angstrom) along rotation vector (see DynDom info file)" },
183 { "-head", FALSE, etRVEC, {head},
184 "First atom of the arrow vector" },
185 { "-tail", FALSE, etRVEC, {tail},
186 "Last atom of the arrow vector" }
188 int i, j, natoms, isize;
190 atom_id *index = NULL, *index_all;
191 char title[256], *grpname;
194 rvec *x, *v, *xout, *vout;
199 { efPDB, "-f", "dyndom", ffREAD },
200 { efTRO, "-o", "rotated", ffWRITE },
201 { efNDX, "-n", "domains", ffREAD }
203 #define NFILE asize(fnm)
205 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
206 asize(desc), desc, 0, NULL, &oenv))
211 get_stx_coordnum (opt2fn("-f", NFILE, fnm), &natoms);
212 init_t_atoms(&atoms, natoms, TRUE);
215 read_stx_conf(opt2fn("-f", NFILE, fnm), title, &atoms, x, v, NULL, box);
219 printf("Select group to rotate:\n");
220 rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
221 printf("Going to rotate %s containg %d atoms\n", grpname, isize);
223 snew(index_all, atoms.nr);
224 for (i = 0; (i < atoms.nr); i++)
229 status = open_trx(opt2fn("-o", NFILE, fnm), "w");
232 for (i = 0; (i < nframes); i++, label++)
234 angle = angle0 + (i*(angle1-angle0))/(nframes-1);
235 trans = trans0*0.1*angle/maxangle;
236 printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
237 i, label, angle, trans);
238 rot_conf(&atoms, x, v, trans, angle, head, tail, isize, index, xout, vout);
244 for (j = 0; (j < atoms.nr); j++)
246 atoms.resinfo[atoms.atom[j].resind].chainid = label;
249 write_trx(status, atoms.nr, index_all, &atoms, i, angle, box, xout, vout, NULL);