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     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              "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
 235              "         based on residue and atom names, since they could not be\n"
 236              "         definitively assigned from the information in your input\n"
 237              "         files. These guessed numbers might deviate from the mass\n"
 238              "         and radius of the atom type. Please check the output\n"
 239              "         files if necessary.\n\n");
 240       aps->bWarned = TRUE;
 241     }
 242   }
 243 }
 244
 245 gmx_atomprop_t gmx_atomprop_init(void)
 246 {
 247   gmx_atomprop *aps;
 248   int p;
 249
 250   snew(aps,1);
 251
 252   gmx_residuetype_init(&aps->restype);
 253   aps->bWarned = FALSE;
 254
 255   return (gmx_atomprop_t)aps;
 256 }
 257
 258 static void destroy_prop(aprop_t *ap)
 259 {
 260   int i;
 261
 262   if (ap->bSet) {
 263     sfree(ap->db);
 264
 265     for(i=0; i<ap->nprop; i++) {
 266       sfree(ap->atomnm[i]);
 267       sfree(ap->resnm[i]);
 268     }
 269     sfree(ap->atomnm);
 270     sfree(ap->resnm);
 271     sfree(ap->bAvail);
 272     sfree(ap->value);
 273   }
 274 }
 275
 276 void gmx_atomprop_destroy(gmx_atomprop_t aps)
 277 {
 278   gmx_atomprop *ap = (gmx_atomprop*) aps;
 279   int p;
 280
 281   if (aps == NULL) {
 282     printf("\nWARNING: gmx_atomprop_destroy called with a NULL pointer\n\n");
 283     return;
 284   }
 285
 286   for(p=0; p<epropNR; p++) {
 287     destroy_prop(&ap->prop[p]);
 288   }
 289
 290   gmx_residuetype_destroy(ap->restype);
 291
 292   sfree(ap);
 293 }
 294
 295 gmx_bool gmx_atomprop_query(gmx_atomprop_t aps,
 296                         int eprop,const char *resnm,const char *atomnm,
 297                         real *value)
 298 {
 299   gmx_atomprop *ap = (gmx_atomprop*) aps;
 300   size_t i;
 301   int  j;
 302 #define MAXQ 32
 303   char atomname[MAXQ],resname[MAXQ];
 304   gmx_bool bExact;
 305
 306   set_prop(aps,eprop);
 307   if ((strlen(atomnm) > MAXQ-1) || (strlen(resnm) > MAXQ-1)) {
 308     if (debug)
 309       fprintf(debug,"WARNING: will only compare first %d characters\n",
 310               MAXQ-1);
 311   }
 312   if (isdigit(atomnm[0])) {
 313     /* put digit after atomname */
 314     for (i=1; (i<min(MAXQ-1,strlen(atomnm))); i++)
 315       atomname[i-1] = atomnm[i];
 316     atomname[i-1] = atomnm[0];
 317     atomname[i]   = '\0';
 318   }
 319   else {
 320     strncpy(atomname,atomnm,MAXQ-1);
 321   }
 322   upstring(atomname);
 323   strncpy(resname,resnm,MAXQ-1);
 324   upstring(resname);
 325
 326   j = get_prop_index(&(ap->prop[eprop]),ap->restype,resname,
 327                      atomname,&bExact);
 328
 329   if (j >= 0) {
 330     *value = ap->prop[eprop].value[j];
 331     return TRUE;
 332   }
 333   else {
 334     *value = ap->prop[eprop].def;
 335     return FALSE;
 336   }
 337 }
 338
 339 char *gmx_atomprop_element(gmx_atomprop_t aps,int atomnumber)
 340 {
 341   gmx_atomprop *ap = (gmx_atomprop*) aps;
 342   int i;
 343
 344   set_prop(aps,epropElement);
 345   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 346     if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber) {
 347       return ap->prop[epropElement].atomnm[i];
 348     }
 349   }
 350   return NULL;
 351 }
 352
 353 int gmx_atomprop_atomnumber(gmx_atomprop_t aps,const char *elem)
 354 {
 355   gmx_atomprop *ap = (gmx_atomprop*) aps;
 356   int i;
 357
 358   set_prop(aps,epropElement);
 359   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 360     if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i],elem) == 0) {
 361       return gmx_nint(ap->prop[epropElement].value[i]);
 362     }
 363   }
 364   return NOTSET;
 365 }