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