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