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