src/tools/g_morph.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Green Red Orange Magenta Azure Cyan Skyblue
  35  */
  36 #include "confio.h"
  37 #include "statutil.h"
  38 #include "smalloc.h"
  39 #include "macros.h"
  40 #include "copyrite.h"
  41 #include "xvgr.h"
  42 #include "index.h"
  43 #include "do_fit.h"
  44
  45 static real dointerp(int n,rvec x1[],rvec x2[],rvec xx[],
  46                     int I,int N,real first,real last)
  47 {
  48   int    i,j;
  49   double fac,fac0,fac1;
  50
  51   fac  = first + (I*(last-first))/(N-1);
  52   fac0 = 1-fac;
  53   fac1 = fac;
  54   for(i=0; (i<n); i++)
  55     for(j=0; (j<DIM); j++)
  56       xx[i][j] = fac0*x1[i][j] + fac1*x2[i][j];
  57
  58   return fac;
  59 }
  60
  61 int gmx_morph(int argc,char *argv[])
  62 {
  63   char *desc[] = {
  64     "g_morph does a linear interpolation of conformations in order to",
  65     "create intermediates. Of course these are completely unphysical, but",
  66     "that you may try to justify yourself. Output is in the form of a ",
  67     "generic trajectory. The number of intermediates can be controlled with",
  68     "the -ninterm flag. The first and last flag correspond to the way of",
  69     "interpolating: 0 corresponds to input structure 1 while",
  70     "1 corresponds to input strucutre 2.",
  71     "If you specify first < 0 or last > 1 extrapolation will be",
  72     "on the path from input structure x1 to x2. In general the coordinates",
  73     "of the intermediate x(i) out of N total intermidates correspond to:[PAR]",
  74     "x(i) = x1 + (first+(i/(N-1))*(last-first))*(x2-x1)[PAR]",
  75     "Finally the RMSD with respect to both input structures can be computed",
  76     "if explicitly selected (-or option). In that case an index file may be",
  77     "read to select what group RMS is computed from."
  78   };
  79   t_filenm fnm[] = {
  80     { efSTX, "-f1", "conf1",  ffREAD },
  81     { efSTX, "-f2", "conf2",  ffREAD },
  82     { efTRX, "-o",  "interm", ffWRITE },
  83     { efXVG, "-or", "rms-interm", ffOPTWR },
  84     { efNDX, "-n",  "index",  ffOPTRD }
  85   };
  86 #define NFILE asize(fnm)
  87   static  int  ninterm = 11;
  88   static  real first   = 0.0;
  89   static  real last    = 1.0;
  90   static  bool bFit    = TRUE;
  91   t_pargs pa [] = {
  92     { "-ninterm", FALSE, etINT,  {&ninterm},
  93       "Number of intermediates" },
  94     { "-first",   FALSE, etREAL, {&first},
  95       "Corresponds to first generated structure (0 is input x0, see above)" },
  96     { "-last",    FALSE, etREAL, {&last},
  97       "Corresponds to last generated structure (1 is input x1, see above)" },
  98     { "-fit",     FALSE, etBOOL, {&bFit},
  99       "Do a least squares fit of the second to the first structure before interpolating" }
 100   };
 101   char     *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
 102   FILE     *fp=NULL;
 103   int      i,isize,is_lsq,status,nat1,nat2;
 104   atom_id  *index,*index_lsq,*index_all,*dummy;
 105   t_atoms  atoms;
 106   rvec     *x1,*x2,*xx,*v;
 107   matrix   box;
 108   real     rms1,rms2,fac,*mass;
 109   char     title[STRLEN],*grpname;
 110   bool     bRMS;
 111
 112   CopyRight(stderr,argv[0]);
 113   parse_common_args(&argc,argv,PCA_CAN_VIEW,
 114                     NFILE,fnm,asize(pa),pa,asize(desc),desc,
 115                     0,NULL);
 116   get_stx_coordnum (opt2fn("-f1",NFILE,fnm),&nat1);
 117   get_stx_coordnum (opt2fn("-f2",NFILE,fnm),&nat2);
 118   if (nat1 != nat2)
 119     fatal_error(0,"Number of atoms in first structure is %d, in second %d",
 120                 nat1,nat2);
 121
 122   init_t_atoms(&atoms,nat1,TRUE);
 123   snew(x1,nat1);
 124   snew(x2,nat1);
 125   snew(xx,nat1);
 126   snew(v,nat1);
 127
 128   read_stx_conf(opt2fn("-f1",NFILE,fnm),title,&atoms,x1,v,box);
 129   read_stx_conf(opt2fn("-f2",NFILE,fnm),title,&atoms,x2,v,box);
 130
 131   snew(mass,nat1);
 132   snew(index_all,nat1);
 133   for(i=0; (i<nat1); i++) {
 134     mass[i] = 1;
 135     index_all[i] = i;
 136   }
 137   if (bFit) {
 138     printf("Select group for LSQ superposition:\n");
 139     get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&is_lsq,&index_lsq,
 140               &grpname);
 141     reset_x(is_lsq,index_lsq,nat1,index_all,x1,mass);
 142     reset_x(is_lsq,index_lsq,nat1,index_all,x2,mass);
 143     do_fit(nat1,mass,x1,x2);
 144   }
 145
 146   bRMS = opt2bSet("-or",NFILE,fnm);
 147   if (bRMS) {
 148     fp = xvgropen(opt2fn("-or",NFILE,fnm),"RMSD","Conf","(nm)");
 149     xvgr_legend(fp,asize(leg),leg);
 150     printf("Select group for RMSD calculation:\n");
 151     get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&isize,&index,&grpname);
 152     printf("You selected group %s, containing %d atoms\n",grpname,isize);
 153     rms1 = rmsdev_ind(isize,index,mass,x1,x2);
 154     fprintf(stderr,"RMSD between input conformations is %g nm\n",rms1);
 155   }
 156
 157   snew(dummy,nat1);
 158   for(i=0; (i<nat1); i++)
 159     dummy[i] = i;
 160   status = open_trx(ftp2fn(efTRX,NFILE,fnm),"w");
 161
 162   for(i=0; (i<ninterm); i++) {
 163     fac = dointerp(nat1,x1,x2,xx,i,ninterm,first,last);
 164     write_trx(status,nat1,dummy,&atoms,i,fac,box,xx,NULL);
 165     if (bRMS) {
 166       rms1 = rmsdev_ind(isize,index,mass,x1,xx);
 167       rms2 = rmsdev_ind(isize,index,mass,x2,xx);
 168       fprintf(fp,"%10g  %10g  %10g\n",fac,rms1,rms2);
 169     }
 170   }
 171
 172   close_trx(status);
 173
 174   if (bRMS) {
 175     fclose(fp);
 176     do_view(opt2fn("-or",NFILE,fnm),"-nxy");
 177   }
 178
 179   thanx(stderr);
 180
 181   return 0;
 182 }