src/tools/genpr.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 <math.h>
  37 #include "sysstuff.h"
  38 #include "statutil.h"
  39 #include "string.h"
  40 #include "copyrite.h"
  41 #include "smalloc.h"
  42 #include "typedefs.h"
  43 #include "confio.h"
  44 #include "futil.h"
  45 #include "macros.h"
  46 #include "index.h"
  47
  48 int main(int argc,char *argv[])
  49 {
  50   static char *desc[] = {
  51     "genpr produces an include file for a topology containing",
  52     "a list of atom numbers and three force constants for the",
  53     "X, Y and Z direction. A single isotropic force constant may",
  54     "be given on the command line instead of three components.[PAR]",
  55     "WARNING: genpr only works for the first molecule.",
  56     "Position restraints are interactions within molecules, therefore",
  57     "they should be included within the correct [TT][ moleculetype ][tt]",
  58     "block in the topology. Since the atom numbers in every moleculetype",
  59     "in the topology start at 1 and the numbers in the input file for",
  60     "genpr number consecutively from 1, genpr will only produce a useful",
  61     "file for the first molecule.[PAR]",
  62     "The -of option produces an index file that can be used for",
  63     "freezing atoms. In this case the input file must be a pdb file."
  64   };
  65   static rvec    fc={1000.0,1000.0,1000.0};
  66   static real    freeze_level;
  67   t_pargs pa[] = {
  68     { "-fc", FALSE, etRVEC, {fc},
  69       "force constants (kJ mol-1 nm-2)" },
  70     { "-freeze", FALSE, etREAL, {&freeze_level},
  71       "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" }
  72   };
  73
  74   t_atoms atoms;
  75   int     i;
  76   FILE    *out;
  77   int          igrp;
  78   atom_id      *ind_grp;
  79   char         *gn_grp;
  80   char         title[STRLEN];
  81   matrix       box;
  82   bool         bFreeze;
  83
  84   t_filenm fnm[] = {
  85     { efSTX, "-f",  NULL,    ffREAD },
  86     { efNDX, "-n",  NULL,    ffOPTRD },
  87     { efITP, "-o",  "posre", ffWRITE },
  88     { efNDX, "-of", "freeze",    ffOPTWR }
  89   };
  90 #define NFILE asize(fnm)
  91
  92   CopyRight(stderr,argv[0]);
  93   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
  94                     asize(desc),desc,0,NULL);
  95
  96   bFreeze = opt2bSet("-of",NFILE,fnm) || opt2parg_bSet("-freeze",asize(pa),pa);
  97
  98   if ( !opt2bSet("-n",NFILE,fnm) ) {
  99     if ( !ftp2bSet(efSTX,NFILE,fnm) )
 100       fatal_error(0,"no index file and no structure file suplied");
 101     else {
 102       rvec *x,*v;
 103
 104       get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&(atoms.nr));
 105       init_t_atoms(&atoms,atoms.nr,TRUE);
 106       snew(x,atoms.nr);
 107       snew(v,atoms.nr);
 108       fprintf(stderr,"\nReading structure file\n");
 109       read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,x,v,box);
 110       sfree(x);
 111       sfree(v);
 112     }
 113   }
 114   if (bFreeze) {
 115     if (atoms.pdbinfo == NULL)
 116       fatal_error(0,"No B-factors in input file %s, use a pdb file next time.",
 117                   ftp2fn(efSTX,NFILE,fnm));
 118
 119     out=opt2FILE("-of",NFILE,fnm,"w");
 120     fprintf(out,"[ freeze ]\n");
 121     for(i=0; (i<atoms.nr); i++) {
 122       if (atoms.pdbinfo[i].bfac <= freeze_level)
 123         fprintf(out,"%d\n",i+1);
 124     }
 125     fclose(out);
 126   }
 127   else {
 128     printf("Select group to position restrain\n");
 129     get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&igrp,&ind_grp,&gn_grp);
 130
 131     out=ftp2FILE(efITP,NFILE,fnm,"w");
 132     fprintf(out,"; position restraints for %s of %s\n\n",gn_grp,title);
 133     fprintf(out,"[ position_restraints ]\n");
 134     fprintf(out,";%3s %5s %9s %10s %10s\n","i","funct","fcx","fcy","fcz");
 135     for(i=0; i<igrp; i++)
 136       fprintf(out,"%4d %4d %10g %10g %10g\n",
 137               ind_grp[i]+1,1,fc[XX],fc[YY],fc[ZZ]);
 138     fclose(out);
 139   }
 140   thanx(stderr);
 141
 142   return 0;
 143 }