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