src/gromacs/gmxana/gmx_rotmat.cpp

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