src/gromacs/gmxana/gmx_rotmat.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmxpre.h"
  36
  37 #include <math.h>
  38 #include <string.h>
  39
  40 #include "gromacs/commandline/pargs.h"
  41 #include "gromacs/fileio/tpxio.h"
  42 #include "gromacs/fileio/trxio.h"
  43 #include "gromacs/fileio/xvgr.h"
  44 #include "gromacs/gmxana/gmx_ana.h"
  45 #include "gromacs/legacyheaders/macros.h"
  46 #include "gromacs/legacyheaders/typedefs.h"
  47 #include "gromacs/legacyheaders/viewit.h"
  48 #include "gromacs/math/do_fit.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/rmpbc.h"
  51 #include "gromacs/topology/index.h"
  52 #include "gromacs/utility/futil.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 static void get_refx(output_env_t oenv, const char *trxfn, int nfitdim, int skip,
  56                      int gnx, int *index,
  57                      gmx_bool bMW, t_topology *top, int ePBC, rvec *x_ref)
  58 {
  59     int          natoms, nfr_all, nfr, i, j, a, r, c, min_fr;
  60     t_trxstatus *status;
  61     real        *ti, min_t;
  62     double       tot_mass, msd, *srmsd, min_srmsd, srmsd_tot;
  63     rvec        *x, **xi;
  64     real         xf;
  65     matrix       box, R;
  66     real        *w_rls;
  67     gmx_rmpbc_t  gpbc = NULL;
  68
  69
  70     nfr_all = 0;
  71     nfr     = 0;
  72     snew(ti, 100);
  73     snew(xi, 100);
  74     natoms = read_first_x(oenv, &status, trxfn, &ti[nfr], &x, box);
  75
  76     snew(w_rls, gnx);
  77     tot_mass = 0;
  78     for (a = 0; a < gnx; a++)
  79     {
  80         if (index[a] >= natoms)
  81         {
  82             gmx_fatal(FARGS, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[a]+1, natoms);
  83         }
  84         w_rls[a]  = (bMW ? top->atoms.atom[index[a]].m : 1.0);
  85         tot_mass += w_rls[a];
  86     }
  87     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
  88
  89     do
  90     {
  91         if (nfr_all % skip == 0)
  92         {
  93             gmx_rmpbc(gpbc, natoms, box, x);
  94             snew(xi[nfr], gnx);
  95             for (i = 0; i < gnx; i++)
  96             {
  97                 copy_rvec(x[index[i]], xi[nfr][i]);
  98             }
  99             reset_x(gnx, NULL, gnx, NULL, xi[nfr], w_rls);
 100             nfr++;
 101             if (nfr % 100 == 0)
 102             {
 103                 srenew(ti, nfr+100);
 104                 srenew(xi, nfr+100);
 105             }
 106         }
 107         nfr_all++;
 108     }
 109     while (read_next_x(oenv, status, &ti[nfr], x, box));
 110     close_trj(status);
 111     sfree(x);
 112
 113     gmx_rmpbc_done(gpbc);
 114
 115     snew(srmsd, nfr);
 116     for (i = 0; i < nfr; i++)
 117     {
 118         printf("\rProcessing frame %d of %d", i, nfr);
 119         for (j = i+1; j < nfr; j++)
 120         {
 121             calc_fit_R(nfitdim, gnx, w_rls, xi[i], xi[j], R);
 122
 123             msd = 0;
 124             for (a = 0; a < gnx; a++)
 125             {
 126                 for (r = 0; r < DIM; r++)
 127                 {
 128                     xf = 0;
 129                     for (c = 0; c < DIM; c++)
 130                     {
 131                         xf += R[r][c]*xi[j][a][c];
 132                     }
 133                     msd += w_rls[a]*sqr(xi[i][a][r] - xf);
 134                 }
 135             }
 136             msd      /= tot_mass;
 137             srmsd[i] += sqrt(msd);
 138             srmsd[j] += sqrt(msd);
 139         }
 140         sfree(xi[i]);
 141     }
 142     printf("\n");
 143     sfree(w_rls);
 144
 145     min_srmsd = GMX_REAL_MAX;
 146     min_fr    = -1;
 147     min_t     = -1;
 148     srmsd_tot = 0;
 149     for (i = 0; i < nfr; i++)
 150     {
 151         srmsd[i] /= (nfr - 1);
 152         if (srmsd[i] < min_srmsd)
 153         {
 154             min_srmsd = srmsd[i];
 155             min_fr    = i;
 156             min_t     = ti[i];
 157         }
 158         srmsd_tot += srmsd[i];
 159     }
 160     sfree(srmsd);
 161
 162     printf("Average RMSD between all structures: %.3f\n", srmsd_tot/nfr);
 163     printf("Structure with lowest RMSD to all others: time %g, av. RMSD %.3f\n",
 164            min_t, min_srmsd);
 165
 166     for (a = 0; a < gnx; a++)
 167     {
 168         copy_rvec(xi[min_fr][a], x_ref[index[a]]);
 169     }
 170
 171     sfree(xi);
 172 }
 173
 174 int gmx_rotmat(int argc, char *argv[])
 175 {
 176     const char     *desc[] = {
 177         "[THISMODULE] plots the rotation matrix required for least squares fitting",
 178         "a conformation onto the reference conformation provided with",
 179         "[TT]-s[tt]. Translation is removed before fitting.",
 180         "The output are the three vectors that give the new directions",
 181         "of the x, y and z directions of the reference conformation,",
 182         "for example: (zx,zy,zz) is the orientation of the reference",
 183         "z-axis in the trajectory frame.",
 184         "[PAR]",
 185         "This tool is useful for, for instance,",
 186         "determining the orientation of a molecule",
 187         "at an interface, possibly on a trajectory produced with",
 188         "[TT]gmx trjconv -fit rotxy+transxy[tt] to remove the rotation",
 189         "in the [IT]x-y[it] plane.",
 190         "[PAR]",
 191         "Option [TT]-ref[tt] determines a reference structure for fitting,",
 192         "instead of using the structure from [TT]-s[tt]. The structure with",
 193         "the lowest sum of RMSD's to all other structures is used.",
 194         "Since the computational cost of this procedure grows with",
 195         "the square of the number of frames, the [TT]-skip[tt] option",
 196         "can be useful. A full fit or only a fit in the [IT]x-y[it] plane can",
 197         "be performed.",
 198         "[PAR]",
 199         "Option [TT]-fitxy[tt] fits in the [IT]x-y[it] plane before determining",
 200         "the rotation matrix."
 201     };
 202     const char     *reffit[] =
 203     { NULL, "none", "xyz", "xy", NULL };
 204     static int      skip   = 1;
 205     static gmx_bool bFitXY = FALSE, bMW = TRUE;
 206     t_pargs         pa[]   = {
 207         { "-ref", FALSE, etENUM, {reffit},
 208           "Determine the optimal reference structure" },
 209         { "-skip", FALSE, etINT, {&skip},
 210           "Use every nr-th frame for [TT]-ref[tt]" },
 211         { "-fitxy", FALSE, etBOOL, {&bFitXY},
 212           "Fit the x/y rotation before determining the rotation" },
 213         { "-mw", FALSE, etBOOL, {&bMW},
 214           "Use mass weighted fitting" }
 215     };
 216     FILE           *out;
 217     t_trxstatus    *status;
 218     t_topology      top;
 219     int             ePBC;
 220     rvec           *x_ref, *x;
 221     matrix          box, R;
 222     real            t;
 223     int             natoms, i;
 224     char           *grpname, title[256];
 225     int             gnx;
 226     gmx_rmpbc_t     gpbc = NULL;
 227     atom_id        *index;
 228     output_env_t    oenv;
 229     real           *w_rls;
 230     const char     *leg[]  = { "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz" };
 231 #define NLEG asize(leg)
 232     t_filenm        fnm[] = {
 233         { efTRX, "-f",   NULL,       ffREAD },
 234         { efTPS, NULL,   NULL,       ffREAD },
 235         { efNDX, NULL,   NULL,       ffOPTRD },
 236         { efXVG, NULL,   "rotmat",   ffWRITE }
 237     };
 238 #define NFILE asize(fnm)
 239
 240     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 241                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 242     {
 243         return 0;
 244     }
 245
 246     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x_ref, NULL, box, bMW);
 247
 248     gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
 249
 250     gmx_rmpbc(gpbc, top.atoms.nr, box, x_ref);
 251
 252     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 253
 254     if (reffit[0][0] != 'n')
 255     {
 256         get_refx(oenv, ftp2fn(efTRX, NFILE, fnm), reffit[0][2] == 'z' ? 3 : 2, skip,
 257                  gnx, index, bMW, &top, ePBC, x_ref);
 258     }
 259
 260     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 261
 262     snew(w_rls, natoms);
 263     for (i = 0; i < gnx; i++)
 264     {
 265         if (index[i] >= natoms)
 266         {
 267             gmx_fatal(FARGS, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[i]+1, natoms);
 268         }
 269         w_rls[index[i]] = (bMW ? top.atoms.atom[index[i]].m : 1.0);
 270     }
 271
 272     if (reffit[0][0] == 'n')
 273     {
 274         reset_x(gnx, index, natoms, NULL, x_ref, w_rls);
 275     }
 276
 277     out = xvgropen(ftp2fn(efXVG, NFILE, fnm),
 278                    "Fit matrix", "Time (ps)", "", oenv);
 279     xvgr_legend(out, NLEG, leg, oenv);
 280
 281     do
 282     {
 283         gmx_rmpbc(gpbc, natoms, box, x);
 284
 285         reset_x(gnx, index, natoms, NULL, x, w_rls);
 286
 287         if (bFitXY)
 288         {
 289             do_fit_ndim(2, natoms, w_rls, x_ref, x);
 290         }
 291
 292         calc_fit_R(DIM, natoms, w_rls, x_ref, x, R);
 293
 294         fprintf(out,
 295                 "%7g %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f\n",
 296                 t,
 297                 R[XX][XX], R[XX][YY], R[XX][ZZ],
 298                 R[YY][XX], R[YY][YY], R[YY][ZZ],
 299                 R[ZZ][XX], R[ZZ][YY], R[ZZ][ZZ]);
 300     }
 301     while (read_next_x(oenv, status, &t, x, box));
 302
 303     gmx_rmpbc_done(gpbc);
 304
 305     close_trj(status);
 306
 307     gmx_ffclose(out);
 308
 309     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
 310
 311     return 0;
 312 }