1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
61 static void get_refx(output_env_t oenv, const char *trxfn, int nfitdim, int skip,
63 gmx_bool bMW, t_topology *top, int ePBC, rvec *x_ref)
65 int natoms, nfr_all, nfr, i, j, a, r, c, min_fr;
68 double tot_mass, msd, *srmsd, min_srmsd, srmsd_tot;
73 gmx_rmpbc_t gpbc = NULL;
80 natoms = read_first_x(oenv, &status, trxfn, &ti[nfr], &x, box);
84 for (a = 0; a < gnx; a++)
86 if (index[a] >= natoms)
88 gmx_fatal(FARGS, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[a]+1, natoms);
90 w_rls[a] = (bMW ? top->atoms.atom[index[a]].m : 1.0);
93 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
97 if (nfr_all % skip == 0)
99 gmx_rmpbc(gpbc, natoms, box, x);
101 for (i = 0; i < gnx; i++)
103 copy_rvec(x[index[i]], xi[nfr][i]);
105 reset_x(gnx, NULL, gnx, NULL, xi[nfr], w_rls);
115 while (read_next_x(oenv, status, &ti[nfr], x, box));
119 gmx_rmpbc_done(gpbc);
122 for (i = 0; i < nfr; i++)
124 printf("\rProcessing frame %d of %d", i, nfr);
125 for (j = i+1; j < nfr; j++)
127 calc_fit_R(nfitdim, gnx, w_rls, xi[i], xi[j], R);
130 for (a = 0; a < gnx; a++)
132 for (r = 0; r < DIM; r++)
135 for (c = 0; c < DIM; c++)
137 xf += R[r][c]*xi[j][a][c];
139 msd += w_rls[a]*sqr(xi[i][a][r] - xf);
143 srmsd[i] += sqrt(msd);
144 srmsd[j] += sqrt(msd);
151 min_srmsd = GMX_REAL_MAX;
155 for (i = 0; i < nfr; i++)
157 srmsd[i] /= (nfr - 1);
158 if (srmsd[i] < min_srmsd)
160 min_srmsd = srmsd[i];
164 srmsd_tot += srmsd[i];
168 printf("Average RMSD between all structures: %.3f\n", srmsd_tot/nfr);
169 printf("Structure with lowest RMSD to all others: time %g, av. RMSD %.3f\n",
172 for (a = 0; a < gnx; a++)
174 copy_rvec(xi[min_fr][a], x_ref[index[a]]);
180 int gmx_rotmat(int argc, char *argv[])
182 const char *desc[] = {
183 "[TT]g_rotmat[tt] plots the rotation matrix required for least squares fitting",
184 "a conformation onto the reference conformation provided with",
185 "[TT]-s[tt]. Translation is removed before fitting.",
186 "The output are the three vectors that give the new directions",
187 "of the x, y and z directions of the reference conformation,",
188 "for example: (zx,zy,zz) is the orientation of the reference",
189 "z-axis in the trajectory frame.",
191 "This tool is useful for, for instance,",
192 "determining the orientation of a molecule",
193 "at an interface, possibly on a trajectory produced with",
194 "[TT]trjconv -fit rotxy+transxy[tt] to remove the rotation",
195 "in the [IT]x-y[it] plane.",
197 "Option [TT]-ref[tt] determines a reference structure for fitting,",
198 "instead of using the structure from [TT]-s[tt]. The structure with",
199 "the lowest sum of RMSD's to all other structures is used.",
200 "Since the computational cost of this procedure grows with",
201 "the square of the number of frames, the [TT]-skip[tt] option",
202 "can be useful. A full fit or only a fit in the [IT]x-y[it] plane can",
205 "Option [TT]-fitxy[tt] fits in the [IT]x-y[it] plane before determining",
206 "the rotation matrix."
208 const char *reffit[] =
209 { NULL, "none", "xyz", "xy", NULL };
211 static gmx_bool bFitXY = FALSE, bMW = TRUE;
213 { "-ref", FALSE, etENUM, {reffit},
214 "Determine the optimal reference structure" },
215 { "-skip", FALSE, etINT, {&skip},
216 "Use every nr-th frame for [TT]-ref[tt]" },
217 { "-fitxy", FALSE, etBOOL, {&bFitXY},
218 "Fit the x/y rotation before determining the rotation" },
219 { "-mw", FALSE, etBOOL, {&bMW},
220 "Use mass weighted fitting" }
230 char *grpname, title[256];
232 gmx_rmpbc_t gpbc = NULL;
236 const char *leg[] = { "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz" };
237 #define NLEG asize(leg)
239 { efTRX, "-f", NULL, ffREAD },
240 { efTPS, NULL, NULL, ffREAD },
241 { efNDX, NULL, NULL, ffOPTRD },
242 { efXVG, NULL, "rotmat", ffWRITE }
244 #define NFILE asize(fnm)
246 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
247 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
252 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x_ref, NULL, box, bMW);
254 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
256 gmx_rmpbc(gpbc, top.atoms.nr, box, x_ref);
258 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
260 if (reffit[0][0] != 'n')
262 get_refx(oenv, ftp2fn(efTRX, NFILE, fnm), reffit[0][2] == 'z' ? 3 : 2, skip,
263 gnx, index, bMW, &top, ePBC, x_ref);
266 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
269 for (i = 0; i < gnx; i++)
271 if (index[i] >= natoms)
273 gmx_fatal(FARGS, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[i]+1, natoms);
275 w_rls[index[i]] = (bMW ? top.atoms.atom[index[i]].m : 1.0);
278 if (reffit[0][0] == 'n')
280 reset_x(gnx, index, natoms, NULL, x_ref, w_rls);
283 out = xvgropen(ftp2fn(efXVG, NFILE, fnm),
284 "Fit matrix", "Time (ps)", "", oenv);
285 xvgr_legend(out, NLEG, leg, oenv);
289 gmx_rmpbc(gpbc, natoms, box, x);
291 reset_x(gnx, index, natoms, NULL, x, w_rls);
295 do_fit_ndim(2, natoms, w_rls, x_ref, x);
298 calc_fit_R(DIM, natoms, w_rls, x_ref, x, R);
301 "%7g %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f\n",
303 R[XX][XX], R[XX][YY], R[XX][ZZ],
304 R[YY][XX], R[YY][YY], R[YY][ZZ],
305 R[ZZ][XX], R[ZZ][YY], R[ZZ][ZZ]);
307 while (read_next_x(oenv, status, &t, x, box));
309 gmx_rmpbc_done(gpbc);
315 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");