src/gromacs/gmxana/gmx_dyndom.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015, 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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40
  41 #include "gromacs/commandline/pargs.h"
  42 #include "gromacs/fileio/confio.h"
  43 #include "gromacs/fileio/trxio.h"
  44 #include "gromacs/gmxana/gmx_ana.h"
  45 #include "gromacs/math/3dtransforms.h"
  46 #include "gromacs/math/units.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/topology/atoms.h"
  49 #include "gromacs/topology/index.h"
  50 #include "gromacs/topology/topology.h"
  51 #include "gromacs/utility/arraysize.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 static void rot_conf(t_atoms *atoms, rvec x[], rvec v[], real trans, real angle,
  56                      rvec head, rvec tail, int isize, int index[],
  57                      rvec xout[], rvec vout[])
  58 {
  59     rvec     arrow, xcm;
  60     real     theta, phi, arrow_len;
  61     mat4     Rx, Ry, Rz, Rinvy, Rinvz, Mtot;
  62     mat4     temp1, temp2, temp3;
  63     vec4     xv;
  64     int      i, ai;
  65
  66     rvec_sub(tail, head, arrow);
  67     arrow_len = norm(arrow);
  68     if (debug)
  69     {
  70         fprintf(debug, "Arrow vector:   %10.4f  %10.4f  %10.4f\n",
  71                 arrow[XX], arrow[YY], arrow[ZZ]);
  72         fprintf(debug, "Effective translation %g nm\n", trans);
  73     }
  74     if (arrow_len == 0.0)
  75     {
  76         gmx_fatal(FARGS, "Arrow vector not given");
  77     }
  78
  79     /* Copy all aoms to output */
  80     for (i = 0; (i < atoms->nr); i++)
  81     {
  82         copy_rvec(x[i], xout[i]);
  83         copy_rvec(v[i], vout[i]);
  84     }
  85
  86     /* Compute center of mass and move atoms there */
  87     clear_rvec(xcm);
  88     for (i = 0; (i < isize); i++)
  89     {
  90         rvec_inc(xcm, x[index[i]]);
  91     }
  92     for (i = 0; (i < DIM); i++)
  93     {
  94         xcm[i] /= isize;
  95     }
  96     if (debug)
  97     {
  98         fprintf(debug, "Center of mass: %10.4f  %10.4f  %10.4f\n",
  99                 xcm[XX], xcm[YY], xcm[ZZ]);
 100     }
 101     for (i = 0; (i < isize); i++)
 102     {
 103         rvec_sub(x[index[i]], xcm, xout[index[i]]);
 104     }
 105
 106     /* Compute theta and phi that describe the arrow */
 107     theta = std::acos(arrow[ZZ]/arrow_len);
 108     phi   = std::atan2(arrow[YY]/arrow_len, arrow[XX]/arrow_len);
 109     if (debug)
 110     {
 111         fprintf(debug, "Phi = %.1f, Theta = %.1f\n", RAD2DEG*phi, RAD2DEG*theta);
 112     }
 113
 114     /* Now the total rotation matrix: */
 115     /* Rotate a couple of times */
 116     gmx_mat4_init_rotation(ZZ, -phi, Rz);
 117     gmx_mat4_init_rotation(YY, M_PI/2-theta, Ry);
 118     gmx_mat4_init_rotation(XX, angle*DEG2RAD, Rx);
 119     Rx[WW][XX] = trans;
 120     gmx_mat4_init_rotation(YY, theta-M_PI/2, Rinvy);
 121     gmx_mat4_init_rotation(ZZ, phi, Rinvz);
 122
 123     gmx_mat4_mmul(temp1, Ry, Rz);
 124     gmx_mat4_mmul(temp2, Rinvy, Rx);
 125     gmx_mat4_mmul(temp3, temp2, temp1);
 126     gmx_mat4_mmul(Mtot, Rinvz, temp3);
 127
 128     if (debug)
 129     {
 130         gmx_mat4_print(debug, "Rz", Rz);
 131         gmx_mat4_print(debug, "Ry", Ry);
 132         gmx_mat4_print(debug, "Rx", Rx);
 133         gmx_mat4_print(debug, "Rinvy", Rinvy);
 134         gmx_mat4_print(debug, "Rinvz", Rinvz);
 135         gmx_mat4_print(debug, "Mtot", Mtot);
 136     }
 137
 138     for (i = 0; (i < isize); i++)
 139     {
 140         ai = index[i];
 141         gmx_mat4_transform_point(Mtot, xout[ai], xv);
 142         rvec_add(xv, xcm, xout[ai]);
 143         gmx_mat4_transform_point(Mtot, v[ai], xv);
 144         copy_rvec(xv, vout[ai]);
 145     }
 146 }
 147
 148 int gmx_dyndom(int argc, char *argv[])
 149 {
 150     const char       *desc[] = {
 151         "[THISMODULE] reads a [REF].pdb[ref] file output from DynDom",
 152         "(http://www.cmp.uea.ac.uk/dyndom/).",
 153         "It reads the coordinates, the coordinates of the rotation axis,",
 154         "and an index file containing the domains.",
 155         "Furthermore, it takes the first and last atom of the arrow file",
 156         "as command line arguments (head and tail) and",
 157         "finally it takes the translation vector (given in DynDom info file)",
 158         "and the angle of rotation (also as command line arguments). If the angle",
 159         "determined by DynDom is given, one should be able to recover the",
 160         "second structure used for generating the DynDom output.",
 161         "Because of limited numerical accuracy this should be verified by",
 162         "computing an all-atom RMSD (using [gmx-confrms]) rather than by file",
 163         "comparison (using diff).[PAR]",
 164         "The purpose of this program is to interpolate and extrapolate the",
 165         "rotation as found by DynDom. As a result unphysical structures with",
 166         "long or short bonds, or overlapping atoms may be produced. Visual",
 167         "inspection, and energy minimization may be necessary to",
 168         "validate the structure."
 169     };
 170     static real       trans0 = 0;
 171     static rvec       head   = { 0, 0, 0 };
 172     static rvec       tail   = { 0, 0, 0 };
 173     static real       angle0 = 0, angle1 = 0, maxangle = 0;
 174     static int        label  = 0, nframes = 11;
 175     t_pargs           pa[]   = {
 176         { "-firstangle",    FALSE, etREAL, {&angle0},
 177           "Angle of rotation about rotation vector" },
 178         { "-lastangle",    FALSE, etREAL, {&angle1},
 179           "Angle of rotation about rotation vector" },
 180         { "-nframe",   FALSE, etINT,  {&nframes},
 181           "Number of steps on the pathway" },
 182         { "-maxangle", FALSE, etREAL, {&maxangle},
 183           "DymDom dtermined angle of rotation about rotation vector" },
 184         { "-trans",    FALSE, etREAL, {&trans0},
 185           "Translation (Angstrom) along rotation vector (see DynDom info file)" },
 186         { "-head",     FALSE, etRVEC, {head},
 187           "First atom of the arrow vector" },
 188         { "-tail",     FALSE, etRVEC, {tail},
 189           "Last atom of the arrow vector" }
 190     };
 191     int               i, j, natoms, isize;
 192     t_trxstatus      *status;
 193     int              *index = NULL, *index_all;
 194     char             *grpname;
 195     real              angle, trans;
 196     rvec             *x, *v, *xout, *vout;
 197     matrix            box;
 198     gmx_output_env_t *oenv;
 199
 200     t_filenm          fnm[] = {
 201         { efPDB, "-f", "dyndom",  ffREAD },
 202         { efTRO, "-o", "rotated", ffWRITE },
 203         { efNDX, "-n", "domains", ffREAD }
 204     };
 205 #define NFILE asize(fnm)
 206
 207     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 208                            asize(desc), desc, 0, NULL, &oenv))
 209     {
 210         return 0;
 211     }
 212
 213     if (maxangle == 0)
 214     {
 215         gmx_fatal(FARGS, "maxangle not given");
 216     }
 217
 218     t_topology *top;
 219     snew(top, 1);
 220     read_tps_conf(opt2fn("-f", NFILE, fnm), top, NULL, &x, &v, box, FALSE);
 221     t_atoms  &atoms = top->atoms;
 222     if (atoms.pdbinfo == NULL)
 223     {
 224         snew(atoms.pdbinfo, atoms.nr);
 225     }
 226     natoms = atoms.nr;
 227     snew(xout, natoms);
 228     snew(vout, natoms);
 229
 230     printf("Select group to rotate:\n");
 231     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 232     printf("Going to rotate %s containg %d atoms\n", grpname, isize);
 233
 234     snew(index_all, atoms.nr);
 235     for (i = 0; (i < atoms.nr); i++)
 236     {
 237         index_all[i] = i;
 238     }
 239
 240     status = open_trx(opt2fn("-o", NFILE, fnm), "w");
 241
 242     label = 'A';
 243     for (i = 0; (i < nframes); i++, label++)
 244     {
 245         angle = angle0 + (i*(angle1-angle0))/(nframes-1);
 246         trans = trans0*0.1*angle/maxangle;
 247         printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
 248                i, label, angle, trans);
 249         rot_conf(&atoms, x, v, trans, angle, head, tail, isize, index, xout, vout);
 250
 251         if (label > 'Z')
 252         {
 253             label -= 26;
 254         }
 255         for (j = 0; (j < atoms.nr); j++)
 256         {
 257             atoms.resinfo[atoms.atom[j].resind].chainid = label;
 258         }
 259
 260         write_trx(status, atoms.nr, index_all, &atoms, i, angle, box, xout, vout, NULL);
 261     }
 262     close_trx(status);
 263
 264     return 0;
 265 }