src/gromacs/gmxana/gmx_morph.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include "gromacs/legacyheaders/macros.h"
  42 #include "gromacs/fileio/confio.h"
  43 #include "gromacs/legacyheaders/viewit.h"
  44 #include "gromacs/topology/index.h"
  45 #include "gmx_ana.h"
  46 #include "gromacs/fileio/trxio.h"
  47
  48 #include "gromacs/commandline/pargs.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/math/do_fit.h"
  51 #include "gromacs/topology/atoms.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 static real dointerp(int n, rvec x1[], rvec x2[], rvec xx[],
  56                      int I, int N, real first, real last)
  57 {
  58     int    i, j;
  59     double fac, fac0, fac1;
  60
  61     fac  = first + (I*(last-first))/(N-1);
  62     fac0 = 1-fac;
  63     fac1 = fac;
  64     for (i = 0; (i < n); i++)
  65     {
  66         for (j = 0; (j < DIM); j++)
  67         {
  68             xx[i][j] = fac0*x1[i][j] + fac1*x2[i][j];
  69         }
  70     }
  71
  72     return fac;
  73 }
  74
  75 int gmx_morph(int argc, char *argv[])
  76 {
  77     const char      *desc[] = {
  78         "[THISMODULE] does a linear interpolation of conformations in order to",
  79         "create intermediates. Of course these are completely unphysical, but",
  80         "that you may try to justify yourself. Output is in the form of a ",
  81         "generic trajectory. The number of intermediates can be controlled with",
  82         "the [TT]-ninterm[tt] flag. The first and last flag correspond to the way of",
  83         "interpolating: 0 corresponds to input structure 1 while",
  84         "1 corresponds to input structure 2.",
  85         "If you specify [TT]-first[tt] < 0 or [TT]-last[tt] > 1 extrapolation will be",
  86         "on the path from input structure x[SUB]1[sub] to x[SUB]2[sub]. In general, the coordinates",
  87         "of the intermediate x(i) out of N total intermediates correspond to:[PAR]",
  88         "x(i) = x[SUB]1[sub] + (first+(i/(N-1))*(last-first))*(x[SUB]2[sub]-x[SUB]1[sub])[PAR]",
  89         "Finally the RMSD with respect to both input structures can be computed",
  90         "if explicitly selected ([TT]-or[tt] option). In that case, an index file may be",
  91         "read to select the group from which the RMS is computed."
  92     };
  93     t_filenm         fnm[] = {
  94         { efSTX, "-f1", "conf1",  ffREAD },
  95         { efSTX, "-f2", "conf2",  ffREAD },
  96         { efTRX, "-o",  "interm", ffWRITE },
  97         { efXVG, "-or", "rms-interm", ffOPTWR },
  98         { efNDX, "-n",  "index",  ffOPTRD }
  99     };
 100 #define NFILE asize(fnm)
 101     static  int      ninterm = 11;
 102     static  real     first   = 0.0;
 103     static  real     last    = 1.0;
 104     static  gmx_bool bFit    = TRUE;
 105     t_pargs          pa []   = {
 106         { "-ninterm", FALSE, etINT,  {&ninterm},
 107           "Number of intermediates" },
 108         { "-first",   FALSE, etREAL, {&first},
 109           "Corresponds to first generated structure (0 is input x[SUB]1[sub], see above)" },
 110         { "-last",    FALSE, etREAL, {&last},
 111           "Corresponds to last generated structure (1 is input x[SUB]2[sub], see above)" },
 112         { "-fit",     FALSE, etBOOL, {&bFit},
 113           "Do a least squares fit of the second to the first structure before interpolating" }
 114     };
 115     const char      *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
 116     FILE            *fp    = NULL;
 117     int              i, isize, is_lsq, nat1, nat2;
 118     t_trxstatus     *status;
 119     atom_id         *index, *index_lsq, *index_all, *dummy;
 120     t_atoms          atoms;
 121     rvec            *x1, *x2, *xx, *v;
 122     matrix           box;
 123     real             rms1, rms2, fac, *mass;
 124     char             title[STRLEN], *grpname;
 125     gmx_bool         bRMS;
 126     output_env_t     oenv;
 127
 128     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
 129                            NFILE, fnm, asize(pa), pa, asize(desc), desc,
 130                            0, NULL, &oenv))
 131     {
 132         return 0;
 133     }
 134     get_stx_coordnum (opt2fn("-f1", NFILE, fnm), &nat1);
 135     get_stx_coordnum (opt2fn("-f2", NFILE, fnm), &nat2);
 136     if (nat1 != nat2)
 137     {
 138         gmx_fatal(FARGS, "Number of atoms in first structure is %d, in second %d",
 139                   nat1, nat2);
 140     }
 141
 142     init_t_atoms(&atoms, nat1, TRUE);
 143     snew(x1, nat1);
 144     snew(x2, nat1);
 145     snew(xx, nat1);
 146     snew(v, nat1);
 147
 148     read_stx_conf(opt2fn("-f1", NFILE, fnm), title, &atoms, x1, v, NULL, box);
 149     read_stx_conf(opt2fn("-f2", NFILE, fnm), title, &atoms, x2, v, NULL, box);
 150
 151     snew(mass, nat1);
 152     snew(index_all, nat1);
 153     for (i = 0; (i < nat1); i++)
 154     {
 155         mass[i]      = 1;
 156         index_all[i] = i;
 157     }
 158     if (bFit)
 159     {
 160         printf("Select group for LSQ superposition:\n");
 161         get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &is_lsq, &index_lsq,
 162                   &grpname);
 163         reset_x(is_lsq, index_lsq, nat1, index_all, x1, mass);
 164         reset_x(is_lsq, index_lsq, nat1, index_all, x2, mass);
 165         do_fit(nat1, mass, x1, x2);
 166     }
 167
 168     bRMS = opt2bSet("-or", NFILE, fnm);
 169     if (bRMS)
 170     {
 171         fp = xvgropen(opt2fn("-or", NFILE, fnm), "RMSD", "Conf", "(nm)", oenv);
 172         xvgr_legend(fp, asize(leg), leg, oenv);
 173         printf("Select group for RMSD calculation:\n");
 174         get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &isize, &index, &grpname);
 175         printf("You selected group %s, containing %d atoms\n", grpname, isize);
 176         rms1 = rmsdev_ind(isize, index, mass, x1, x2);
 177         fprintf(stderr, "RMSD between input conformations is %g nm\n", rms1);
 178     }
 179
 180     snew(dummy, nat1);
 181     for (i = 0; (i < nat1); i++)
 182     {
 183         dummy[i] = i;
 184     }
 185     status = open_trx(ftp2fn(efTRX, NFILE, fnm), "w");
 186
 187     for (i = 0; (i < ninterm); i++)
 188     {
 189         fac = dointerp(nat1, x1, x2, xx, i, ninterm, first, last);
 190         write_trx(status, nat1, dummy, &atoms, i, fac, box, xx, NULL, NULL);
 191         if (bRMS)
 192         {
 193             rms1 = rmsdev_ind(isize, index, mass, x1, xx);
 194             rms2 = rmsdev_ind(isize, index, mass, x2, xx);
 195             fprintf(fp, "%10g  %10g  %10g\n", fac, rms1, rms2);
 196         }
 197     }
 198
 199     close_trx(status);
 200
 201     if (bRMS)
 202     {
 203         gmx_ffclose(fp);
 204         do_view(oenv, opt2fn("-or", NFILE, fnm), "-nxy");
 205     }
 206
 207     return 0;
 208 }