src/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   bool   bSet;
  55   int    nprop,maxprop;
  56   char   *db;
  57   double def;
  58   char   **atomnm;
  59   char   **resnm;
  60   bool   *bAvail;
  61   real   *value;
  62 } aprop_t;
  63
  64 typedef struct gmx_atomprop {
  65   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                           bool *bExact)
  93 {
  94   int  i,j=NOTFOUND;
  95   long int alen,rlen;
  96   long int malen,mrlen;
  97   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   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     upstring(atomnm);
 162     upstring(resnm);
 163     ap->atomnm[ap->nprop] = strdup(atomnm);
 164     ap->resnm[ap->nprop]  = strdup(resnm);
 165     j = ap->nprop;
 166     ap->nprop++;
 167   }
 168   if (ap->bAvail[j]) {
 169     if (ap->value[j] == p)
 170       fprintf(stderr,"Warning double identical entries for %s %s %g on line %d in file %s\n",
 171               resnm,atomnm,p,line,ap->db);
 172     else {
 173       fprintf(stderr,"Warning double different entries %s %s %g and %g on line %d in file %s\n"
 174               "Using last entry (%g)\n",
 175               resnm,atomnm,p,ap->value[j],line,ap->db,p);
 176       ap->value[j] = p;
 177     }
 178   }
 179   else {
 180     ap->bAvail[j] = TRUE;
 181     ap->value[j]  = p;
 182   }
 183 }
 184
 185 static void read_prop(gmx_atomprop_t aps,int eprop,double factor)
 186 {
 187   gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 188   FILE   *fp;
 189   char   line[STRLEN],resnm[32],atomnm[32];
 190   double pp;
 191   int    line_no;
 192   aprop_t *ap;
 193
 194   ap = &ap2->prop[eprop];
 195
 196   fp      = libopen(ap->db);
 197   line_no = 0;
 198   while(get_a_line(fp,line,STRLEN)) {
 199     line_no++;
 200     if (sscanf(line,"%s %s %lf",resnm,atomnm,&pp) == 3) {
 201       pp *= factor;
 202       add_prop(ap,aps->restype,resnm,atomnm,pp,line_no);
 203     }
 204     else
 205       fprintf(stderr,"WARNING: Error in file %s at line %d ignored\n",
 206               ap->db,line_no);
 207   }
 208
 209   /* for libraries we can use the low-level close routines */
 210   ffclose(fp);
 211
 212   ap->bSet = TRUE;
 213 }
 214
 215 static void set_prop(gmx_atomprop_t aps,int eprop)
 216 {
 217   gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 218   const char *fns[epropNR]  = { "atommass.dat", "vdwradii.dat", "dgsolv.dat", "electroneg.dat", "elements.dat" };
 219   double fac[epropNR] = { 1.0,    1.0,  418.4, 1.0, 1.0 };
 220   double def[epropNR] = { 12.011, 0.14, 0.0, 2.2, -1 };
 221   aprop_t *ap;
 222
 223   ap = &ap2->prop[eprop];
 224   if (!ap->bSet) {
 225     ap->db  = strdup(fns[eprop]);
 226     ap->def = def[eprop];
 227     read_prop(aps,eprop,fac[eprop]);
 228
 229     if (debug)
 230       fprintf(debug,"Entries in %s: %d\n",ap->db,ap->nprop);
 231
 232     if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW)) {
 233       printf("\n");
 234       printf("WARNING: masses and atomic (Van der Waals) radii will be determined\n");
 235       printf("         based on residue and atom names. These numbers can deviate\n");
 236       printf("         from the correct mass and radius of the atom type.\n");
 237       printf("\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 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   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++] = atomnm[0];
 315     atomname[i]   = '\0';
 316   }
 317   else {
 318     strncpy(atomname,atomnm,MAXQ-1);
 319   }
 320   upstring(atomname);
 321   strncpy(resname,resnm,MAXQ-1);
 322   upstring(resname);
 323
 324   j = get_prop_index(&(ap->prop[eprop]),ap->restype,resname,
 325                      atomname,&bExact);
 326
 327   if (j >= 0) {
 328     *value = ap->prop[eprop].value[j];
 329     return TRUE;
 330   }
 331   else {
 332     *value = ap->prop[eprop].def;
 333     return FALSE;
 334   }
 335 }
 336
 337 char *gmx_atomprop_element(gmx_atomprop_t aps,int atomnumber)
 338 {
 339   gmx_atomprop *ap = (gmx_atomprop*) aps;
 340   int i;
 341
 342   set_prop(aps,epropElement);
 343   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 344     if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber) {
 345       return ap->prop[epropElement].atomnm[i];
 346     }
 347   }
 348   return NULL;
 349 }
 350
 351 int gmx_atomprop_atomnumber(gmx_atomprop_t aps,const char *elem)
 352 {
 353   gmx_atomprop *ap = (gmx_atomprop*) aps;
 354   int i;
 355
 356   set_prop(aps,epropElement);
 357   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 358     if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i],elem) == 0) {
 359       return gmx_nint(ap->prop[epropElement].value[i]);
 360     }
 361   }
 362   return NOTSET;
 363 }