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 "confio.h"
  43 #include "copyrite.h"
  44 #include "xvgr.h"
  45 #include "index.h"
  46 #include "do_fit.h"
  47 #include "gmx_ana.h"
  48 #include "gmx_fatal.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     parse_common_args(&argc, argv, PCA_CAN_VIEW,
 125                       NFILE, fnm, asize(pa), pa, asize(desc), desc,
 126                       0, NULL, &oenv);
 127     get_stx_coordnum (opt2fn("-f1", NFILE, fnm), &nat1);
 128     get_stx_coordnum (opt2fn("-f2", NFILE, fnm), &nat2);
 129     if (nat1 != nat2)
 130     {
 131         gmx_fatal(FARGS, "Number of atoms in first structure is %d, in second %d",
 132                   nat1, nat2);
 133     }
 134
 135     init_t_atoms(&atoms, nat1, TRUE);
 136     snew(x1, nat1);
 137     snew(x2, nat1);
 138     snew(xx, nat1);
 139     snew(v, nat1);
 140
 141     read_stx_conf(opt2fn("-f1", NFILE, fnm), title, &atoms, x1, v, NULL, box);
 142     read_stx_conf(opt2fn("-f2", NFILE, fnm), title, &atoms, x2, v, NULL, box);
 143
 144     snew(mass, nat1);
 145     snew(index_all, nat1);
 146     for (i = 0; (i < nat1); i++)
 147     {
 148         mass[i]      = 1;
 149         index_all[i] = i;
 150     }
 151     if (bFit)
 152     {
 153         printf("Select group for LSQ superposition:\n");
 154         get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &is_lsq, &index_lsq,
 155                   &grpname);
 156         reset_x(is_lsq, index_lsq, nat1, index_all, x1, mass);
 157         reset_x(is_lsq, index_lsq, nat1, index_all, x2, mass);
 158         do_fit(nat1, mass, x1, x2);
 159     }
 160
 161     bRMS = opt2bSet("-or", NFILE, fnm);
 162     if (bRMS)
 163     {
 164         fp = xvgropen(opt2fn("-or", NFILE, fnm), "RMSD", "Conf", "(nm)", oenv);
 165         xvgr_legend(fp, asize(leg), leg, oenv);
 166         printf("Select group for RMSD calculation:\n");
 167         get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &isize, &index, &grpname);
 168         printf("You selected group %s, containing %d atoms\n", grpname, isize);
 169         rms1 = rmsdev_ind(isize, index, mass, x1, x2);
 170         fprintf(stderr, "RMSD between input conformations is %g nm\n", rms1);
 171     }
 172
 173     snew(dummy, nat1);
 174     for (i = 0; (i < nat1); i++)
 175     {
 176         dummy[i] = i;
 177     }
 178     status = open_trx(ftp2fn(efTRX, NFILE, fnm), "w");
 179
 180     for (i = 0; (i < ninterm); i++)
 181     {
 182         fac = dointerp(nat1, x1, x2, xx, i, ninterm, first, last);
 183         write_trx(status, nat1, dummy, &atoms, i, fac, box, xx, NULL, NULL);
 184         if (bRMS)
 185         {
 186             rms1 = rmsdev_ind(isize, index, mass, x1, xx);
 187             rms2 = rmsdev_ind(isize, index, mass, x2, xx);
 188             fprintf(fp, "%10g  %10g  %10g\n", fac, rms1, rms2);
 189         }
 190     }
 191
 192     close_trx(status);
 193
 194     if (bRMS)
 195     {
 196         ffclose(fp);
 197         do_view(oenv, opt2fn("-or", NFILE, fnm), "-nxy");
 198     }
 199
 200     thanx(stderr);
 201
 202     return 0;
 203 }