src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "statutil.h"
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include "confio.h"
  46 #include "copyrite.h"
  47 #include "xvgr.h"
  48 #include "index.h"
  49 #include "do_fit.h"
  50 #include "gmx_ana.h"
  51 #include "gmx_fatal.h"
  52
  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     for(j=0; (j<DIM); j++)
  65       xx[i][j] = fac0*x1[i][j] + fac1*x2[i][j];
  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   CopyRight(stderr,argv[0]);
 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     gmx_fatal(FARGS,"Number of atoms in first structure is %d, in second %d",
 131                 nat1,nat2);
 132
 133   init_t_atoms(&atoms,nat1,TRUE);
 134   snew(x1,nat1);
 135   snew(x2,nat1);
 136   snew(xx,nat1);
 137   snew(v,nat1);
 138
 139   read_stx_conf(opt2fn("-f1",NFILE,fnm),title,&atoms,x1,v,NULL,box);
 140   read_stx_conf(opt2fn("-f2",NFILE,fnm),title,&atoms,x2,v,NULL,box);
 141
 142   snew(mass,nat1);
 143   snew(index_all,nat1);
 144   for(i=0; (i<nat1); i++) {
 145     mass[i] = 1;
 146     index_all[i] = i;
 147   }
 148   if (bFit) {
 149     printf("Select group for LSQ superposition:\n");
 150     get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&is_lsq,&index_lsq,
 151               &grpname);
 152     reset_x(is_lsq,index_lsq,nat1,index_all,x1,mass);
 153     reset_x(is_lsq,index_lsq,nat1,index_all,x2,mass);
 154     do_fit(nat1,mass,x1,x2);
 155   }
 156
 157   bRMS = opt2bSet("-or",NFILE,fnm);
 158   if (bRMS) {
 159     fp = xvgropen(opt2fn("-or",NFILE,fnm),"RMSD","Conf","(nm)",oenv);
 160     xvgr_legend(fp,asize(leg),leg,oenv);
 161     printf("Select group for RMSD calculation:\n");
 162     get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&isize,&index,&grpname);
 163     printf("You selected group %s, containing %d atoms\n",grpname,isize);
 164     rms1 = rmsdev_ind(isize,index,mass,x1,x2);
 165     fprintf(stderr,"RMSD between input conformations is %g nm\n",rms1);
 166   }
 167
 168   snew(dummy,nat1);
 169   for(i=0; (i<nat1); i++)
 170     dummy[i] = i;
 171   status = open_trx(ftp2fn(efTRX,NFILE,fnm),"w");
 172
 173   for(i=0; (i<ninterm); i++) {
 174     fac = dointerp(nat1,x1,x2,xx,i,ninterm,first,last);
 175     write_trx(status,nat1,dummy,&atoms,i,fac,box,xx,NULL,NULL);
 176     if (bRMS) {
 177       rms1 = rmsdev_ind(isize,index,mass,x1,xx);
 178       rms2 = rmsdev_ind(isize,index,mass,x2,xx);
 179       fprintf(fp,"%10g  %10g  %10g\n",fac,rms1,rms2);
 180     }
 181   }
 182
 183   close_trx(status);
 184
 185   if (bRMS) {
 186     ffclose(fp);
 187     do_view(oenv,opt2fn("-or",NFILE,fnm),"-nxy");
 188   }
 189
 190   thanx(stderr);
 191
 192   return 0;
 193 }