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  gmx_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   const char *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
 108   FILE     *fp=NULL;
 109   int      i,isize,is_lsq,nat1,nat2;
 110   t_trxstatus *status;
 111   atom_id  *index,*index_lsq,*index_all,*dummy;
 112   t_atoms  atoms;
 113   rvec     *x1,*x2,*xx,*v;
 114   matrix   box;
 115   real     rms1,rms2,fac,*mass;
 116   char     title[STRLEN],*grpname;
 117   gmx_bool     bRMS;
 118   output_env_t oenv;
 119
 120   CopyRight(stderr,argv[0]);
 121   parse_common_args(&argc,argv,PCA_CAN_VIEW,
 122                     NFILE,fnm,asize(pa),pa,asize(desc),desc,
 123                     0,NULL,&oenv);
 124   get_stx_coordnum (opt2fn("-f1",NFILE,fnm),&nat1);
 125   get_stx_coordnum (opt2fn("-f2",NFILE,fnm),&nat2);
 126   if (nat1 != nat2)
 127     gmx_fatal(FARGS,"Number of atoms in first structure is %d, in second %d",
 128                 nat1,nat2);
 129
 130   init_t_atoms(&atoms,nat1,TRUE);
 131   snew(x1,nat1);
 132   snew(x2,nat1);
 133   snew(xx,nat1);
 134   snew(v,nat1);
 135
 136   read_stx_conf(opt2fn("-f1",NFILE,fnm),title,&atoms,x1,v,NULL,box);
 137   read_stx_conf(opt2fn("-f2",NFILE,fnm),title,&atoms,x2,v,NULL,box);
 138
 139   snew(mass,nat1);
 140   snew(index_all,nat1);
 141   for(i=0; (i<nat1); i++) {
 142     mass[i] = 1;
 143     index_all[i] = i;
 144   }
 145   if (bFit) {
 146     printf("Select group for LSQ superposition:\n");
 147     get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&is_lsq,&index_lsq,
 148               &grpname);
 149     reset_x(is_lsq,index_lsq,nat1,index_all,x1,mass);
 150     reset_x(is_lsq,index_lsq,nat1,index_all,x2,mass);
 151     do_fit(nat1,mass,x1,x2);
 152   }
 153
 154   bRMS = opt2bSet("-or",NFILE,fnm);
 155   if (bRMS) {
 156     fp = xvgropen(opt2fn("-or",NFILE,fnm),"RMSD","Conf","(nm)",oenv);
 157     xvgr_legend(fp,asize(leg),leg,oenv);
 158     printf("Select group for RMSD calculation:\n");
 159     get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&isize,&index,&grpname);
 160     printf("You selected group %s, containing %d atoms\n",grpname,isize);
 161     rms1 = rmsdev_ind(isize,index,mass,x1,x2);
 162     fprintf(stderr,"RMSD between input conformations is %g nm\n",rms1);
 163   }
 164
 165   snew(dummy,nat1);
 166   for(i=0; (i<nat1); i++)
 167     dummy[i] = i;
 168   status = open_trx(ftp2fn(efTRX,NFILE,fnm),"w");
 169
 170   for(i=0; (i<ninterm); i++) {
 171     fac = dointerp(nat1,x1,x2,xx,i,ninterm,first,last);
 172     write_trx(status,nat1,dummy,&atoms,i,fac,box,xx,NULL,NULL);
 173     if (bRMS) {
 174       rms1 = rmsdev_ind(isize,index,mass,x1,xx);
 175       rms2 = rmsdev_ind(isize,index,mass,x2,xx);
 176       fprintf(fp,"%10g  %10g  %10g\n",fac,rms1,rms2);
 177     }
 178   }
 179
 180   close_trx(status);
 181
 182   if (bRMS) {
 183     ffclose(fp);
 184     do_view(oenv,opt2fn("-or",NFILE,fnm),"-nxy");
 185   }
 186
 187   thanx(stderr);
 188
 189   return 0;
 190 }