src/gmxlib/atomprop.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <ctype.h>
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "string2.h"
  46 #include "futil.h"
  47 #include "maths.h"
  48 #include "gmx_fatal.h"
  49 #include "atomprop.h"
  50 #include "maths.h"
  51 #include "macros.h"
  52 #include "index.h"
  53 #include "strdb.h"
  54 #include "copyrite.h"
  55
  56 typedef struct {
  57   gmx_bool   bSet;
  58   int    nprop,maxprop;
  59   char   *db;
  60   double def;
  61   char   **atomnm;
  62   char   **resnm;
  63   gmx_bool   *bAvail;
  64   real   *value;
  65 } aprop_t;
  66
  67 typedef struct gmx_atomprop {
  68   gmx_bool       bWarned;
  69   aprop_t    prop[epropNR];
  70   gmx_residuetype_t restype;
  71 } gmx_atomprop;
  72
  73
  74
  75 /* NOTFOUND should be smallest, others larger in increasing priority */
  76 enum { NOTFOUND=-4, WILDCARD, WILDPROT, PROTEIN };
  77
  78 /* return number of matching characters,
  79    or NOTFOUND if not at least all characters in char *database match */
  80 static int dbcmp_len(char *search, char *database)
  81 {
  82   int i;
  83
  84   i=0;
  85   while(search[i] && database[i] && (search[i]==database[i]) )
  86     i++;
  87
  88   if (database[i])
  89     i=NOTFOUND;
  90   return i;
  91 }
  92
  93 static int get_prop_index(aprop_t *ap,gmx_residuetype_t restype,
  94                           char *resnm,char *atomnm,
  95                           gmx_bool *bExact)
  96 {
  97   int  i,j=NOTFOUND;
  98   long int alen,rlen;
  99   long int malen,mrlen;
 100   gmx_bool bProtein,bProtWild;
 101
 102   bProtein  = gmx_residuetype_is_protein(restype,resnm);
 103   bProtWild = (strcmp(resnm,"AAA")==0);
 104   malen = NOTFOUND;
 105   mrlen = NOTFOUND;
 106   for(i=0; (i<ap->nprop); i++) {
 107     rlen = dbcmp_len(resnm, ap->resnm[i]);
 108     if (rlen == NOTFOUND) {
 109       if ( (strcmp(ap->resnm[i],"*")==0) ||
 110            (strcmp(ap->resnm[i],"???")==0) )
 111         rlen=WILDCARD;
 112       else if (strcmp(ap->resnm[i],"AAA")==0)
 113         rlen=WILDPROT;
 114     }
 115     alen = dbcmp_len(atomnm, ap->atomnm[i]);
 116     if ( (alen > NOTFOUND) && (rlen > NOTFOUND)) {
 117       if ( ( (alen > malen) && (rlen >= mrlen)) ||
 118            ( (rlen > mrlen) && (alen >= malen) ) ) {
 119         malen = alen;
 120         mrlen = rlen;
 121         j     = i;
 122       }
 123     }
 124   }
 125
 126   *bExact = ((malen == (long int)strlen(atomnm)) &&
 127              ((mrlen == (long int)strlen(resnm)) ||
 128               ((mrlen == WILDPROT) && bProtWild) ||
 129               ((mrlen == WILDCARD) && !bProtein && !bProtWild)));
 130
 131   if (debug) {
 132     fprintf(debug,"searching residue: %4s atom: %4s\n",resnm,atomnm);
 133     if (j == NOTFOUND)
 134       fprintf(debug," not successful\n");
 135     else
 136       fprintf(debug," match: %4s %4s\n",ap->resnm[j],ap->atomnm[j]);
 137   }
 138   return j;
 139 }
 140
 141 static void add_prop(aprop_t *ap,gmx_residuetype_t restype,
 142                      char *resnm,char *atomnm,
 143                      real p,int line)
 144 {
 145   int  i,j;
 146   gmx_bool bExact;
 147
 148   j = get_prop_index(ap,restype,resnm,atomnm,&bExact);
 149
 150   if (!bExact) {
 151     if (ap->nprop >= ap->maxprop) {
 152       ap->maxprop += 10;
 153       srenew(ap->resnm,ap->maxprop);
 154       srenew(ap->atomnm,ap->maxprop);
 155       srenew(ap->value,ap->maxprop);
 156       srenew(ap->bAvail,ap->maxprop);
 157       for(i=ap->nprop; (i<ap->maxprop); i++) {
 158         ap->atomnm[i] = NULL;
 159         ap->resnm[i]  = NULL;
 160         ap->value[i]  = 0;
 161         ap->bAvail[i] = FALSE;
 162       }
 163     }
 164     upstring(atomnm);
 165     upstring(resnm);
 166     ap->atomnm[ap->nprop] = strdup(atomnm);
 167     ap->resnm[ap->nprop]  = strdup(resnm);
 168     j = ap->nprop;
 169     ap->nprop++;
 170   }
 171   if (ap->bAvail[j]) {
 172     if (ap->value[j] == p)
 173       fprintf(stderr,"Warning double identical entries for %s %s %g on line %d in file %s\n",
 174               resnm,atomnm,p,line,ap->db);
 175     else {
 176       fprintf(stderr,"Warning double different entries %s %s %g and %g on line %d in file %s\n"
 177               "Using last entry (%g)\n",
 178               resnm,atomnm,p,ap->value[j],line,ap->db,p);
 179       ap->value[j] = p;
 180     }
 181   }
 182   else {
 183     ap->bAvail[j] = TRUE;
 184     ap->value[j]  = p;
 185   }
 186 }
 187
 188 static void read_prop(gmx_atomprop_t aps,int eprop,double factor)
 189 {
 190   gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 191   FILE   *fp;
 192   char   line[STRLEN],resnm[32],atomnm[32];
 193   double pp;
 194   int    line_no;
 195   aprop_t *ap;
 196
 197   ap = &ap2->prop[eprop];
 198
 199   fp      = libopen(ap->db);
 200   line_no = 0;
 201   while(get_a_line(fp,line,STRLEN)) {
 202     line_no++;
 203     if (sscanf(line,"%s %s %lf",resnm,atomnm,&pp) == 3) {
 204       pp *= factor;
 205       add_prop(ap,aps->restype,resnm,atomnm,pp,line_no);
 206     }
 207     else
 208       fprintf(stderr,"WARNING: Error in file %s at line %d ignored\n",
 209               ap->db,line_no);
 210   }
 211
 212   /* for libraries we can use the low-level close routines */
 213   ffclose(fp);
 214
 215   ap->bSet = TRUE;
 216 }
 217
 218 static void set_prop(gmx_atomprop_t aps,int eprop)
 219 {
 220   gmx_atomprop *ap2 = (gmx_atomprop*) aps;
 221   const char *fns[epropNR]  = { "atommass.dat", "vdwradii.dat", "dgsolv.dat", "electroneg.dat", "elements.dat" };
 222   double fac[epropNR] = { 1.0,    1.0,  418.4, 1.0, 1.0 };
 223   double def[epropNR] = { 12.011, 0.14, 0.0, 2.2, -1 };
 224   aprop_t *ap;
 225
 226   ap = &ap2->prop[eprop];
 227   if (!ap->bSet) {
 228     ap->db  = strdup(fns[eprop]);
 229     ap->def = def[eprop];
 230     read_prop(aps,eprop,fac[eprop]);
 231
 232     if (debug)
 233       fprintf(debug,"Entries in %s: %d\n",ap->db,ap->nprop);
 234
 235     if ( ( (!aps->bWarned) && (eprop == epropMass) ) || (eprop == epropVDW)) {
 236       printf("\n"
 237              "WARNING: Masses and atomic (Van der Waals) radii will be guessed\n"
 238              "         based on residue and atom names, since they could not be\n"
 239              "         definitively assigned from the information in your input\n"
 240              "         files. These guessed numbers might deviate from the mass\n"
 241              "         and radius of the atom type. Please check the output\n"
 242              "         files if necessary.\n\n");
 243       aps->bWarned = TRUE;
 244     }
 245   }
 246 }
 247
 248 gmx_atomprop_t gmx_atomprop_init(void)
 249 {
 250   gmx_atomprop *aps;
 251   int p;
 252
 253   snew(aps,1);
 254
 255   gmx_residuetype_init(&aps->restype);
 256   aps->bWarned = FALSE;
 257
 258   return (gmx_atomprop_t)aps;
 259 }
 260
 261 static void destroy_prop(aprop_t *ap)
 262 {
 263   int i;
 264
 265   if (ap->bSet) {
 266     sfree(ap->db);
 267
 268     for(i=0; i<ap->nprop; i++) {
 269       sfree(ap->atomnm[i]);
 270       sfree(ap->resnm[i]);
 271     }
 272     sfree(ap->atomnm);
 273     sfree(ap->resnm);
 274     sfree(ap->bAvail);
 275     sfree(ap->value);
 276   }
 277 }
 278
 279 void gmx_atomprop_destroy(gmx_atomprop_t aps)
 280 {
 281   gmx_atomprop *ap = (gmx_atomprop*) aps;
 282   int p;
 283
 284   if (aps == NULL) {
 285     printf("\nWARNING: gmx_atomprop_destroy called with a NULL pointer\n\n");
 286     return;
 287   }
 288
 289   for(p=0; p<epropNR; p++) {
 290     destroy_prop(&ap->prop[p]);
 291   }
 292
 293   gmx_residuetype_destroy(ap->restype);
 294
 295   sfree(ap);
 296 }
 297
 298 gmx_bool gmx_atomprop_query(gmx_atomprop_t aps,
 299                         int eprop,const char *resnm,const char *atomnm,
 300                         real *value)
 301 {
 302   gmx_atomprop *ap = (gmx_atomprop*) aps;
 303   size_t i;
 304   int  j;
 305 #define MAXQ 32
 306   char atomname[MAXQ],resname[MAXQ];
 307   gmx_bool bExact;
 308
 309   set_prop(aps,eprop);
 310   if ((strlen(atomnm) > MAXQ-1) || (strlen(resnm) > MAXQ-1)) {
 311     if (debug)
 312       fprintf(debug,"WARNING: will only compare first %d characters\n",
 313               MAXQ-1);
 314   }
 315   if (isdigit(atomnm[0])) {
 316     /* put digit after atomname */
 317     for (i=1; (i<min(MAXQ-1,strlen(atomnm))); i++)
 318       atomname[i-1] = atomnm[i];
 319     atomname[i-1] = atomnm[0];
 320     atomname[i]   = '\0';
 321   }
 322   else {
 323     strncpy(atomname,atomnm,MAXQ-1);
 324   }
 325   upstring(atomname);
 326   strncpy(resname,resnm,MAXQ-1);
 327   upstring(resname);
 328
 329   j = get_prop_index(&(ap->prop[eprop]),ap->restype,resname,
 330                      atomname,&bExact);
 331
 332   if (j >= 0) {
 333     *value = ap->prop[eprop].value[j];
 334     return TRUE;
 335   }
 336   else {
 337     *value = ap->prop[eprop].def;
 338     return FALSE;
 339   }
 340 }
 341
 342 char *gmx_atomprop_element(gmx_atomprop_t aps,int atomnumber)
 343 {
 344   gmx_atomprop *ap = (gmx_atomprop*) aps;
 345   int i;
 346
 347   set_prop(aps,epropElement);
 348   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 349     if (gmx_nint(ap->prop[epropElement].value[i]) == atomnumber) {
 350       return ap->prop[epropElement].atomnm[i];
 351     }
 352   }
 353   return NULL;
 354 }
 355
 356 int gmx_atomprop_atomnumber(gmx_atomprop_t aps,const char *elem)
 357 {
 358   gmx_atomprop *ap = (gmx_atomprop*) aps;
 359   int i;
 360
 361   set_prop(aps,epropElement);
 362   for(i=0; (i<ap->prop[epropElement].nprop); i++) {
 363     if (gmx_strcasecmp(ap->prop[epropElement].atomnm[i],elem) == 0) {
 364       return gmx_nint(ap->prop[epropElement].value[i]);
 365     }
 366   }
 367   return NOTSET;
 368 }