src/gromacs/gmxana/gmx_dyndom.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "3dview.h"
  40 #include "statutil.h"
  41 #include "smalloc.h"
  42 #include "copyrite.h"
  43 #include "index.h"
  44 #include "confio.h"
  45 #include "gmx_fatal.h"
  46 #include "vec.h"
  47 #include "physics.h"
  48 #include "random.h"
  49 #include "gmx_ana.h"
  50 #include "macros.h"
  51
  52
  53 static void rot_conf(t_atoms *atoms, rvec x[], rvec v[], real trans, real angle,
  54                      rvec head, rvec tail, matrix box, int isize, atom_id index[],
  55                      rvec xout[], rvec vout[])
  56 {
  57     rvec     arrow, center, xcm;
  58     real     theta, phi, arrow_len;
  59     mat4     Rx, Ry, Rz, Rinvy, Rinvz, Mtot, Tcm, Tinvcm, Tx;
  60     mat4     temp1, temp2, temp3, temp4, temp21, temp43;
  61     vec4     xv;
  62     int      i, j, ai;
  63
  64     rvec_sub(tail, head, arrow);
  65     arrow_len = norm(arrow);
  66     if (debug)
  67     {
  68         fprintf(debug, "Arrow vector:   %10.4f  %10.4f  %10.4f\n",
  69                 arrow[XX], arrow[YY], arrow[ZZ]);
  70         fprintf(debug, "Effective translation %g nm\n", trans);
  71     }
  72     if (arrow_len == 0.0)
  73     {
  74         gmx_fatal(FARGS, "Arrow vector not given");
  75     }
  76
  77     /* Copy all aoms to output */
  78     for (i = 0; (i < atoms->nr); i++)
  79     {
  80         copy_rvec(x[i], xout[i]);
  81         copy_rvec(v[i], vout[i]);
  82     }
  83
  84     /* Compute center of mass and move atoms there */
  85     clear_rvec(xcm);
  86     for (i = 0; (i < isize); i++)
  87     {
  88         rvec_inc(xcm, x[index[i]]);
  89     }
  90     for (i = 0; (i < DIM); i++)
  91     {
  92         xcm[i] /= isize;
  93     }
  94     if (debug)
  95     {
  96         fprintf(debug, "Center of mass: %10.4f  %10.4f  %10.4f\n",
  97                 xcm[XX], xcm[YY], xcm[ZZ]);
  98     }
  99     for (i = 0; (i < isize); i++)
 100     {
 101         rvec_sub(x[index[i]], xcm, xout[index[i]]);
 102     }
 103
 104     /* Compute theta and phi that describe the arrow */
 105     theta = acos(arrow[ZZ]/arrow_len);
 106     phi   = atan2(arrow[YY]/arrow_len, arrow[XX]/arrow_len);
 107     if (debug)
 108     {
 109         fprintf(debug, "Phi = %.1f, Theta = %.1f\n", RAD2DEG*phi, RAD2DEG*theta);
 110     }
 111
 112     /* Now the total rotation matrix: */
 113     /* Rotate a couple of times */
 114     rotate(ZZ, -phi, Rz);
 115     rotate(YY, M_PI/2-theta, Ry);
 116     rotate(XX, angle*DEG2RAD, Rx);
 117     Rx[WW][XX] = trans;
 118     rotate(YY, theta-M_PI/2, Rinvy);
 119     rotate(ZZ, phi, Rinvz);
 120
 121     mult_matrix(temp1, Ry, Rz);
 122     mult_matrix(temp2, Rinvy, Rx);
 123     mult_matrix(temp3, temp2, temp1);
 124     mult_matrix(Mtot, Rinvz, temp3);
 125
 126     print_m4(debug, "Rz", Rz);
 127     print_m4(debug, "Ry", Ry);
 128     print_m4(debug, "Rx", Rx);
 129     print_m4(debug, "Rinvy", Rinvy);
 130     print_m4(debug, "Rinvz", Rinvz);
 131     print_m4(debug, "Mtot", Mtot);
 132
 133     for (i = 0; (i < isize); i++)
 134     {
 135         ai = index[i];
 136         m4_op(Mtot, xout[ai], xv);
 137         rvec_add(xv, xcm, xout[ai]);
 138         m4_op(Mtot, v[ai], xv);
 139         copy_rvec(xv, vout[ai]);
 140     }
 141 }
 142
 143 int gmx_dyndom(int argc, char *argv[])
 144 {
 145     const char  *desc[] = {
 146         "[TT]g_dyndom[tt] reads a [TT].pdb[tt] file output from DynDom",
 147         "(http://www.cmp.uea.ac.uk/dyndom/).",
 148         "It reads the coordinates, the coordinates of the rotation axis,",
 149         "and an index file containing the domains.",
 150         "Furthermore, it takes the first and last atom of the arrow file",
 151         "as command line arguments (head and tail) and",
 152         "finally it takes the translation vector (given in DynDom info file)",
 153         "and the angle of rotation (also as command line arguments). If the angle",
 154         "determined by DynDom is given, one should be able to recover the",
 155         "second structure used for generating the DynDom output.",
 156         "Because of limited numerical accuracy this should be verified by",
 157         "computing an all-atom RMSD (using [TT]g_confrms[tt]) rather than by file",
 158         "comparison (using diff).[PAR]",
 159         "The purpose of this program is to interpolate and extrapolate the",
 160         "rotation as found by DynDom. As a result unphysical structures with",
 161         "long or short bonds, or overlapping atoms may be produced. Visual",
 162         "inspection, and energy minimization may be necessary to",
 163         "validate the structure."
 164     };
 165     static real  trans0 = 0;
 166     static rvec  head   = { 0, 0, 0 };
 167     static rvec  tail   = { 0, 0, 0 };
 168     static real  angle0 = 0, angle1 = 0, maxangle = 0;
 169     static int   label  = 0, nframes = 11;
 170     t_pargs      pa[]   = {
 171         { "-firstangle",    FALSE, etREAL, {&angle0},
 172           "Angle of rotation about rotation vector" },
 173         { "-lastangle",    FALSE, etREAL, {&angle1},
 174           "Angle of rotation about rotation vector" },
 175         { "-nframe",   FALSE, etINT,  {&nframes},
 176           "Number of steps on the pathway" },
 177         { "-maxangle", FALSE, etREAL, {&maxangle},
 178           "DymDom dtermined angle of rotation about rotation vector" },
 179         { "-trans",    FALSE, etREAL, {&trans0},
 180           "Translation (Angstrom) along rotation vector (see DynDom info file)" },
 181         { "-head",     FALSE, etRVEC, {head},
 182           "First atom of the arrow vector" },
 183         { "-tail",     FALSE, etRVEC, {tail},
 184           "Last atom of the arrow vector" }
 185     };
 186     int          i, j, natoms, isize;
 187     t_trxstatus *status;
 188     atom_id     *index = NULL, *index_all;
 189     char         title[256], *grpname;
 190     t_atoms      atoms;
 191     real         angle, trans;
 192     rvec        *x, *v, *xout, *vout;
 193     matrix       box;
 194     output_env_t oenv;
 195
 196     t_filenm     fnm[] = {
 197         { efPDB, "-f", "dyndom",  ffREAD },
 198         { efTRO, "-o", "rotated", ffWRITE },
 199         { efNDX, "-n", "domains", ffREAD }
 200     };
 201 #define NFILE asize(fnm)
 202
 203     parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 204                       asize(desc), desc, 0, NULL, &oenv);
 205
 206     get_stx_coordnum (opt2fn("-f", NFILE, fnm), &natoms);
 207     init_t_atoms(&atoms, natoms, TRUE);
 208     snew(x, natoms);
 209     snew(v, natoms);
 210     read_stx_conf(opt2fn("-f", NFILE, fnm), title, &atoms, x, v, NULL, box);
 211     snew(xout, natoms);
 212     snew(vout, natoms);
 213
 214     printf("Select group to rotate:\n");
 215     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 216     printf("Going to rotate %s containg %d atoms\n", grpname, isize);
 217
 218     snew(index_all, atoms.nr);
 219     for (i = 0; (i < atoms.nr); i++)
 220     {
 221         index_all[i] = i;
 222     }
 223
 224     status = open_trx(opt2fn("-o", NFILE, fnm), "w");
 225
 226     label = 'A';
 227     for (i = 0; (i < nframes); i++, label++)
 228     {
 229         angle = angle0 + (i*(angle1-angle0))/(nframes-1);
 230         trans = trans0*0.1*angle/maxangle;
 231         printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
 232                i, label, angle, trans);
 233         rot_conf(&atoms, x, v, trans, angle, head, tail, box, isize, index, xout, vout);
 234
 235         if (label > 'Z')
 236         {
 237             label -= 26;
 238         }
 239         for (j = 0; (j < atoms.nr); j++)
 240         {
 241             atoms.resinfo[atoms.atom[j].resind].chainid = label;
 242         }
 243
 244         write_trx(status, atoms.nr, index_all, &atoms, i, angle, box, xout, vout, NULL);
 245     }
 246     close_trx(status);
 247
 248     thanx(stderr);
 249
 250     return 0;
 251 }