src/tools/gmx_genpr.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 <math.h>
  40 #include "sysstuff.h"
  41 #include "statutil.h"
  42 #include "string.h"
  43 #include "copyrite.h"
  44 #include "smalloc.h"
  45 #include "typedefs.h"
  46 #include "confio.h"
  47 #include "futil.h"
  48 #include "macros.h"
  49 #include "vec.h"
  50 #include "index.h"
  51 #include "gmx_fatal.h"
  52 #include "gmx_ana.h"
  53
  54 int gmx_genpr(int argc,char *argv[])
  55 {
  56   const char *desc[] = {
  57     "genrestr produces an include file for a topology containing",
  58     "a list of atom numbers and three force constants for the",
  59     "X, Y and Z direction. A single isotropic force constant may",
  60     "be given on the command line instead of three components.[PAR]",
  61     "WARNING: position restraints only work for the one molecule at a time.",
  62     "Position restraints are interactions within molecules, therefore",
  63     "they should be included within the correct [TT][ moleculetype ][tt]",
  64     "block in the topology. Since the atom numbers in every moleculetype",
  65     "in the topology start at 1 and the numbers in the input file for",
  66     "genpr number consecutively from 1, genpr will only produce a useful",
  67     "file for the first molecule.[PAR]",
  68     "The -of option produces an index file that can be used for",
  69     "freezing atoms. In this case the input file must be a pdb file.[PAR]",
  70     "With the [TT]-disre[tt] option half a matrix of distance restraints",
  71     "is generated instead of position restraints. With this matrix, that",
  72     "one typically would apply to C-alpha atoms in a protein, one can",
  73     "maintain the overall conformation of a protein without tieing it to",
  74     "a specific position (as with position restraints)."
  75   };
  76   static rvec    fc={1000.0,1000.0,1000.0};
  77   static real    freeze_level = 0.0;
  78   static real    disre_dist = 0.1;
  79   static real    disre_frac = 0.0;
  80   static real    disre_up2  = 1.0;
  81   static bool    bDisre=FALSE;
  82   static bool    bConstr=FALSE;
  83   static real    cutoff = -1.0;
  84
  85   t_pargs pa[] = {
  86     { "-fc", FALSE, etRVEC, {fc},
  87       "force constants (kJ mol-1 nm-2)" },
  88     { "-freeze", FALSE, etREAL, {&freeze_level},
  89       "if the -of option or this one is given an index file will be written containing atom numbers of all atoms that have a B-factor less than the level given here" },
  90     { "-disre", FALSE, etBOOL, {&bDisre},
  91       "Generate a distance restraint matrix for all the atoms in index" },
  92     { "-disre_dist", FALSE, etREAL, {&disre_dist},
  93       "Distance range around the actual distance for generating distance restraints" },
  94     { "-disre_frac", FALSE, etREAL, {&disre_frac},
  95       "Fraction of distance to be used as interval rather than a fixed distance. If the fraction of the distance that you specify here is less than the distance given in the previous option, that one is used instead." },
  96     { "-disre_up2", FALSE, etREAL, {&disre_up2},
  97       "Distance between upper bound for distance restraints, and the distance at which the force becomes constant (see manual)" },
  98     { "-cutoff", FALSE, etREAL, {&cutoff},
  99       "Only generate distance restraints for atoms pairs within cutoff (nm)" },
 100     { "-constr", FALSE, etBOOL, {&bConstr},
 101       "Generate a constraint matrix rather than distance restraints" }
 102   };
 103 #define npargs asize(pa)
 104
 105   output_env_t oenv;
 106   t_atoms      *atoms=NULL;
 107   int          i,j,k;
 108   FILE         *out;
 109   int          igrp;
 110   real         d,dd,lo,hi;
 111   atom_id      *ind_grp;
 112   const char   *xfn,*nfn;
 113   char         *gn_grp;
 114   char         title[STRLEN];
 115   matrix       box;
 116   bool         bFreeze;
 117   rvec         dx,*x=NULL,*v=NULL;
 118
 119   t_filenm fnm[] = {
 120     { efSTX, "-f",  NULL,    ffREAD },
 121     { efNDX, "-n",  NULL,    ffOPTRD },
 122     { efITP, "-o",  "posre", ffWRITE },
 123     { efNDX, "-of", "freeze",    ffOPTWR }
 124   };
 125 #define NFILE asize(fnm)
 126
 127   CopyRight(stderr,argv[0]);
 128   parse_common_args(&argc,argv,0,NFILE,fnm,npargs,pa,
 129                     asize(desc),desc,0,NULL,&oenv);
 130
 131   bFreeze = opt2bSet("-of",NFILE,fnm) || opt2parg_bSet("-freeze",asize(pa),pa);
 132   bDisre  = bDisre || opt2parg_bSet("-disre_dist",npargs,pa);
 133   xfn     = opt2fn_null("-f",NFILE,fnm);
 134   nfn     = opt2fn_null("-n",NFILE,fnm);
 135
 136   if (( nfn == NULL ) && ( xfn == NULL))
 137     gmx_fatal(FARGS,"no index file and no structure file suplied");
 138
 139   if ((disre_frac < 0) || (disre_frac >= 1))
 140     gmx_fatal(FARGS,"disre_frac should be between 0 and 1");
 141   if (disre_dist < 0)
 142     gmx_fatal(FARGS,"disre_dist should be >= 0");
 143
 144   if (xfn != NULL) {
 145     snew(atoms,1);
 146     get_stx_coordnum(xfn,&(atoms->nr));
 147     init_t_atoms(atoms,atoms->nr,TRUE);
 148     snew(x,atoms->nr);
 149     snew(v,atoms->nr);
 150     fprintf(stderr,"\nReading structure file\n");
 151     read_stx_conf(xfn,title,atoms,x,v,NULL,box);
 152   }
 153
 154   if (bFreeze) {
 155     if (atoms && atoms->pdbinfo)
 156       gmx_fatal(FARGS,"No B-factors in input file %s, use a pdb file next time.",
 157                 xfn);
 158
 159     out=opt2FILE("-of",NFILE,fnm,"w");
 160     fprintf(out,"[ freeze ]\n");
 161     for(i=0; (i<atoms->nr); i++) {
 162       if (atoms->pdbinfo[i].bfac <= freeze_level)
 163         fprintf(out,"%d\n",i+1);
 164     }
 165     ffclose(out);
 166   }
 167   else if ((bDisre || bConstr) && x) {
 168     printf("Select group to generate %s matrix from\n",
 169            bConstr ? "constraint" : "distance restraint");
 170     get_index(atoms,nfn,1,&igrp,&ind_grp,&gn_grp);
 171
 172     out=ftp2FILE(efITP,NFILE,fnm,"w");
 173     if (bConstr) {
 174       fprintf(out,"; constraints for %s of %s\n\n",gn_grp,title);
 175       fprintf(out,"[ constraints ]\n");
 176       fprintf(out,";%4s %5s %1s %10s\n","i","j","1","dist");
 177     }
 178     else {
 179       fprintf(out,"; distance restraints for %s of %s\n\n",gn_grp,title);
 180       fprintf(out,"[ distance_restraints ]\n");
 181       fprintf(out,";%4s %5s %1s %5s %10s %10s %10s %10s %10s\n","i","j","?",
 182               "label","funct","lo","up1","up2","weight");
 183     }
 184     for(i=k=0; i<igrp; i++)
 185       for(j=i+1; j<igrp; j++,k++) {
 186         rvec_sub(x[ind_grp[i]],x[ind_grp[j]],dx);
 187         d = norm(dx);
 188         if (bConstr)
 189           fprintf(out,"%5d %5d %1d %10g\n",ind_grp[i]+1,ind_grp[j]+1,1,d);
 190         else {
 191           if (cutoff < 0 || d < cutoff)
 192           {
 193             if (disre_frac > 0)
 194               dd = min(disre_dist,disre_frac*d);
 195             else
 196               dd = disre_dist;
 197             lo = max(0,d-dd);
 198             hi = d+dd;
 199             fprintf(out,"%5d %5d %1d %5d %10d %10g %10g %10g %10g\n",
 200                   ind_grp[i]+1,ind_grp[j]+1,1,k,1,
 201                   lo,hi,hi+1,1.0);
 202                 }
 203         }
 204       }
 205     ffclose(out);
 206   }
 207   else {
 208     printf("Select group to position restrain\n");
 209     get_index(atoms,nfn,1,&igrp,&ind_grp,&gn_grp);
 210
 211     out=ftp2FILE(efITP,NFILE,fnm,"w");
 212     fprintf(out,"; position restraints for %s of %s\n\n",gn_grp,title);
 213     fprintf(out,"[ position_restraints ]\n");
 214     fprintf(out,";%3s %5s %9s %10s %10s\n","i","funct","fcx","fcy","fcz");
 215     for(i=0; i<igrp; i++)
 216       fprintf(out,"%4d %4d %10g %10g %10g\n",
 217               ind_grp[i]+1,1,fc[XX],fc[YY],fc[ZZ]);
 218     ffclose(out);
 219   }
 220   if (xfn) {
 221     sfree(x);
 222     sfree(v);
 223   }
 224
 225   thanx(stderr);
 226
 227   return 0;
 228 }