src/gromacs/gmxpreprocess/h_db.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  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <string.h>
  41 #include "string2.h"
  42 #include "sysstuff.h"
  43 #include "smalloc.h"
  44 #include "futil.h"
  45 #include "symtab.h"
  46 #include "h_db.h"
  47 #include "gmxfio.h"
  48 #include "fflibutil.h"
  49 #include "gmx_fatal.h"
  50 #include "macros.h"
  51
  52 /* There are 11 types of adding hydrogens, numbered from
  53  * 1 thru 11. Each of these has a specific number of
  54  * control atoms, that determine how the hydrogens are added.
  55  * Here these number are given. Because arrays start at 0 an
  56  * extra dummy for index 0 is added
  57  */
  58 const int ncontrol[] = { -1, 3, 3, 3, 3, 4, 3, 1, 3, 3, 1, 1 };
  59 #define maxcontrol asize(ncontrol)
  60
  61 int compaddh(const void *a, const void *b)
  62 {
  63     t_hackblock *ah, *bh;
  64
  65     ah = (t_hackblock *)a;
  66     bh = (t_hackblock *)b;
  67     return gmx_strcasecmp(ah->name, bh->name);
  68 }
  69
  70 void print_ab(FILE *out, t_hack *hack, char *nname)
  71 {
  72     int i;
  73
  74     fprintf(out, "%d\t%d\t%s", hack->nr, hack->tp, nname);
  75     for (i = 0; (i < hack->nctl); i++)
  76     {
  77         fprintf(out, "\t%s", hack->a[i]);
  78     }
  79     fprintf(out, "\n");
  80 }
  81
  82
  83 void read_ab(char *line, const char *fn, t_hack *hack)
  84 {
  85     int  i, nh, tp, ns;
  86     char a[4][12];
  87     char hn[32];
  88
  89     ns = sscanf(line, "%d%d%s%s%s%s%s", &nh, &tp, hn, a[0], a[1], a[2], a[3]);
  90     if (ns < 4)
  91     {
  92         gmx_fatal(FARGS, "wrong format in input file %s on line\n%s\n", fn, line);
  93     }
  94
  95     hack->nr = nh;
  96     hack->tp = tp;
  97     if ((tp < 1) || (tp >= maxcontrol))
  98     {
  99         gmx_fatal(FARGS, "Error in hdb file %s:\nH-type should be in 1-%d. Offending line:\n%s", fn, maxcontrol-1, line);
 100     }
 101
 102     hack->nctl = ns - 3;
 103     if ((hack->nctl != ncontrol[hack->tp]) && (ncontrol[hack->tp] != -1))
 104     {
 105         gmx_fatal(FARGS, "Error in hdb file %s:\nWrong number of control atoms (%d iso %d) on line:\n%s\n", fn, hack->nctl, ncontrol[hack->tp], line);
 106     }
 107     for (i = 0; (i < hack->nctl); i++)
 108     {
 109         hack->a[i] = strdup(a[i]);
 110     }
 111     for (; i < 4; i++)
 112     {
 113         hack->a[i] = NULL;
 114     }
 115     hack->oname = NULL;
 116     hack->nname = strdup(hn);
 117     hack->atom  = NULL;
 118     hack->cgnr  = NOTSET;
 119     hack->bXSet = FALSE;
 120     for (i = 0; i < DIM; i++)
 121     {
 122         hack->newx[i] = NOTSET;
 123     }
 124 }
 125
 126 static void read_h_db_file(const char *hfn, int *nahptr, t_hackblock **ah)
 127 {
 128     FILE        *in;
 129     char         filebase[STRLEN], line[STRLEN], buf[STRLEN];
 130     int          i, n, nab, nah;
 131     t_hackblock *aah;
 132
 133     if (debug)
 134     {
 135         fprintf(debug, "Hydrogen Database (%s):\n", hfn);
 136     }
 137
 138     fflib_filename_base(hfn, filebase, STRLEN);
 139     /* Currently filebase is read and set, but not used.
 140      * hdb entries from any hdb file and be applied to rtp entries
 141      * in any rtp file.
 142      */
 143
 144     in = fflib_open(hfn);
 145
 146     nah = *nahptr;
 147     aah = *ah;
 148     while (fgets2(line, STRLEN-1, in))
 149     {
 150         if (sscanf(line, "%s%n", buf, &n) != 1)
 151         {
 152             fprintf(stderr, "Error in hdb file: nah = %d\nline = '%s'\n",
 153                     nah, line);
 154             break;
 155         }
 156         if (debug)
 157         {
 158             fprintf(debug, "%s", buf);
 159         }
 160         srenew(aah, nah+1);
 161         clear_t_hackblock(&aah[nah]);
 162         aah[nah].name     = strdup(buf);
 163         aah[nah].filebase = strdup(filebase);
 164
 165         if (sscanf(line+n, "%d", &nab) == 1)
 166         {
 167             if (debug)
 168             {
 169                 fprintf(debug, "  %d\n", nab);
 170             }
 171             snew(aah[nah].hack, nab);
 172             aah[nah].nhack = nab;
 173             for (i = 0; (i < nab); i++)
 174             {
 175                 if (feof(in))
 176                 {
 177                     gmx_fatal(FARGS, "Expected %d lines of hydrogens, found only %d "
 178                               "while reading Hydrogen Database %s residue %s",
 179                               nab, i-1, aah[nah].name, hfn);
 180                 }
 181                 if (NULL == fgets(buf, STRLEN, in))
 182                 {
 183                     gmx_fatal(FARGS, "Error reading from file %s", hfn);
 184                 }
 185                 read_ab(buf, hfn, &(aah[nah].hack[i]));
 186             }
 187         }
 188         nah++;
 189     }
 190     ffclose(in);
 191
 192     /* Sort the list (necessary to be able to use bsearch */
 193     qsort(aah, nah, (size_t)sizeof(**ah), compaddh);
 194
 195     /*
 196        if (debug)
 197        dump_h_db(hfn,nah,aah);
 198      */
 199
 200     *nahptr = nah;
 201     *ah     = aah;
 202 }
 203
 204 int read_h_db(const char *ffdir, t_hackblock **ah)
 205 {
 206     int    nhdbf, f;
 207     char **hdbf;
 208     int    nah;
 209     FILE  *fp;
 210
 211     /* Read the hydrogen database file(s).
 212      * Do not generate an error when no files are found.
 213      */
 214     nhdbf = fflib_search_file_end(ffdir, ".hdb", FALSE, &hdbf);
 215     nah   = 0;
 216     *ah   = NULL;
 217     for (f = 0; f < nhdbf; f++)
 218     {
 219         read_h_db_file(hdbf[f], &nah, ah);
 220         sfree(hdbf[f]);
 221     }
 222     sfree(hdbf);
 223
 224     return nah;
 225 }
 226
 227 t_hackblock *search_h_db(int nh, t_hackblock ah[], char *key)
 228 {
 229     t_hackblock ahkey, *result;
 230
 231     if (nh <= 0)
 232     {
 233         return NULL;
 234     }
 235
 236     ahkey.name = key;
 237
 238     result = (t_hackblock *)bsearch(&ahkey, ah, nh, (size_t)sizeof(ah[0]), compaddh);
 239
 240     return result;
 241 }