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