src/gromacs/gmxana/gmx_dyndom.c

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