src/gromacs/gmxana/gmx_morph.cpp

   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,2015,2017, 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 "gromacs/commandline/pargs.h"
  40 #include "gromacs/commandline/viewit.h"
  41 #include "gromacs/fileio/confio.h"
  42 #include "gromacs/fileio/trxio.h"
  43 #include "gromacs/fileio/xvgr.h"
  44 #include "gromacs/gmxana/gmx_ana.h"
  45 #include "gromacs/math/do_fit.h"
  46 #include "gromacs/topology/atoms.h"
  47 #include "gromacs/topology/index.h"
  48 #include "gromacs/topology/topology.h"
  49 #include "gromacs/utility/arraysize.h"
  50 #include "gromacs/utility/cstringutil.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/smalloc.h"
  53
  54 static real dointerp(int n, rvec x1[], rvec x2[], rvec xx[],
  55                      int I, int N, real first, real last)
  56 {
  57     int    i, j;
  58     double fac, fac0, fac1;
  59
  60     fac  = first + (I*(last-first))/(N-1);
  61     fac0 = 1-fac;
  62     fac1 = fac;
  63     for (i = 0; (i < n); i++)
  64     {
  65         for (j = 0; (j < DIM); j++)
  66         {
  67             xx[i][j] = fac0*x1[i][j] + fac1*x2[i][j];
  68         }
  69     }
  70
  71     return fac;
  72 }
  73
  74 int gmx_morph(int argc, char *argv[])
  75 {
  76     const char       *desc[] = {
  77         "[THISMODULE] does a linear interpolation of conformations in order to",
  78         "create intermediates. Of course these are completely unphysical, but",
  79         "that you may try to justify yourself. Output is in the form of a ",
  80         "generic trajectory. The number of intermediates can be controlled with",
  81         "the [TT]-ninterm[tt] flag. The first and last flag correspond to the way of",
  82         "interpolating: 0 corresponds to input structure 1 while",
  83         "1 corresponds to input structure 2.",
  84         "If you specify [TT]-first[tt] < 0 or [TT]-last[tt] > 1 extrapolation will be",
  85         "on the path from input structure x[SUB]1[sub] to x[SUB]2[sub]. In general, the coordinates",
  86         "of the intermediate x(i) out of N total intermediates correspond to:[PAR]",
  87         "x(i) = x[SUB]1[sub] + (first+(i/(N-1))*(last-first))*(x[SUB]2[sub]-x[SUB]1[sub])[PAR]",
  88         "Finally the RMSD with respect to both input structures can be computed",
  89         "if explicitly selected ([TT]-or[tt] option). In that case, an index file may be",
  90         "read to select the group from which the RMS is computed."
  91     };
  92     t_filenm          fnm[] = {
  93         { efSTX, "-f1", "conf1",  ffREAD },
  94         { efSTX, "-f2", "conf2",  ffREAD },
  95         { efTRX, "-o",  "interm", ffWRITE },
  96         { efXVG, "-or", "rms-interm", ffOPTWR },
  97         { efNDX, "-n",  "index",  ffOPTRD }
  98     };
  99 #define NFILE asize(fnm)
 100     static  int       ninterm = 11;
 101     static  real      first   = 0.0;
 102     static  real      last    = 1.0;
 103     static  gmx_bool  bFit    = TRUE;
 104     t_pargs           pa []   = {
 105         { "-ninterm", FALSE, etINT,  {&ninterm},
 106           "Number of intermediates" },
 107         { "-first",   FALSE, etREAL, {&first},
 108           "Corresponds to first generated structure (0 is input x[SUB]1[sub], see above)" },
 109         { "-last",    FALSE, etREAL, {&last},
 110           "Corresponds to last generated structure (1 is input x[SUB]2[sub], see above)" },
 111         { "-fit",     FALSE, etBOOL, {&bFit},
 112           "Do a least squares fit of the second to the first structure before interpolating" }
 113     };
 114     const char       *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
 115     FILE             *fp    = nullptr;
 116     int               i, isize, is_lsq, nat1, nat2;
 117     t_trxstatus      *status;
 118     int              *index, *index_lsq, *index_all, *dummy;
 119     rvec             *x1, *x2, *xx;
 120     matrix            box;
 121     real              rms1, rms2, fac, *mass;
 122     char             *grpname;
 123     gmx_bool          bRMS;
 124     gmx_output_env_t *oenv;
 125
 126     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
 127                            NFILE, fnm, asize(pa), pa, asize(desc), desc,
 128                            0, nullptr, &oenv))
 129     {
 130         return 0;
 131     }
 132
 133     t_topology *top;
 134     snew(top, 1);
 135     read_tps_conf(opt2fn("-f1", NFILE, fnm), top, nullptr, &x1, nullptr, box, FALSE);
 136     nat1 = top->atoms.nr;
 137     read_tps_conf(opt2fn("-f2", NFILE, fnm), top, nullptr, &x2, nullptr, box, FALSE);
 138     nat2 = top->atoms.nr;
 139     if (nat1 != nat2)
 140     {
 141         gmx_fatal(FARGS, "Number of atoms in first structure is %d, in second %d",
 142                   nat1, nat2);
 143     }
 144     snew(xx, nat1);
 145     t_atoms  &atoms = top->atoms;
 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, nullptr, nullptr);
 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         xvgrclose(fp);
 200         do_view(oenv, opt2fn("-or", NFILE, fnm), "-nxy");
 201     }
 202
 203     return 0;
 204 }