src/gromacs/gmxlib/atomprop.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  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "sysstuff.h"
  41 #include "smalloc.h"
  42 #include "string2.h"
  43 #include "futil.h"
  44 #include "maths.h"
  45 #include "gmx_fatal.h"
  46 #include "atomprop.h"
  47 #include "maths.h"
  48 #include "macros.h"
  49 #include "index.h"
  50 #include "strdb.h"
  51 #include "copyrite.h"
  52
  53 typedef struct {
  54   gmx_bool   bSet;
  55   int    nprop,maxprop;
  56   char   *db;
  57   double def;
  58   char   **atomnm;
  59   char   **resnm;
  60   gmx_bool   *bAvail;
  61   real   *value;
  62 } aprop_t;
  63
  64 typedef struct gmx_atomprop {
  65   gmx_bool       bWarned;
  66   aprop_t    prop[epropNR];
  67   gmx_residuetype_t restype;
  68 } gmx_atomprop;
  69
  70
  71
  72 /* NOTFOUND should be smallest, others larger in increasing priority */
  73 enum { NOTFOUND=-4, WILDCARD, WILDPROT, PROTEIN };
  74
  75 /* return number of matching characters,
  76    or NOTFOUND if not at least all characters in char *database match */
  77 static int dbcmp_len(char *search, char *database)
  78 {
  79   int i;
  80
  81   i=0;
  82   while(search[i] && database[i] && (search[i]==database[i]) )
  83     i++;
  84
  85   if (database[i])
  86     i=NOTFOUND;
  87   return i;
  88 }
  89
  90 static int get_prop_index(aprop_t *ap,gmx_residuetype_t restype,
  91                           char *resnm,char *atomnm,
  92                           gmx_bool *bExact)
  93 {
  94   int  i,j=NOTFOUND;
  95   long int alen,rlen;
  96   long int malen,mrlen;
  97   gmx_bool bProtein,bProtWild;
  98
  99   bProtein  = gmx_residuetype_is_protein(restype,resnm);
 100   bProtWild = (strcmp(resnm,"AAA")==0);
 101   malen = NOTFOUND;
 102   mrlen = NOTFOUND;
 103   for(i=0; (i<ap->nprop); i++) {
 104     rlen = dbcmp_len(resnm, ap->resnm[i]);
 105     if (rlen == NOTFOUND) {
 106       if ( (strcmp(ap->resnm[i],"*")==0) ||
 107            (strcmp(ap->resnm[i],"???")==0) )
 108         rlen=WILDCARD;
 109       else if (strcmp(ap->resnm[i],"AAA")==0)
 110         rlen=WILDPROT;
 111     }
 112     alen = dbcmp_len(atomnm, ap->atomnm[i]);
 113     if ( (alen > NOTFOUND) && (rlen > NOTFOUND)) {
 114       if ( ( (alen > malen) && (rlen >= mrlen)) ||
 115            ( (rlen > mrlen) && (alen >= malen) ) ) {
 116         malen = alen;
 117         mrlen = rlen;
 118         j     = i;
 119       }
 120     }
 121   }
 122
 123   *bExact = ((malen == (long int)strlen(atomnm)) &&
 124              ((mrlen == (long int)strlen(resnm)) ||
 125               ((mrlen == WILDPROT) && bProtWild) ||
 126               ((mrlen == WILDCARD) && !bProtein && !bProtWild)));
 127
 128   if (debug) {
 129     fprintf(debug,"searching residue: %4s atom: %4s\n",resnm,atomnm);
 130     if (j == NOTFOUND)
 131       fprintf(debug," not successful\n");
 132     else
 133       fprintf(debug," match: %4s %4s\n",ap->resnm[j],ap->atomnm[j]);
 134   }
 135   return j;
 136 }
 137
 138 static void add_prop(aprop_t *ap,gmx_residuetype_t restype,
 139                      char *resnm,char *atomnm,
 140                      real p,int line)
 141 {
 142   int  i,j;
 143   gmx_bool bExact;
 144
 145   j = get_prop_index(ap,restype,resnm,atomnm,&bExact);
 146
 147   if (!bExact) {
 148     if (ap->nprop >= ap->maxprop) {
 149       ap->maxprop += 10;
 150       srenew(ap->resnm,ap->maxprop);
 151       srenew(ap->atomnm,ap->maxprop);
 152       srenew(ap->value,ap->maxprop);
 153       srenew(ap->bAvail,ap->maxprop);
 154       for(i=ap->nprop; (i<ap->maxprop); i++) {
 155         ap->atomnm[i] = NULL;
 156         ap->resnm[i]  = NULL;
 157         ap->value[i]  = 0;
 158         ap->bAvail[i] = FALSE;
 159       }
 160     }
 161     ap->atomnm[ap->nprop] = strdup(atomnm);
 162     ap->resnm[ap->nprop]  = strdup(resnm);
 163     j = ap->nprop;
 164     ap->nprop++;
 165   }
 166   if (ap->bAvail[j]) {
 167     if (ap->value[j] == p)
 168       fprintf(stderr,"Warning double identical entries for %s %s %g on line %d in file %s\n",
 169               resnm,atomnm,p,line,ap->db);
 170     else {
 171       fprintf(stderr,"Warning double different entries %s %s %g and %g on line %d in file %s\n"
 172               "Using last entry (%g)\n",
 173               resnm,atomnm,p,ap->value[j],line,ap->db,p);
 174       ap->value[j] = p;
 175     }
 176   }
 177   else {
 178     ap->bAvail[j] = TRUE;
 179     ap->value[j]  = p;
 180   }
 181 }
 182
 183 static void read_prop(gmx_atomprop_t aps,int eprop,double factor)
 184 {
 185   gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 186   FILE   *fp;
 187   char   line[STRLEN],resnm[32],atomnm[32];
 188   double pp;
 189   int    line_no;
 190   aprop_t *ap;
 191
 192   ap = &ap2->prop[eprop];
 193
 194   fp      = libopen(ap->db);
 195   line_no = 0;
 196   while(get_a_line(fp,line,STRLEN)) {
 197     line_no++;
 198     if (sscanf(line,"%s %s %lf",resnm,atomnm,&pp) == 3) {
 199       pp *= factor;
 200       add_prop(ap,aps->restype,resnm,atomnm,pp,line_no);
 201     }
 202     else
 203       fprintf(stderr,"WARNING: Error in file %s at line %d ignored\n",
 204               ap->db,line_no);
 205   }
 206
 207   /* for libraries we can use the low-level close routines */
 208   ffclose(fp);
 209
 210   ap->bSet = TRUE;
 211 }
 212
 213 static void set_prop(gmx_atomprop_t aps,int eprop)
 214 {
 215   gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 216   const char *fns[epropNR]  = { "atommass.dat", "vdwradii.dat", "dgsolv.dat", "electroneg.dat", "elements.dat" };
 217   double fac[epropNR] = { 1.0,    1.0,  418.4, 1.0, 1.0 };
 218   double def[epropNR] = { 12.011, 0.14, 0.0, 2.2, -1 };
 219   aprop_t *ap;
 220
 221   ap = &ap2->prop[eprop];
 222   if (!ap->bSet) {
 223     ap->db  = strdup(fns[eprop]);
 224     ap->def = def[eprop];
 225     read_prop(aps,eprop,fac[eprop]);
 226
 227     if (debug)
 228       fprintf(debug,"Entries in %s: %d\n",ap->db,ap->nprop);
 229
 230     if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW)) {
 231       printf("\n"
 232              "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
 233              "         based on residue and atom names, since they could not be\n"
 234              "         definitively assigned from the information in your input\n"
 235              "         files. These guessed numbers might deviate from the mass\n"
 236              "         and radius of the atom type. Please check the output\n"
 237              "         files if necessary.\n\n");
 238       aps->bWarned = TRUE;
 239     }
 240   }
 241 }
 242
 243 gmx_atomprop_t gmx_atomprop_init(void)
 244 {
 245   gmx_atomprop *aps;
 246   int p;
 247
 248   snew(aps,1);
 249
 250   gmx_residuetype_init(&aps->restype);
 251   aps->bWarned = FALSE;
 252
 253   return (gmx_atomprop_t)aps;
 254 }
 255
 256 static void destroy_prop(aprop_t *ap)
 257 {
 258   int i;
 259
 260   if (ap->bSet) {
 261     sfree(ap->db);
 262
 263     for(i=0; i<ap->nprop; i++) {
 264       sfree(ap->atomnm[i]);
 265       sfree(ap->resnm[i]);
 266     }
 267     sfree(ap->atomnm);
 268     sfree(ap->resnm);
 269     sfree(ap->bAvail);
 270     sfree(ap->value);
 271   }
 272 }
 273
 274 void gmx_atomprop_destroy(gmx_atomprop_t aps)
 275 {
 276   gmx_atomprop *ap = (gmx_atomprop*) aps;
 277   int p;
 278
 279   if (aps == NULL) {
 280     printf("\nWARNING: gmx_atomprop_destroy called with a NULL pointer\n\n");
 281     return;
 282   }
 283
 284   for(p=0; p<epropNR; p++) {
 285     destroy_prop(&ap->prop[p]);
 286   }
 287
 288   gmx_residuetype_destroy(ap->restype);
 289
 290   sfree(ap);
 291 }
 292
 293 gmx_bool gmx_atomprop_query(gmx_atomprop_t aps,
 294                         int eprop,const char *resnm,const char *atomnm,
 295                         real *value)
 296 {
 297   gmx_atomprop *ap = (gmx_atomprop*) aps;
 298   size_t i;
 299   int  j;
 300 #define MAXQ 32
 301   char atomname[MAXQ],resname[MAXQ];
 302   gmx_bool bExact;
 303
 304   set_prop(aps,eprop);
 305   if ((strlen(atomnm) > MAXQ-1) || (strlen(resnm) > MAXQ-1)) {
 306     if (debug)
 307       fprintf(debug,"WARNING: will only compare first %d characters\n",
 308               MAXQ-1);
 309   }
 310   if (isdigit(atomnm[0])) {
 311     /* put digit after atomname */
 312     for (i=1; (i<min(MAXQ-1,strlen(atomnm))); i++)
 313       atomname[i-1] = atomnm[i];
 314     atomname[i-1] = atomnm[0];
 315     atomname[i]   = '\0';
 316   }
 317   else {
 318     strncpy(atomname,atomnm,MAXQ-1);
 319   }
 320   strncpy(resname,resnm,MAXQ-1);
 321
 322   j = get_prop_index(&(ap->prop[eprop]),ap->restype,resname,
 323                      atomname,&bExact);
 324
 325   if (j >= 0) {
 326     *value = ap->prop[eprop].value[j];
 327     return TRUE;
 328   }
 329   else {
 330     *value = ap->prop[eprop].def;
 331     return FALSE;
 332   }
 333 }
 334
 335 char *gmx_atomprop_element(gmx_atomprop_t aps,int atomnumber)
 336 {
 337   gmx_atomprop *ap = (gmx_atomprop*) aps;
 338   int i;
 339
 340   set_prop(aps,epropElement);
 341   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 342     if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber) {
 343       return ap->prop[epropElement].atomnm[i];
 344     }
 345   }
 346   return NULL;
 347 }
 348
 349 int gmx_atomprop_atomnumber(gmx_atomprop_t aps,const char *elem)
 350 {
 351   gmx_atomprop *ap = (gmx_atomprop*) aps;
 352   int i;
 353
 354   set_prop(aps,epropElement);
 355   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 356     if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i],elem) == 0) {
 357       return gmx_nint(ap->prop[epropElement].value[i]);
 358     }
 359   }
 360   return NOTSET;
 361 }