src/gromacs/gmxana/gmx_morph.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 "statutil.h"
  40 #include "smalloc.h"
  41 #include "macros.h"
  42 #include "gromacs/fileio/confio.h"
  43 #include "xvgr.h"
  44 #include "index.h"
  45 #include "do_fit.h"
  46 #include "gmx_ana.h"
  47 #include "gmx_fatal.h"
  48 #include "gromacs/fileio/trxio.h"
  49
  50
  51 static real dointerp(int n, rvec x1[], rvec x2[], rvec xx[],
  52                      int I, int N, real first, real last)
  53 {
  54     int    i, j;
  55     double fac, fac0, fac1;
  56
  57     fac  = first + (I*(last-first))/(N-1);
  58     fac0 = 1-fac;
  59     fac1 = fac;
  60     for (i = 0; (i < n); i++)
  61     {
  62         for (j = 0; (j < DIM); j++)
  63         {
  64             xx[i][j] = fac0*x1[i][j] + fac1*x2[i][j];
  65         }
  66     }
  67
  68     return fac;
  69 }
  70
  71 int gmx_morph(int argc, char *argv[])
  72 {
  73     const char      *desc[] = {
  74         "[TT]g_morph[tt] does a linear interpolation of conformations in order to",
  75         "create intermediates. Of course these are completely unphysical, but",
  76         "that you may try to justify yourself. Output is in the form of a ",
  77         "generic trajectory. The number of intermediates can be controlled with",
  78         "the [TT]-ninterm[tt] flag. The first and last flag correspond to the way of",
  79         "interpolating: 0 corresponds to input structure 1 while",
  80         "1 corresponds to input structure 2.",
  81         "If you specify [TT]-first[tt] < 0 or [TT]-last[tt] > 1 extrapolation will be",
  82         "on the path from input structure x[SUB]1[sub] to x[SUB]2[sub]. In general, the coordinates",
  83         "of the intermediate x(i) out of N total intermediates correspond to:[PAR]",
  84         "x(i) = x[SUB]1[sub] + (first+(i/(N-1))*(last-first))*(x[SUB]2[sub]-x[SUB]1[sub])[PAR]",
  85         "Finally the RMSD with respect to both input structures can be computed",
  86         "if explicitly selected ([TT]-or[tt] option). In that case, an index file may be",
  87         "read to select the group from which the RMS is computed."
  88     };
  89     t_filenm         fnm[] = {
  90         { efSTX, "-f1", "conf1",  ffREAD },
  91         { efSTX, "-f2", "conf2",  ffREAD },
  92         { efTRX, "-o",  "interm", ffWRITE },
  93         { efXVG, "-or", "rms-interm", ffOPTWR },
  94         { efNDX, "-n",  "index",  ffOPTRD }
  95     };
  96 #define NFILE asize(fnm)
  97     static  int      ninterm = 11;
  98     static  real     first   = 0.0;
  99     static  real     last    = 1.0;
 100     static  gmx_bool bFit    = TRUE;
 101     t_pargs          pa []   = {
 102         { "-ninterm", FALSE, etINT,  {&ninterm},
 103           "Number of intermediates" },
 104         { "-first",   FALSE, etREAL, {&first},
 105           "Corresponds to first generated structure (0 is input x[SUB]1[sub], see above)" },
 106         { "-last",    FALSE, etREAL, {&last},
 107           "Corresponds to last generated structure (1 is input x[SUB]2[sub], see above)" },
 108         { "-fit",     FALSE, etBOOL, {&bFit},
 109           "Do a least squares fit of the second to the first structure before interpolating" }
 110     };
 111     const char      *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
 112     FILE            *fp    = NULL;
 113     int              i, isize, is_lsq, nat1, nat2;
 114     t_trxstatus     *status;
 115     atom_id         *index, *index_lsq, *index_all, *dummy;
 116     t_atoms          atoms;
 117     rvec            *x1, *x2, *xx, *v;
 118     matrix           box;
 119     real             rms1, rms2, fac, *mass;
 120     char             title[STRLEN], *grpname;
 121     gmx_bool         bRMS;
 122     output_env_t     oenv;
 123
 124     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
 125                            NFILE, fnm, asize(pa), pa, asize(desc), desc,
 126                            0, NULL, &oenv))
 127     {
 128         return 0;
 129     }
 130     get_stx_coordnum (opt2fn("-f1", NFILE, fnm), &nat1);
 131     get_stx_coordnum (opt2fn("-f2", NFILE, fnm), &nat2);
 132     if (nat1 != nat2)
 133     {
 134         gmx_fatal(FARGS, "Number of atoms in first structure is %d, in second %d",
 135                   nat1, nat2);
 136     }
 137
 138     init_t_atoms(&atoms, nat1, TRUE);
 139     snew(x1, nat1);
 140     snew(x2, nat1);
 141     snew(xx, nat1);
 142     snew(v, nat1);
 143
 144     read_stx_conf(opt2fn("-f1", NFILE, fnm), title, &atoms, x1, v, NULL, box);
 145     read_stx_conf(opt2fn("-f2", NFILE, fnm), title, &atoms, x2, v, NULL, box);
 146
 147     snew(mass, nat1);
 148     snew(index_all, nat1);
 149     for (i = 0; (i < nat1); i++)
 150     {
 151         mass[i]      = 1;
 152         index_all[i] = i;
 153     }
 154     if (bFit)
 155     {
 156         printf("Select group for LSQ superposition:\n");
 157         get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &is_lsq, &index_lsq,
 158                   &grpname);
 159         reset_x(is_lsq, index_lsq, nat1, index_all, x1, mass);
 160         reset_x(is_lsq, index_lsq, nat1, index_all, x2, mass);
 161         do_fit(nat1, mass, x1, x2);
 162     }
 163
 164     bRMS = opt2bSet("-or", NFILE, fnm);
 165     if (bRMS)
 166     {
 167         fp = xvgropen(opt2fn("-or", NFILE, fnm), "RMSD", "Conf", "(nm)", oenv);
 168         xvgr_legend(fp, asize(leg), leg, oenv);
 169         printf("Select group for RMSD calculation:\n");
 170         get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &isize, &index, &grpname);
 171         printf("You selected group %s, containing %d atoms\n", grpname, isize);
 172         rms1 = rmsdev_ind(isize, index, mass, x1, x2);
 173         fprintf(stderr, "RMSD between input conformations is %g nm\n", rms1);
 174     }
 175
 176     snew(dummy, nat1);
 177     for (i = 0; (i < nat1); i++)
 178     {
 179         dummy[i] = i;
 180     }
 181     status = open_trx(ftp2fn(efTRX, NFILE, fnm), "w");
 182
 183     for (i = 0; (i < ninterm); i++)
 184     {
 185         fac = dointerp(nat1, x1, x2, xx, i, ninterm, first, last);
 186         write_trx(status, nat1, dummy, &atoms, i, fac, box, xx, NULL, NULL);
 187         if (bRMS)
 188         {
 189             rms1 = rmsdev_ind(isize, index, mass, x1, xx);
 190             rms2 = rmsdev_ind(isize, index, mass, x2, xx);
 191             fprintf(fp, "%10g  %10g  %10g\n", fac, rms1, rms2);
 192         }
 193     }
 194
 195     close_trx(status);
 196
 197     if (bRMS)
 198     {
 199         ffclose(fp);
 200         do_view(oenv, opt2fn("-or", NFILE, fnm), "-nxy");
 201     }
 202
 203     return 0;
 204 }